Fluid-mediated sources of granular temperature at finite Reynolds numbers

Abstract

We derive analytical solutions for hydrodynamic sources and sinks to granular temperature in moderately dense suspensions of elastic particles at finite Reynolds numbers. Modelling the neighbour-induced drag disturbances with a Langevin equation allows an exact solution for the joint fluctuating acceleration–velocity distribution function $P(v^{\prime },a^{\prime };t)$. Quadrant-conditioned covariance integrals of $P(v^{\prime },a^{\prime };t)$ yield the hydrodynamic source and sink that dictate the evolution of granular temperature that can be used in Eulerian two-fluid models. Analytical predictions agree with benchmark data from particle-resolved direct numerical simulations and show promise as a general theory from gas–solid to bubbly flows.

1. Introduction

Hydrodynamic interactions between viscous fluids and a disperse phase (particles, droplets or bubbles) give rise to complex dynamics that are important to many engineering and environmental applications. From the production of biofuels to post-combustion carbon capture, multiphase reactors are at the heart of nearly all chemical transformation processes. Environmental systems such as gravity currents, debris flows and sand dunes are also of great societal interest. Broadly speaking, the aforementioned examples are characterized by turbulent fluid flow and moderate to high solids volume fractions.

The kinetic theory of rapid granular flows is now generally accepted as a valid description of moderately dense particulate systems. Namely, Enskog theory connects rigorously the evolution of the single-particle distribution function

(1.1)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{x}} f + \boldsymbol{\nabla}_{\boldsymbol{v}} \boldsymbol{\cdot} \left( \boldsymbol{a}_h f \right) = \boldsymbol{\varOmega} \end{equation}

to continuum equations for the solids-phase moments, e.g. mass, momentum and kinetic energy (granular temperature) (Lun et al. 1984; Garzó & Dufty 1999). In (1.1), $f$ is the average number density of particles at location $\boldsymbol {x}$ with velocity $\boldsymbol {v}$ at time $t$, subscripts on the $\boldsymbol {\nabla }$ operator denote the sample space variables being differentiated, $\boldsymbol {a}_h$ is the acceleration due to hydrodynamic interactions (other external forces are neglected), and $\boldsymbol {\varOmega }$ is the collision operator. For granular flows, one would ignore the term involving $\boldsymbol {a}_h$ in (1.1) and complete a Chapman–Enskog expansion to arrive at macroscopic balance relations (e.g. granular temperature) and transport coefficients. However, neglecting the presence of a viscous fluid makes granular kinetic theories incapable of predicting solids motion in fluidized systems.

To address this deficiency, driven systems with stochastic velocity fluctuations have been theoretically (Puglisi et al. 1998; Cafiero, Luding & Jürgen Herrmann 2000; Pagonabarraga et al. 2001; Srebro & Levine 2004) and experimentally (Yu, Schröter & Sperl 2020) examined and incorporated into Chapman–Enskog expansions (Garzó et al. 2012; Khalil & Garzó 2020). Additionally, closures for hydrodynamic sources and sinks to granular temperature have been proposed from phenomenological scaling arguments and multipole simulations (Koch & Sangani 1999). However, fluid-mediated sources of granular temperature, resulting from a statistical description of hydrodynamic forces, have not been validated at finite Reynolds numbers and solids volume fraction. We show, for the first time, consistency between analytical solutions for sources and sinks to granular temperature and data obtained from particle-resolved direct numerical simulation (PR–DNS) at finite Reynolds numbers and solids volume fraction.

2. Homogeneous fluidization of spherical particles

In this study, we consider homogeneous fluidization of elastic smooth spheres (see figure 1). Elastic contacts conserve kinetic energy, thus collisions may redistribute particle velocity fluctuations $\boldsymbol {v}_p^\prime$ without providing a direct sink to granular temperature, defined as $T=\langle \boldsymbol {v}_{p}^{\prime } \boldsymbol {\cdot } \boldsymbol {v}_{p}^{\prime }\rangle /3$. Therefore, capturing granular temperature in the present system requires appropriate specification of the particle acceleration due to hydrodynamic interactions. Here, angled brackets denote an ensemble average, and a single prime denotes a fluctuation from the ensemble average. A constant mean pressure gradient is imposed on the fluid to obtain a uniform slip velocity with mean Reynolds number ${Re}_m = (1-\langle \phi \rangle )\rho _fd_p| \langle \boldsymbol {w} \rangle |/\mu _f$, which drives the particles to a steady granular temperature, where $\langle \phi \rangle$ is the mean solids volume fraction, $\rho _f$ is the fluid density, $d_p$ is the particle diameter, $\langle \boldsymbol {w} \rangle = \langle \boldsymbol {u}_f \rangle - \langle \boldsymbol {v}_p \rangle$ is the mean slip velocity between the fluid $\langle \boldsymbol {u}_f \rangle$ and particles $\langle \boldsymbol {v}_p \rangle$, and $\mu _f$ is the dynamic viscosity.

Figure 1. PR–DNS of homogeneous fluidization at ${Re}_m=20$, $\rho _p/\rho _f=1000$ and $\langle \phi \rangle =0.1$. Arrows denote fluid streamlines; velocity magnitude shown in colour.

The evolution of granular temperature in this system is dictated by the covariance of fluctuating hydrodynamic acceleration $\boldsymbol {a}_{h}^{\prime }$ and particle velocity $\boldsymbol {v}_{p}^{\prime }$:

(2.1)\begin{equation} \frac{{\rm d}T}{{\rm d}t} \equiv S-\varGamma =\frac{2}{3}\,\langle \boldsymbol{a}_{h}^{\prime}\boldsymbol{\cdot}\boldsymbol{v}_{p}^{\prime}\rangle, \end{equation}

which contains sources $S$ and sinks $\varGamma$ (Tenneti, Mehrabadi & Subramaniam 2016). Therefore, accurate descriptions for $\boldsymbol {a}_{h}^{\prime }$ are crucial to capturing $T(t)$ in homogeneous fluidization of elastic spheres.

For a particle subjected to hydrodynamic forces, its total translational velocity $\boldsymbol {v}_p = \langle \boldsymbol {v}_p \rangle + \boldsymbol {v}_p^{\prime }$ follows from

(2.2)\begin{equation} m_p\,\frac{{\rm d}\boldsymbol{v}_p}{{\rm d}t}=\int_{\partial \varOmega} \boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S, \end{equation}

where $m_p$ is the particle mass, and surface integration of the fluid stress tensor $\boldsymbol {\tau }$ (comprised of pressure and viscous stress) gives the total hydrodynamic acceleration $\boldsymbol {a}_h$. A model for $\boldsymbol {a}_h(t)$ has been obtained in the limit of an isolated sphere and Stokes flow by Maxey & Riley (1983); it involves a superposition of forces from the undisturbed fluid flow, quasi-steady drag, added mass and Basset history. At finite solids volume fractions and Reynolds numbers, analytical evaluation of the fluid stress integral is not tractable. Alternatively, correlations are often obtained from PR–DNS by ensemble averaging the net hydrodynamic force acting on a suspension. However, as shown in figure 1, particles interact with fluid wakes generated by their neighbours – referred to as pseudo-turbulent kinetic energy (PTKE) – leading to a distribution of hydrodynamic forces (Akiki, Jackson & Balachandar 2016) that drive relative motion between particles. The application of a drag correlation, obtained from ensemble averaging, cannot capture the distribution in hydrodynamic force. We describe this PTKE-induced drag distribution as a stochastic process in § 3. The predictive capability of the theory is emphasized in § 4 for $\{{Re}_m \in [10, 100 ]; \langle \phi \rangle \in [0.1, 0.4 ]; \rho _p/\rho _f \in [100, 1000 ] \}$, where $\rho _p$ is the particle density, and considerations are given to low density ratio (bubbly) flows $\rho _p/\rho _f \ll 1$ in § 4.3.

3. Stochastic process and solutions

At the conditions considered here, hydrodynamic disturbances are attributed to PTKE generated by fluid boundary layers of neighbouring particles (Mehrabadi et al. 2015; Shallcross, Fox & Capecelatro 2020). Thus the drag force is decomposed into a mean contribution and a stochastic fluctuation $\boldsymbol {a}_h^{\prime \prime }$ according to

(3.1)\begin{equation} \boldsymbol{a}_h\equiv\frac{1}{m_p}\int_{\partial \varOmega} \boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n}\,{\rm d}S = \frac{1}{\tau_d} \left(\langle \boldsymbol{u}_f \rangle - \boldsymbol{v}_p \right) + \boldsymbol{a}_h^{\prime\prime}, \end{equation}

where $\tau _d =\tau _p/ F$ is the hydrodynamic time scale, $\tau _p=\rho _p d_p^2/(18 \mu _f)$ is the Stokes time scale, and $F(\langle \phi \rangle,{Re}_m)$ is a correction to Stokes drag that may be obtained from PR–DNS. Ensemble averaging (3.1) gives $\langle \boldsymbol {a}_h \rangle$. Removal of $\langle \boldsymbol {a}_h \rangle$ yields a reference frame that moves with the mean particle velocity and a fluctuating acceleration $\boldsymbol {a}_h^{\prime } = -\boldsymbol {v}_p^{\prime }/\tau _d + \boldsymbol {a}_h^{\prime \prime }$, where it is reiterated that single primes denote a fluctuation from the ensemble average, while double primes denote a stochastic fluctuation. We describe $\boldsymbol {a}_h^{\prime \prime }$ with an appropriate acceleration Langevin equation (Lattanzi et al. 2021):

(3.2a)\begin{gather} \mathrm{d}\boldsymbol{v}_p^{\prime}={-}\frac{1}{\tau_d}\,\boldsymbol{v}_p^{\prime} \, \mathrm{d}t + \boldsymbol{a}_{h}^{\prime\prime} \, \mathrm{d}t, \end{gather}

(3.2b)\begin{gather}\mathrm{d}\boldsymbol{a}_{h}^{\prime\prime} ={-}\frac{1}{\tau_{a^{\prime\prime}}}\, \boldsymbol{a}_{h}^{\prime\prime} \, \mathrm{d}t + \sqrt{\frac{2}{\tau_{a^{\prime\prime}}}}\,\sigma_{a^{\prime\prime}} \,\mathrm{d}\boldsymbol{W}_t, \end{gather}

where $\tau _{a^{\prime \prime }}$ is the integral time scale of the stochastic acceleration, $\sigma _{a^{\prime \prime }}$ is the standard deviation, and $\mathrm {d}\boldsymbol {W}_t$ is a Wiener process increment.

At this point, we consider some subtleties associated with (3.1). When applying drag correlations to PR–DNS data, Tenneti et al. (2010) showed that incorrect dynamics for the granular temperature were obtained if the instantaneous particle velocity was utilized. In the absence of the stochastic fluctuation, $\boldsymbol {a}_h^{\prime }$ agrees with the observations of Tenneti et al. (2010) since velocity fluctuations can be dissipated only via $-\boldsymbol {v}_p^{\prime }/\tau _d$. However, Lattanzi et al. (2020) showed that $\boldsymbol {a}_h^{\prime \prime }$ leads to a covariance $\langle \boldsymbol {a}_{h}^{\prime \prime }\boldsymbol {\cdot }\boldsymbol {v}_{p}^{\prime }\rangle$ that acts as a net source to granular temperature when averaged over all of $\boldsymbol {a}_{h}^{\prime \prime }$–$\boldsymbol {v}_{p}^{\prime }$ phase space. Therefore, the fluctuating hydrodynamic acceleration $\boldsymbol {a}_h^{\prime }$ is comprised of a solely dissipative component $-\boldsymbol {v}_p^{\prime }/\tau _d$ as well as a contribution from the stochastic process $\boldsymbol {a}_h^{\prime \prime }$, which may act as a source or sink at the particle level but is overall a source when the entire suspension is considered.

Tenneti et al. (2010) showed that sources occur when the fluctuating acceleration is aligned with the fluctuating particle velocity (quadrants 1 and 3 in $a_h^{\prime }$–$v_p^{\prime }$ phase space), and sinks occur for the converse case (quadrants 2 and 4). To quantify $S$ and $\varGamma$, one must evaluate the quadrant-conditioned acceleration–velocity covariance; see (2.1). Thus we seek the joint probability distribution resulting from (3.2a) and (3.2b). Considering isotropic fluctuations and constant coefficients, we derive the acceleration–velocity distribution in Appendix A. Here, we report the salient result that the probability $P(v^{\prime },a^{\prime };t)=\mathcal {N}(0, \bar {\boldsymbol {\varSigma }}^{-1}(t) )$ is a normal distribution with time-dependent covariance tensor $\bar {\boldsymbol {\varSigma }}(t)$ given by

(3.3)\begin{equation} \left.\begin{gathered} \bar{\boldsymbol{\varSigma}}(t) = \left[\begin{array}{cc} \dfrac{1}{\tau_d^2}\,\sigma_{v^{\prime}}^2(t) - \dfrac{1}{\tau_d}\,\sigma_{v^{\prime}a^{\prime\prime}}(t) + \sigma_{a^{\prime\prime}}^2 & -\dfrac{1}{\tau_d}\,\sigma_{v^{\prime}}^2(t) + \sigma_{v^{\prime}a^{\prime\prime}}(t) \\ -\dfrac{1}{\tau_d}\,\sigma_{v^{\prime}}^2(t) + \sigma_{v^{\prime}a^{\prime\prime}}(t) & \sigma_{v^{\prime}}^2(t) \end{array} \right], \\ \sigma_{v^{\prime}}^2(t) = \sigma_{a^{\prime\prime}}^2 \hat{\tau}^{+} \left[ \tau_d \mathbb{E}_1 + 2 \hat{\tau}^{-} \left( \mathbb{E}_1 - \mathbb{E}_2 \right) + \mathbb{C}_0 \tau_d \left(1 - \mathbb{E}_1 \right) - 2 \rho_0\,\sqrt{\frac{\mathbb{C}_0 \tau_d}{\hat{\tau}^{+}}}\, \hat{\tau}^{-} \left( \mathbb{E}_1 - \mathbb{E}_2 \right) \right], \\ \sigma_{v^{\prime}a^{\prime\prime}}(t) = \sigma_{a^{\prime\prime}}^2 \hat{\tau}^{+} \left[ \mathbb{E}_2 + \rho_0\, \sqrt{\frac{\mathbb{C}_0 \tau_d}{\hat{\tau}^{+}}} \left( 1 - \mathbb{E}_2 \right)\right], \end{gathered}\right\} \end{equation}

where $\mathbb {E}_1 = [ 1 - \exp ( -2t/\tau _d ) ]$, $\mathbb {E}_2 = [ 1 - \exp ( -t/\hat {\tau }^{+} ) ]$, $\hat {\tau }^{+} = \tau _d \tau _{a^{\prime \prime }} /(\tau _d + \tau _{a^{\prime \prime }})$, $\hat {\tau }^{-} = \tau _d \tau _{a^{\prime \prime }} /(\tau _d - \tau _{a^{\prime \prime }})$, $\rho _0 = \sigma _{v^{\prime }a^{\prime \prime },0}/(\sigma _{v^{\prime },0} \, \sigma _{a^{\prime \prime }})$ is the initial correlation coefficient between the fluctuating velocity and stochastic acceleration, and $\mathbb {C}_{0} \geqslant 0$ is a proportionality constant that specifies the initial velocity variance $\sigma _{v^{\prime },0}^2$ as a fraction of the steady velocity variance $\sigma _{v^{\prime },\infty }^2 = \sigma _{a^{\prime \prime }}^2 \hat {\tau }^{+} \tau _d$. Following Stuart & Ord (2010), we derive relations for the quadrant-covariance of a joint-normal (see Appendix B):

(3.4a)\begin{gather} S(t) = \frac{2\sqrt{\varSigma_{11} \varSigma_{22}}}{\rm \pi} \left(\rho \arcsin{\rho} + \sqrt{1-\rho^2} + \frac{\rm \pi}{2} \rho \right), \end{gather}

(3.4b)\begin{gather}\varGamma(t) = \frac{2\sqrt{\varSigma_{11} \varSigma_{22}}}{\rm \pi} \left(\rho \arcsin{\rho} + \sqrt{1-\rho^2} - \frac{\rm \pi}{2} \rho \right), \end{gather}

where $\rho (t)=\varSigma _{12}/\sqrt {\varSigma _{11} \varSigma _{22}}$ is the correlation coefficient evaluated at time $t$.

Evaluation of (3.4) requires closure for the stochastic process. We approximate $\tau _{a^{\prime \prime }}$ with the mean-free time between successive collisions (Chapman, Cowling & Burnett 1970):

(3.5)\begin{equation} \tau_{a^{\prime\prime}} \approx \tau_{col} = \frac{d_p}{24 \langle \phi \rangle g_0}\,\sqrt{\frac{\rm \pi}{T}}, \end{equation}

where $g_0$ is the radial distribution function at contact (Ma & Ahmadi 1988). Details regarding motivation for $\tau _{a^{\prime \prime }} \approx \tau _{col}$ are provided in Lattanzi et al. (2021). Here, we summarize that the aforementioned study considers the time scale analysis of Wylie, Koch & Ladd (2003) and Mehrabadi et al. (2015), which characterizes a ratio between the mean particle response time to hydrodynamic force fluctuations and the mean time between collisions. For the hydrodynamic conditions emphasized in this study, results presented in Lattanzi et al. (2021) suggest that the mean collisional time scale is smaller than the mean hydrodynamic time scale, i.e. the flow is collisionally dominated. Furthermore, the steady granular temperature solution derived here, and given in (4.7), leads to quantitative agreement between theory and PR–DNS data for the hydrodynamic time scale (Lattanzi et al. 2021).

The drag time scale $\tau _d({Re}_m,\langle \phi \rangle )$ in (3.2a) is given by Tenneti, Garg & Subramaniam (2011), while the standard deviation $\sigma _{a^{\prime \prime }}$ in (3.2b) is closed with a correlation obtained from PR–DNS of freely evolving particle assemblies, given by

(3.6a)\begin{gather} \sigma_{a^{\prime\prime}} = \beta\,f_{\langle \phi \rangle}^{\sigma_F}\,f_{iso}\, \frac{(1-\langle \phi \rangle)\left| \langle \boldsymbol{w} \rangle \right|}{\tau_p}, \end{gather}

(3.6b)\begin{gather}f_{\langle \phi \rangle}^{\sigma_F} = 5.39\langle \phi \rangle - 4.00 \langle \phi \rangle^2 + 24.93 \langle \phi \rangle^3, \end{gather}

(3.6c)\begin{gather}f_{iso} = \left( 1 + 0.15\,{Re}_m^{0.687} \right). \end{gather}

Here, $f_{iso}$ represents the drag correction for an isolated particle using the classical correlation of Schiller & Naumann (1933). We note that the correlation provided in (3.6) differs from that of stationary particle assemblies reported in Lattanzi et al. (2021). Here, the $\beta$ factor is introduced to account for anisotropy. It is reiterated that stochastic fluctuations in the present theory are treated as isotropic, and the role of collisions is to redistribute components of the particle velocity fluctuations. To quantify hydrodynamic sources and sinks to granular temperature with an isotropic model, we average the force variance in all three directions. We note that $\beta$ does not drive anisotropy in the granular temperature since the model is isotropic; future work along the lines discussed in § 5 will be needed to capture such phenomena. For conditions considered here, the force variance extracted from PR–DNS is observed to be ${\sim }3$ times larger in the streamwise direction than the transverse directions (Lattanzi et al. 2021), yielding $\beta = \sqrt {5/9}$ in the present work. We note that force anisotropy ($\beta < 1$) may depend upon ${Re}_m$ or $\langle \phi \rangle$, and future efforts to parametrize this quantity would be beneficial.

The present correlation for drag variation is fit to PR–DNS with freely evolving particles, while the correlation provided in Lattanzi et al. (2021) is fit to PR–DNS with fixed particles. It has been shown that the variance in hydrodynamic force grows with ${Re}_T = \rho _f d_p \sqrt {T}/\mu _f$ (Huang et al. 2017). To gauge the role of particle mobility on drag variation, we consider the normalized standard deviation in streamwise hydrodynamic force for freely evolving and fixed particles (see figure 2). For $\langle \phi \rangle \geqslant 0.2$, the standard deviation in drag force is under-predicted by PR–DNS with fixed particles, when compared to freely evolving suspensions. Thus (3.6) takes the same form as that reported by Lattanzi et al. (2021), but with an anisotropy factor $\beta$ and a solids volume fraction correction $f_{\phi }^{\sigma _F}$ that is fit to data from PR–DNS with freely evolving particles. To this end, quantifying $\sigma _{a^{\prime \prime }}$ over a wide range of conditions would also be of great value to the present theory.

Figure 2. Standard deviation in drag force normalized by the drag force on an isolated sphere, $F_{single} = 3 {\rm \pi}\mu _f d_p f_{iso} (1-\langle \phi \rangle )| \langle \boldsymbol {w} \rangle |$ (Schiller & Naumann 1933). Symbols denote PR–DNS data with freely evolving (square) and fixed (circle) particles. The solid line denotes (3.6b), while the dashed line denotes $f_{\phi }^{\sigma _F}$ in Lattanzi et al. (2021); ${Re}_m= 10$ (red), ${Re}_m=20$ (blue), ${Re}_m=50$ (green), ${Re}_m=100$ (magenta). Freely evolving PR–DNS data are at $\rho _p/\rho _f = 100$.

Examining the theoretical inputs provided in (3.5)–(3.6) shows that $\sigma _{a^{\prime \prime }}$ is constant for specified $\{ {Re}_m; \langle \phi \rangle ; \tau _p \}$, but $\tau _{a^{\prime \prime }} = f(T)$ varies in time. To account for the variable memory time scale, the constant-coefficient solutions in (3.3) are integrated forward in time by applying them to a time step $\Delta t$. More specifically, given a specified state at time $s$, the state at time $s + \Delta t$ may be obtained by formally replacing $t=\Delta t$ in (3.3). Repeating the procedure yields $\bar {\boldsymbol {\varSigma }}(t)$, which may be utilized to compute $S(t)$ and $\varGamma (t)$ (see Appendix C).

4. Results

Analytical results are first compared with PR–DNS benchmark data for two canonical flows at ${Re}_m=20$, $\langle \phi \rangle =0.1$ and $\rho _p/\rho _f=1000$. Namely, homogeneous heating systems (HHS) and homogeneous cooling systems (HCS) are examined in detail. In the former systems, particles are given an initial granular temperature $T_0 = 0$, while the latter is initialized as $T_{0} > T_{\infty }$, with $T_{\infty }$ the steady granular temperature. Considering the dimensionless source $\hat {S}(t) = S(t)\,\tau _p/((1-\langle \phi \rangle )^2 | \langle \boldsymbol {w} \rangle |^2 )$ (This normalization is not appropriate for the case of zero mean slip where the mean hydrodynamic drag is zero but there is a purely fluctuating part; such systems could be an interesting canonical case for model development.) and analogous dimensionless sink $\hat {\varGamma }(t)$, we observe acceptable agreement between analytical predictions and PR–DNS data (see figure 3). We note that $T(t)$ is a result of integrating the difference $S(t)-\varGamma (t)$. Therefore, adding a constant to the source and sink at each time $t$ preserves the difference, rendering $T(t)$ invariant. For this reason, differences between $S(t)$ and $\varGamma (t)$ from PR–DNS and theory do not yield significant errors in the granular temperature. Physically speaking, departure between theory and simulations may be due to underlying assumptions in the model process – e.g. isotropic and constant $\sigma _{a^{\prime \prime }}$ or the Markovian Ornstein–Uhlenbeck process with exponential autocorrelation function.

Figure 3. Non-dimensional source (red) and sink (blue) in (a) HHS with $\rho _0=1$, and (b) HCS with $\rho _0=-0.75$, at $\{ {Re}_m=20; \langle \phi \rangle =0.1; \rho _p/\rho _f=1000 \}$. (c) Evolution of non-dimensional granular temperature for HHS (red) and HCS (blue). Analytic solution obtained from (3.4) (lines), PR–DNS data (symbols); insets show the same data with linear scaling.

A direct consequence of characterizing accurately $S(t)$ and $\varGamma (t)$ is that the theory captures the temporal evolution of ${Re}_T$ (see figure 3c). The ability of (3.3) to capture the granular temperature dynamics speaks to the statistical equivalence between the acceleration Langevin model and homogeneous fluidization of elastic particles at finite Reynolds numbers. To demonstrate equivalence at the level of joint statistics, we compare the acceleration–velocity probability distribution to scatter plots obtained from PR–DNS (see figure 4). Again, data extracted from PR–DNS are well characterized by $P(v^{\prime },a^{\prime };t)$ obtained from integration of (3.3). In HHS (figures 4a–c), granular temperature evolution is largely dominated by hydrodynamic sources (quadrants 1 and 3). By contrast, hydrodynamic sinks (quadrants 2 and 4) are dominant in HCS (figures 4d–f) since the systems are initialized with an over-prescribed velocity variance. We emphasize that both systems must satisfy a fluctuation–dissipation relation and converge to a sustained particle velocity variance that results from a balance of source and sink.

Figure 4. Joint p.d.f.s of the fluctuating particle acceleration and fluctuating particle velocity for HHS (a–c) and HCS (d–f). Analytic solution shown by colour, PR–DNS shown by symbols. The first and third quadrants correspond to sources of granular temperature, and the second and fourth quadrants correspond to dissipation.

4.1. Reynolds number and solids volume fraction effects

To examine further the scaling and robustness of solutions presented herein, we consider ${Re}_T(t)$ at ${Re}_m \in [ 10, 100 ]$ with $\langle \phi \rangle =0.1$ and $\rho _p/\rho _f=100$ (see figure 5a). For the conditions considered, the theory predicts granular temperature evolution that is in accordance with PR–DNS results. However, we note that the ${Re}_m=100$ case shows a noticeable over-prediction for ${Re}_T$ at steady state, with relative error ${\approx }18\,\%$. It should be noted that the that the anisotropy in fluid phase kinetic energy has an anomalous dependence at $\langle \phi \rangle = 0.1$, when compared to higher solids volume fractions (Mehrabadi et al. 2015, Tavanashad et al. 2019). The anomalous anisotropy in fluid velocity fluctuations (see figures 2–3 and discussion in § 4.2 of Tavanashad et al. 2019) may in turn affect ${Re}_T$ through anisotropy in the force fluctuations, which are not accounted for in the current model.

Figure 5. (a) Evolution of non-dimensional granular temperature in HHS for varying Reynolds number at $\langle \phi \rangle =0.1$, $\rho _p/\rho _f=100$, with ${Re}_m = 10$ (red), ${Re}_m = 20$ (blue), ${Re}_m = 50$ (green), ${Re}_m = 100$ (black). Evolution of non-dimensional granular temperature for varying solids volume fraction at (b) ${Re}_m =20$, $\rho _p/\rho _f=100$, and (c) $\rho _p/\rho _f=1000$, with $\langle \phi \rangle = 0.1$ (red), $\langle \phi \rangle =0.2$ (blue), $\langle \phi \rangle =0.3$ (green), $\langle \phi \rangle =0.4$ (black). Analytic solution shown by lines, PR–DNS shown by symbols.

Probing the effect of solids volume fraction, we consider ${Re}_T(t)$ at ${Re}_m=20$ with $\langle \phi \rangle \in [ 0.1, 0.4 ]$ and $\rho _p/\rho _f =100$ (see figure 5b). Again, we observe reasonable agreement between theory and PR–DNS. Furthermore, the moderate growth in granular temperature for $\langle \phi \rangle \geqslant 0.2$ is replicated appropriately with (3.6b). Repeating the analysis for $\rho _p/\rho _f =1000$, which yields smaller steady granular temperature, again shows acceptable agreement with PR–DNS (see figure 5c). Thus the present results motivate further the notion that accurate and robust descriptions for $\sigma _{a^{\prime \prime }}$ are of significant value to the present theory.

4.2. Comparison with Stokes flow theory

In this section, we consider comparison with the theory of Koch & Sangani (1999) (KS99), which was developed for inertial particles at moderate solids volume fractions and low Reynolds numbers (Stokes flow regime). The granular temperature balance for homogeneous fluidization of elastic particles (see (2.1)) is closed with KS99 theory according to

(4.1)\begin{equation} \frac{{\rm d} T}{{\rm d} t} \equiv S - \varGamma = \frac{d_p |\langle \boldsymbol{w} \rangle |^2 R_s R_{drag}^2 }{6\sqrt{{\rm \pi} } \tau_p^2 }\,T^{{-}1/2} - \frac{2 R_{diss} }{\tau_p}\,T, \end{equation}

where the coefficients appearing in (4.1) are given by

(4.2)\begin{gather} R_s = \frac{1}{g_0 \left(1 + 3.5 \langle \phi \rangle^{1/2} + 5.9 \langle \phi \rangle \right)}, \end{gather}

(4.3)\begin{gather}R_{drag} = \frac{1 + 3(\langle \phi \rangle/2)^{1/2} + (135/64) \langle \phi \rangle \ln \langle \phi \rangle + 17.14 \langle \phi \rangle}{1 + 0.681 \langle \phi \rangle - 8.48 \langle \phi \rangle^2 + 8.16 \langle \phi \rangle^3}, \end{gather}

(4.4)\begin{align} R_{diss} &= 1 + 3 (\langle \phi \rangle/2)^{1/2} + (135/64) \langle \phi \rangle \ln \langle \phi \rangle + 11.26 \langle \phi \rangle \nonumber\\ & \quad \times \left(1 - 5.1 \langle \phi \rangle + 16.57 \langle \phi \rangle^2 - 21.77 \langle \phi \rangle^3 \right) - \langle \phi \rangle g_0 \ln \varepsilon_m, \end{align}

where $\varepsilon _m=0.01$ is the lubrication breakdown length. For homogeneous fluidization at a fixed Reynolds number, (4.1) yields a Bernoulli differential equation with constant coefficients, and may be solved analytically to obtain

(4.5)\begin{equation} T(t) = \left[ \frac{a}{b} + \left(T_0^{3/2} - \frac{a}{b} \right) \exp \left(- \frac{3 b t }{2} \right) \right]^{2/3}, \end{equation}

where $a=d_p |\langle \boldsymbol {w} \rangle |^2 R_s R_{drag}^2 /( 6\sqrt {{\rm \pi} } \tau _p^2 )$ is the source coefficient, and $b = 2 R_{diss}/ \tau _p$ is the sink coefficient. In the long time limit, (4.5) yields a steady value

(4.6)\begin{equation} T_{\infty} = \left(|\langle \boldsymbol{w} \rangle | R_{drag}\right)^2 \left[ \frac{R_s}{6 \sqrt{\rm \pi} \,{St} \,R_{diss}}\right]^{2/3}, \end{equation}

where we emphasize that the Stokes number ${St} = 2 \tau _p |\langle \boldsymbol {w} \rangle | R_{drag}/d_p$ is computed with the terminal velocity for an isolated sphere subject to Stokes drag, which is related to the slip velocity as $U_t = |\langle \boldsymbol {w} \rangle | R_{drag}$. Equation (4.6) reveals that $T_{\infty } \sim (\rho _p/\rho _f)^{-2/3}$ with KS99 theory.

Employing (4.5) and (4.1), the evolution of sources, sinks, and dimensionless granular temperature may be obtained readily for the conditions considered in figure 3 (see figure 6). Inspection of the source term in (4.1) shows that $S\sim T^{-1/2}$ diverges as $T \rightarrow 0$. For the source-dominated heating system, KS99 theory predicts qualitatively incorrect evolutions for $S(t)$ when compared to PR–DNS (see figure 6a). In fact, the slopes of the lines in figure 6(a) are given by the $T$ exponent in (4.1) (namely, $-1/2, 1$) since the plot is in logarithmic scale. Since the balance in (4.1) may be formulated as an inhomogeneous first-order differential for $U = T^{3/2}$, the weakly singular source for $T$ is mapped to a constant source for $U$, thereby yielding $T(t)$ that is well-behaved in HHS. Therefore, KS99 theory predicts reasonable ${Re}_T (t)$ in HHS, but does so by treating the fluid-mediated source term as a weakly singular impulse, which is physically incorrect. These results show that the hydrodynamic acceleration vector $\boldsymbol {a}_h$, which plays a role in the Enskog equation (1.1), is not properly defined with KS99 theory. For HCS, $S$ and $\varGamma$ behaviour dramatically improves (see figure 6b); this implies that KS99 theory characterizes the hydrodynamic sink better than the source. The leading-order approximation to $\varGamma$ is expected to follow the mean drag closure (a well-characterized quantity). Since KS99 theory employs a mean drag closure for $\varGamma$, improved predictions for dissipation-dominated flows are to be expected.

Figure 6. Same as figure 3, but dashed lines correspond to the theory of Koch & Sangani (1999).

Repeating the parameter sweeps in figure 5 allows us to draw direct comparisons with KS99 theory over a wide range of conditions (see figure 7). At larger Reynolds numbers ${Re}_m > 20$, KS99 theory leads to under-predictions for the steady granular temperature, which is likely a consequence of extending the theory so far from Stokes flow (see figure 7a). At ${Re}_m = 20$, KS99 theory predicts the correct qualitative trends with increasing solids volume fraction and is in reasonable agreement with the present theory (see figures 7b,c).

Figure 7. Same as figure 5, but dashed lines correspond to the theory of Koch & Sangani (1999).

4.3. Density ratio scaling

While the present theory is designed around inertial particles at high density ratios $\rho _p / \rho _f \gg 1$, it is beneficial to consider the behaviour across density ratio. We first note that previous PR–DNS studies at $\rho _p/\rho _f \geqslant 100$ have concluded that $T_{\infty } \sim ( \rho _p/\rho _f )^{-1}$ (Tang, Peters & Kuipers 2016; Tenneti et al. 2016), while the theory of KS99 predicts $T_{\infty } \sim (\rho _p/\rho _f)^{-2/3}$ for inertial particles subject to Stokes flow. The aforementioned scalings from PR–DNS and theory will lead to unrealistic granular temperature for low density particles $\rho _p/\rho _f \ll 1$. The PR–DNS study of Tavanashad et al. (2019) showed that $T_{\infty }$ levels off at low density ratio due to the additional dissipation associated with unsteady hydrodynamic forces (e.g. added mass and Basset history). Furthermore, for ${Re}_m \in [10, 100]$ and $\langle \phi \rangle =0.1$, Tavanashad, Passalacqua & Subramaniam (2021) showed that the correlation of Tenneti et al. (2011) replicated the mean drag forces obtained from freely evolving PR–DNS simulations of buoyant particles down to $\rho _p/\rho _f = 0.01$. Therefore, at ${Re}_m=20, 50, 100$ and $\langle \phi \rangle =0.1$, we expect the $\tau _d$ closure employed here to be reasonably accurate across density ratio. Taking $\lim _{t\to \infty } \sigma _{v^{\prime }}^2(t)$ in (3.3), we obtain an algebraic relation for the steady granular temperature resulting from the acceleration Langevin model:

(4.7)\begin{equation} T_{\infty} = \frac{\sigma_{a^{\prime\prime}}^2 \tau_d^2\, \tau_{a^{\prime\prime}}(T_{\infty})}{\tau_d + \tau_{a^{\prime\prime}} (T_{\infty})}. \end{equation}

Solving (4.7) for $T_{\infty }$ allows us to ascertain the behaviour of the present theory across a wide range of density ratios (see figure 8). Rather than a single power-law scaling, as deduced from inertial particle systems, the present theory predicts an exponent, $n$, that is itself a function of the density ratio, i.e. $T_{\infty } \sim ( \rho _p/\rho _f )^{n(\rho _p/\rho _f)}$. Considering the asymptotic limit of $\rho _p/\rho _f \rightarrow \infty$, one may recognize that $\tau _d \gg \tau _{a^{\prime \prime }}$, thus $T_{\infty } \sim \sigma _{a^{\prime \prime }}^2 \tau _d \tau _{a^{\prime \prime }}(T_{\infty })$. Therefore, in the limit of infinitely large density ratio, the present theory predicts $T_{\infty } \sim (\rho _p/\rho _f)^{-2/3}$, which agrees with KS99. Furthermore, it is stressed that the $n=-2/3$ scaling replicates PR–DNS data for $\rho _p/\rho _f \in [10, 1000]$ better than the correlation of Tenneti et al. (2016) with $n=-1$ scaling. Examining figure 7 in Tenneti et al. (2016), one may also observe the departure from $n=-1$ scaling out to $\rho _p/\rho _f = 2000$. For the asymptotic limit of $\rho _p/\rho _f \rightarrow 0$, one may recognize that $\tau _d \ll \tau _{a^{\prime \prime }}$, thus $T_{\infty } \sim \sigma _{a^{\prime \prime }}^2 \tau _d^2$. Therefore, in the limit of infinitely small density ratio, the present theory predicts $T_{\infty } \sim (\rho _p/\rho _f)^{0}$, which agrees with the PR–DNS of Tavanashad et al. (2019). In summary, the theory described herein exhibits correct scaling for the granular temperature dependence upon density ratio in the asymptotically inertial and buoyant particle limits, and captures qualitatively the magnitude of the steady granular temperature across the range of density ratios considered.

Figure 8. Non-dimensional granular temperature at steady state as a function of density ratio at ${Re}_m = 20$ (red), ${Re}_m = 50$ (blue), ${Re}_m = 100$ (black), with $\langle \phi \rangle = 0.1$ and $\rho _p/\rho _f \in [0.001, 1000]$. Analytical solutions are shown by solid lines, PR–DNS of Tavanashad et al. (2019) are shown by symbols, and $n=-2/3$ scaling is shown by dotted lines.

5. Conclusion

Analytical solutions for fluid-mediated sources of granular temperature were derived from an acceleration Langevin equation, which models drag disturbances in moderately dense suspensions of elastic particles at finite Reynolds numbers. Quadrant-conditioned covariance integrals of the joint fluctuating acceleration–velocity distribution function $P(v^{\prime },a^{\prime };t)$ yield directly the hydrodynamic source and sink. Analytical predictions for the evolution of source, sink and granular temperature are in agreement with data obtained from PR–DNS. Furthermore, the derived $P(v^{\prime },a^{\prime };t)$ agrees with acceleration–velocity scatter plots obtained from PR–DNS, and the theory predicts qualitatively correct steady granular temperature scaling for inertial and buoyant particles.

A variety of approaches may be taken when combining the present work with kinetic theory. The first, and most rigorous, approach would be to follow Garzó et al. (2012) and start off with the Enskog equation in (1.1). The closure for hydrodynamic acceleration then has implications on the equilibrium distribution, the solids transport coefficients and the granular temperature balance (inclusion of $S-\varGamma$). We also emphasize that anisotropic drag forces (observed here) may then lead to an equilibrium distribution that is an anisotropic Gaussian (Vié, Doisneau & Massot 2015), where anisotropy in the particle Reynolds stresses becomes a consequence of drag anisotropy and collisions. A deep dive into this first approach is of great value but beyond the scope of the present work. A second, and far simpler, approach is to append the granular temperature balance to include $S - \varGamma$. Since acceleration–velocity covariances are not transported quantities in standard two-fluid models, evaluation of $\rho _0$ in the present theory poses a challenge for inhomogeneous flows. We consider a local equilibrium approximation for this term in Appendix C.

While results presented here at low density ratio are incipient, the modelling paradigm shows promise as a viable starting place for developing a theory across particle density ratio. In this spirit, improved $\sigma _{a^{\prime \prime }}$ correlations may be obtained and utilized alongside drag correlations that are valid across density ratio (Tavanashad et al. 2021). These efforts may allow extension of the present theory beyond inertial particles, i.e. a general description for fluid-mediated sources to granular temperature that is valid from gas–solid to bubbly flows appears possible.

Finally, it is noted that the present theory is built upon the assumption that hydrodynamic forces within a suspension are normally distributed. While normal distributions have been observed by the authors at many of the considered conditions $\{{Re}_m; \langle \phi \rangle ; \rho _p/\rho _f \}$, as well as by others at similar conditions (Huang et al. 2017; Esteghamatian et al. 2018; Balachandar 2020), the mechanism(s) leading to normally distributed drag fluctuations pose an outstanding question. Future work that expounds upon the normal drag distribution would provide valuable insight into the regions of validity for the present theory.