2. Theory

In any statistical closure problem an important modelling question is the adequacy of the mathematical representation to capture physical phenomena. The work of Hrenya & Sinclair (Reference Hrenya and Sinclair1997) shows that closure only at the level of the mean momentum is not adequate to capture the profiles of mean solid volume fraction and mean particle velocity in core-annular flow, but a closure at the level of the second moment of particle velocities is necessary. The starting point in this work for understanding the generation of particle velocity fluctuations is to consider the evolution of the one-particle distribution function, also referred to as the number density function (NDF).

We consider monodisperse particles moving in a gas under the influence of both hydrodynamic and collisional forces. We consider the one-particle distribution function $f(\boldsymbol{x},\boldsymbol{v},t)$ that is the density of the average number of particles $\langle N\rangle$ in $(\boldsymbol{x},\boldsymbol{v})$ space, such that $\int f(\boldsymbol{x},\boldsymbol{v},t)\,\text{d}\boldsymbol{x}\,\text{d}\boldsymbol{v}=\langle N\rangle$ . If the system is examined on time scales smaller than the duration of each collision (as is the case for the soft-sphere collision model used in the PR-DNS), then the corresponding evolution equation for the one-particle distribution function $f(\boldsymbol{x},\boldsymbol{v},t)$  is given by (Subramaniam Reference Subramaniam2001; Liboff Reference Liboff2003)

where the conditional average of the particle acceleration $\left\langle \boldsymbol{A}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle =\left\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle +\left\langle \boldsymbol{A}_{c}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle$ has been decomposed into contributions arising from both hydrodynamic $\left\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle$ and collisional forces $\left\langle \boldsymbol{A}_{c}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle$ . If a soft-sphere collision model is considered then each collision occurs over a finite time interval, and this is the case for the PR-DNS computations where the duration of each soft-sphere collision is fully resolved. The PR-DNS results presented later in the paper correspond to this soft-sphere NDF equation (2.1).

If a mesoscale time averaging is performed over a time scale that is larger than the duration of each collisional contact, but smaller than the time between collisions, then collision events can be modelled as binary and instantaneous. This leads to the Boltzmann model for a dilute system of hard spheres immersed in a fluid (or the Boltzmann–Enskog model for higher particle solid volume fraction). Note that in these hard sphere models, only the change in impulse can be defined for an instantaneous collision, whereas in the soft-sphere model the particle acceleration due to collision can be meaningfully defined. We now examine the Boltzmann–Enskog transport equation for the one-particle distribution function in gas–solid flow.

For monodisperse particles undergoing binary, instantaneous collisions in a fluid, the one-particle distribution function $f(\boldsymbol{x},\boldsymbol{v},t)$  is defined in a position–velocity space (Williams Reference Williams1958; Koch Reference Koch1990; Subramaniam Reference Subramaniam2000, Reference Subramaniam2013; Liboff Reference Liboff2003; Garzo, Dufty & Hrenya Reference Garzo, Dufty and Hrenya2007; Marchisio & Fox Reference Marchisio and Fox2013), and evolves by the following transport equation:

where $\boldsymbol{{\rm\nabla}}_{\boldsymbol{x}}$ and $\boldsymbol{{\rm\nabla}}_{\boldsymbol{v}}$ denote the gradient operators in the position–velocity space, respectively, $\langle \tilde{\boldsymbol{A}}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\rangle$ is the mesoscale conditional average of the particle acceleration due to hydrodynamic forces, and ${\dot{f}}_{\mathit{coll}}$ is the collisional term that depends on the two-particle density. (Angle brackets here represent averaging over the two-particle density in the kinetic theory of gas–solid flow (Koch Reference Koch1990).) The mesoscale conditional average of the particle acceleration due to hydrodynamic forces $\langle \tilde{\boldsymbol{A}}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\rangle$ in this equation is the mesoscale time-averaged counterpart of the corresponding term in (2.1), but because the hydrodynamic fields change negligibly over the duration of a collision, as a first approximation they can be taken to be identical. So from this point we drop the tilde that distinguishes this quantity between the two equations. A closure model for the collisional term results in a kinetic equation. This well-known equation has been extensively studied in the context of granular flows where collisions are inelastic. Extensions to non-dilute cases that follow the Enskog approach have also been pursued.

The principal difference between the kinetic theory of molecular (or granular) gases and the kinetic theory of gas–solid flow is that in the latter, the conditional particle acceleration due to the hydrodynamic force $\left\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle$ depends on the particle velocity through slip with respect to the fluid, and hence is retained inside the velocity derivative in the velocity transport term. This conditional particle acceleration is caused by the hydrodynamic force $\boldsymbol{F}_{\mathit{fluid}}=m\left\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle$ experienced by the particle due to the pressure and velocity gradient fields at the particle surface. In the transport equation for the distribution function (cf. (2.2)), $\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\rangle$ depends on the instantaneous particle velocity $\boldsymbol{v}$ and can vary in both space and time. For the spatially homogeneous case with monodisperse particles it can be interpreted as the average acceleration experienced by a particle with velocity $\boldsymbol{v}$ . The averaging operator $\langle \cdot \rangle$ represents integration over all higher-order multiparticle distribution functions (Koch Reference Koch1990; Subramaniam Reference Subramaniam2000) that can be defined on the basis of the ensemble of particles with position and velocity $\{\boldsymbol{X}^{(n)}(t),\boldsymbol{V}^{(n)}(t),n=1,\ldots ,N\}$ . In particular, the conditional acceleration $\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\rangle$ is obtained by integrating out its dependence on the two-particle density (pair correlation function). In other words, the conditional acceleration $\langle \boldsymbol{A}_{h}\mid \boldsymbol{x},\boldsymbol{v};t\rangle$ is not completely determined by the particle velocity, but may be affected by the presence of neighbour particles.

Returning to (2.1) that applies to time-resolved collisions as in the PR-DNS, for simplicity we denote the conditional acceleration term $\left\langle \boldsymbol{A}\mid \boldsymbol{x},\boldsymbol{v};t\right\rangle$ representing both hydrodynamic and collisional effects as $\boldsymbol{A}$ in the rest of the paper. An important observation is that a model for $\boldsymbol{A}$ implies a model for the mean interphase momentum transfer as well as a model for the source and dissipation of granular temperature. Note that if we write $\boldsymbol{A}=\left\langle \boldsymbol{A}\right\rangle +\boldsymbol{A}^{\prime \prime }$ , then $\left\langle \boldsymbol{A}\right\rangle$ represents the mean interphase momentum transfer (drag) that has been studied extensively, as noted earlier. Since collisions conserve linear momentum, the collisional contribution to mean interphase momentum transfer is zero. The focus of this paper is to develop a model for the fluctuating particle acceleration that closes the evolution equation of particle granular temperature. We now discuss the evolution of particle granular temperature and identify the source and dissipation terms in statistically homogeneous gas–solid flows.

3. Numerical method

In this study, PR-DNS of freely-evolving gas–solid suspensions are performed using the particle-resolved uncontaminated-fluid reconcilable immersed-boundary method (PUReIBM) that has been developed to simulate flow past fixed particle assemblies (Garg et al. Reference Garg, Tenneti, Mohd-Yusof, Subramaniam, Pannala, Syamlal and O’Brien2011) and freely-evolving suspensions (Tenneti et al. Reference Tenneti, Garg, Hrenya, Fox and Subramaniam2010). In PUReIBM, the entire physical domain is discretized using a uniform Cartesian grid and the governing equations of fluid flow are solved on all the grid points (including those lying inside the particles). The governing equations that are solved in the PUReIBM method are the continuity equation

and the Navier–Stokes equations

Here $\boldsymbol{u}$ is the instantaneous velocity field and $\boldsymbol{g}_{\mathit{IBM}}$ is the pressure gradient. The third term on the right-hand side of (3.2) is the immersed boundary (IB) force $\boldsymbol{f}$ that accounts for the presence of solid particles by ensuring the no slip and no-penetration boundary conditions at the surface of the particles. The final term on the right-hand side of (3.2) is the force per unit volume due to acceleration $\boldsymbol{A}_{f}$ of the reference frame. The governing equations in PUReIBM are formulated in an accelerating frame of reference that moves with the mean velocity of the particles. This formulation enables the simulation of not only sedimenting suspensions, but also suspensions with arbitrary values of mean slip velocity, while maintaining average particle motion at rest in the accelerating frame motion. The full details of the formulation in an accelerating frame of reference are given in the work of Tenneti et al. (Reference Tenneti, Garg, Hrenya, Fox and Subramaniam2010).

As described earlier, the IB force $\boldsymbol{f}$ in (3.2) is used to impose the no slip and no-penetration conditions on the fluid velocity at the fluid–particle interface. Complete details of the computation of the IB force can be found in the work of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011). The surface of the spherical particle is represented by a discrete number of points, called boundary points, that are parametrized in spherical coordinates. Two additional sets of points, termed the exterior and interior points are generated by projecting the boundary points onto spheres of radii $r+{\rm\Delta}r$ and $r-{\rm\Delta}r$ respectively, with ${\rm\Delta}r$ chosen to be equal to the grid spacing. In PUReIBM, the IB force is computed at the interior points so that a desired velocity $\boldsymbol{u}^{(d,k)}$ is obtained at the $k$ th interior point. Following the direct forcing method proposed by Yusof (Reference Yusof1996) the IB force $\boldsymbol{f}^{(k)}$ at the $k$ th interior point is computed as:

In the above expression, $\boldsymbol{S}^{(k,n)}$ is the convective term evaluated at the $k$ th interior point. When the expression for the IB force given in (3.3) is substituted in the discretized version of (3.2), the velocity at the $k$ th interior point is forced to the desired value $\boldsymbol{u}^{(k,d)}$ . Note that all the terms in the above equation are evaluated at the $k$ th interior point. The desired velocity $\boldsymbol{u}^{(k,d)}$ depends on the velocity of the particle. For instance, for a fixed particle the desired velocity at the interior point is equal in magnitude but opposite in direction of the fluid velocity at the corresponding exterior point so that the velocity at the boundary point is zero. The IB force so computed at all the interior points is interpolated to the neighbouring grid nodes that do not include those lying in the fluid phase to get  $\boldsymbol{f}$ .

It has been shown by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) that volume-averaged equations in a periodic domain can be reconciled with the ensemble-averaged equations of statistically homogeneous gas–solid suspensions. Hence the governing equations (3.1) and (3.2) are solved numerically in a cubic domain with periodic boundaries, since the primary goal is to simulate statistically-homogeneous suspensions. Similar to DNS of single-phase turbulence, it is simpler to solve for appropriately defined fluctuating variables. The velocity field is decomposed into a spatially uniform mean flow $\left\langle \boldsymbol{u}\right\rangle _{\mathscr{V}}\left(t\right)$ and a fluctuating velocity field $\boldsymbol{u}^{\prime }$ , i.e.

where the volumetric mean velocity is obtained by averaging the velocity field over the entire computational domain. Similar decompositions are written for the convective term, the pressure gradient $\boldsymbol{g}_{IBM}$ and IB force. Substituting the above decompositions in (3.1) and (3.2), followed by averaging over the entire computational domain, yields the volume-averaged mass and momentum conservation equations (detailed equations are given by Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011)). Evolution equations for the fluctuating variables are derived by subtracting the volume-averaged equations from their instantaneous counterparts. The resulting equations for the fluctuating variables are solved using a pseudo-spectral method, with the Crank–Nicolson scheme for the viscous terms and an Adams–Bashforth scheme for the convective terms. A fractional time-stepping method that is based on Kim and Moin’s approach (Kim & Moin Reference Kim and Moin1985) is used to advance the fluctuating velocity fields in time.

In freely-evolving gas–solid suspensions, the particles move under the influence of hydrodynamic and collisional forces. In PUReIBM, the particles are represented in a Lagrangian frame of reference. The position and translational velocity of the $i$ th particle at time $t$ are denoted $\boldsymbol{X}^{(i)}(t)$ and $\boldsymbol{V}^{(i)}(t)$ , respectively, and they evolve in time as:

where $\boldsymbol{B}$ is any external body force, $\boldsymbol{F}_{h}^{(i)}$ is the hydrodynamic force (calculated from the velocity and pressure fields at the particle surface) and $\boldsymbol{F}_{ij}^{(c)}$ is the contact force on the $i$ th particle as a result of collision with $j$ th particle. Here $N_{p}$ is the total number of particles used in the simulation. Particle collisions are treated using a soft-sphere model originally proposed by Cundall & Strack (Reference Cundall and Strack1979). In the soft-sphere approach, particles are allowed to overlap during a collision, and the contact mechanics between the overlapping particles is modelled by a system of springs and dashpots in both normal and tangential directions. The spring causes the colliding particles to rebound, and the dashpot mimics the dissipation of kinetic energy due to inelastic collisions. The spring stiffness coefficients in the tangential and normal directions are $k_{t}$ and $k_{n}$ , respectively. Similarly, the dashpot damping coefficients in the tangential and normal directions are ${\it\eta}_{t}$ and ${\it\eta}_{n}$ , respectively. The spring stiffness and dashpot damping coefficients are related to the coefficient of restitution and the coefficient of friction (Garg et al. Reference Garg, Galvin, Li and Pannala2010).

The particles considered in this study are assumed to be frictionless. This implies that the tangential component of the contact force is zero, and thus rotation of particles is not accounted for in this study. Therefore, only the normal component of the contact force $\boldsymbol{F}_{\boldsymbol{n}_{ij}}$ is considered at time $t$ and it is given by

where ${\it\delta}_{ij}$ is the overlap between the particles computed using the relation

and $\boldsymbol{V}_{n_{ij}}$ is the relative velocity in the normal direction that is defined using

The normal vector $\hat{\boldsymbol{r}}_{ij}$ is the unit vector along the line of contact pointing from particle $i$ to particle $j$ . The hydrodynamic and contact forces computed at each time step are then used to evolve the position and translational velocity of particles by (3.5) and (3.6), respectively.

The salient feature that distinguishes PUReIBM from other immersed-boundary method approaches is that the IB forcing is restricted to the grid points that lie inside the solid. Consequently the velocity and pressure fields in the fluid phase are a solution to the unmodified Navier–Stokes equations with the appropriate boundary conditions. The hydrodynamic force acting on the particles is computed directly from the stress tensor at the particle surface. Another consequence of the non-contamination of the fluid flow solution by the IB forcing is that the evolution equation of the volume-averaged fluid velocity exactly corresponds to that of the ensemble-averaged fluid-phase velocity. Therefore it is possible to reconcile the solution obtained from PR-DNS with any random field theory of multiphase flow.

4. Validation

The PUReIBM methodology has been validated in a comprehensive suite of test cases: (i) Stokes flow past simple cubic (SC) and face centred cubic (FCC) arrangements (ranging from dilute to close-packed limit) with the boundary-integral method of Zick & Homsy (Reference Zick and Homsy1982), (ii) Stokes flow past random arrays of monodisperse spheres with lattice Boltzmann (LB) simulations of van der Hoef et al. (Reference van der Hoef, Beetstra and Kuipers2005), (iii) moderate to high Reynolds numbers ( $\mathit{Re}_{m}\leqslant 300$ ) in SC and FCC arrangements with LB simulations of Hill et al. (Reference Hill, Koch and Ladd2001b ) and (iv) high Reynolds number flow past random arrays of monodisperse spheres with the ANSYS-FLUENT CFD package where body-fitted grids were used to resolve the flow around particles. It has also been extended to study passive scalar transport, and validated for heat transfer in flow past fixed particle assemblies of monodisperse spheres (Tenneti et al. Reference Tenneti, Sun, Garg and Subramaniam2013).

In addition to the comprehensive validation of the PUReIBM method (Garg et al. Reference Garg, Tenneti, Mohd-Yusof, Subramaniam, Pannala, Syamlal and O’Brien2011; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011), we present additional results here to validate the PUReIBM methodology in the case of a particle sedimenting under gravity in a quiescent fluid. The simulation is performed in a cubic periodic box, whose side is ten times the particle diameter. The sides of the cube form a Cartesian coordinate system with the axes labelled $x$ , $y$ and $z$ . The acceleration due to gravity is confined to the $x$ – $y$ plane and is directed at an angle ${\rm\pi}/16$ to the $x$ -axis, so that the particle does not interact with the wake of its nearest periodic images. The simulations are carried out for various Archimedes numbers $Ar$ defined as

The sedimenting particle attains its terminal velocity when the hydrodynamic force acting on it balances its weight. For a given Archimedes number, the output of the simulation is a Reynolds number based on the terminal velocity and particle diameter. The terminal Reynolds number can be predicted using the Richardson–Zaki (RZ) correlation for the terminal velocity (Richardson & Zaki Reference Richardson and Zaki1954). The results presented in table 1 show that the terminal velocity obtained in PUReIBM simulations are in good agreement with that predicted by the RZ correlation. These results further establish the capability of PUReIBM to perform simulations of freely-evolving gas–solid suspensions. A summary of various applications of the PUReIBM PR-DNS strategy is given in Tenneti & Subramaniam (Reference Tenneti and Subramaniam2014). The simulation set-up that is used to perform PR-DNS of freely-evolving suspensions to quantify the steady-state source and dissipation of granular temperature, is discussed in the next section.

Table 1.

Comparison of the terminal Reynolds number obtained from PUReIBM simulations with that predicted by the RZ correlation for different Archimedes numbers.

5. Simulation set-up

The independent dimensionless parameters relevant to the study of freely-evolving gas–solid suspensions are the solid volume fraction ${\it\phi}$ , the solid-to-fluid density ratio ${\it\rho}_{p}/{\it\rho}_{f}$ and the Reynolds number based on the mean slip velocity $\mathit{Re}_{m}$ . The mean-slip Reynolds number is defined as

where ${\it\mu}_{f}$ is the dynamic viscosity of the fluid phase. The mean slip velocity $\left\langle \boldsymbol{W}\right\rangle$ is defined as $\langle \boldsymbol{W}\rangle =\langle \boldsymbol{v}\rangle -\langle \boldsymbol{u}^{(f)}\rangle$ , and is the relative velocity between the solid-phase mean velocity and the fluid-phase mean velocity. Averaged quantities in these PR-DNS of statistically homogeneous flow are computed by first volume averaging the flow variable (e.g. $\langle \boldsymbol{u}^{(f)}\rangle$ ) for one realization (corresponding to an initial particle configuration), and subsequently ensemble averaging over different particle configurations (corresponding to a specified average solid volume fraction, mean-slip Reynolds number and initial pair correlation function).

For a given triplet $\left({\it\phi},\mathit{Re}_{m},{\it\rho}_{p}/{\it\rho}_{f}\right)$ , the gas–solid suspension attains a steady granular temperature. The simulation set-up in the present study has been chosen to supply arbitrary values of $\left({\it\phi},\mathit{Re}_{m},{\it\rho}_{p}/{\it\rho}_{f}\right)$ as input, and keep them fixed during the course of the simulation. The task of maintaining a given $\mathit{Re}_{m}$ can be accomplished by specifying a mean pressure gradient that does not exactly balance the weight of the particles, but exerts the requisite body force to maintain the desired slip velocity between the particles and fluid. However, now both the mean particle velocity and the mean fluid velocity change in time because there is no steady solution in the laboratory frame to the mean momentum balance in each phase. Note that even though the mean phasic velocities are evolving in time, their difference – the mean slip velocity – attains a steady value.

The difficulty in simulating this flow set-up in the laboratory frame with periodic boundary conditions is that the continuous increase in fluid and particle velocities places unnecessary restrictions on the time step through the Courant condition. To circumvent this problem we use a different simulation set-up that performs the PR-DNS in an accelerating reference frame such that the particles have a zero mean velocity with respect to the computational grid. The advantage of this set-up is that the desired mean-slip Reynolds number is specified as an input parameter, and we are able to solve the problem with reasonable time steps that resolve the flow.

A homogeneous initial particle configuration is generated for the given solid volume fraction ${\it\phi}$ . The initial positions of the particles are obtained by equilibrating an ensemble of particles undergoing elastic collisions in the absence of interstitial fluid. The particles are advanced on a time step ${\rm\Delta}t_{\mathit{coll}}$ that is determined by the spring stiffness and the dashpot coefficients. The flow fields are updated on a time step ${\rm\Delta}t_{\mathit{fluid}}$ , which ensures that both the convective and viscous time scales are well resolved. At the start of a flow time step the forces acting on the particles are computed based on the flow fields obtained at the end of the previous flow time step. If ${\rm\Delta}t_{\mathit{coll}}$ is smaller than ${\rm\Delta}t_{\mathit{fluid}}$ the particles are stepped by ${\rm\Delta}t_{\mathit{coll}}$ until the end of the flow time step, otherwise both the particles and the fluid are stepped by ${\rm\Delta}t_{\mathit{coll}}$ . The simulation is continued until the granular temperature reaches a steady state.

The computational domain used is a cube with sides of length $L$ which is discretized using a regular Cartesian grid with $M$ grid cells in each direction so that ${\rm\Delta}x=L/M$ is the size of each grid cell. The spatial resolution is represented by the number of grid cells across the diameter of a particle, which is denoted $D_{m}=d_{p}/{\rm\Delta}x$ . For random arrangements of particles, the ratio $L/d_{p}$ is an independent parameter. The minimum box length is determined by the criterion that the spatial autocorrelation of flow statistics must decay to zero within the box. This is to prevent the periodicity of the numerical solution from leading to unphysical flow fields. The numerical parameter $L/d_{p}$ also determines the number of particles $N_{p}$ in the box such that for a given volume fraction ${\it\phi}_{}$ it is given by

where the square brackets denote the nearest integer. The various numerical parameters used in the simulations are reported in table 2.

Table 2.

Numerical parameters (the ratio of the length of the box to the particle diameter $L/d$ , particle diameter in grid units $D_{m}$ , and number of particles $N_{p}$ ) used in the simulations.

In the present study, PR-DNS of freely-evolving gas–solid suspensions are performed for solid volume fractions between 0.1 and 0.4 and for mean-slip Reynolds numbers between 10 and 100. The solid-to-fluid density ratio has been varied from 100 to 2000. Note that the lower limit of the density ratio corresponds to cork particles in air, while the upper limit is characteristic of glass beads suspended in air. The actual values of the physical parameters used in the simulations are provided in table 3. Simulations have been performed for each set of $\left({\it\phi},\mathit{Re}_{m},{\it\rho}_{p}/{\it\rho}_{f}\right)$ values, resulting in a total of 80 cases. For several cases, multiple independent simulations have been performed to ensure statistical convergence. Thus, more than 200 simulations are used to quantify the steady values of granular temperature, and the source and dissipation of granular temperature. These results are described in the subsequent sections.

Table 3.

Physical parameters (solid volume fraction, mean-slip Reynolds number, initial Reynolds number based on granular temperature and solid–fluid density ratio) used in PUReIBM simulations.

6. Quantification of steady-state source and dissipation of granular temperature

The simulation set-up described in the previous section has been used to perform PR-DNS of freely-evolving gas–solid suspensions initialized with different values of the granular temperature. When characterizing the effect of particle velocity fluctuations it is useful to define a Reynolds number based on the granular temperature $\mathit{Re}_{T}$ as:

Estimation of the granular temperature from PR-DNS is discussed in detail by Tenneti et al. (Reference Tenneti, Garg, Hrenya, Fox and Subramaniam2010). The temporal evolution of $T/\!\left|\left\langle \boldsymbol{W}\right\rangle \right|^{2}$ and $\mathit{Re}_{T}$ for a freely-evolving gas–solid suspension at a volume fraction of 0.1 and solid-to-fluid density ratio of 100, with different mean-slip Reynolds numbers is shown in figures 4(a) and 4(b), respectively. In these simulations all the particles start from rest and hence the initial $\mathit{Re}_{T}$ is zero. We observe that the granular temperature reaches a statistical steady state for all values of the mean-slip Reynolds number. In these simulations, a mean pressure gradient is applied to obtain the desired $\mathit{Re}_{m}$ , which results in a mean force $\left\langle \boldsymbol{F}\right\rangle$ acting on the suspension. The energy generated by the application of the mean pressure gradient, $\left\langle \boldsymbol{F}\right\rangle \boldsymbol{\cdot }\left\langle \boldsymbol{W}\right\rangle$ , is responsible for generation of velocity fluctuations in the fluid as well as in the particle phase. It has been shown by Mehrabadi et al. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015) that the fluid velocity fluctuations also attain a steady value.

Figure 4.

(a) Temporal evolution of granular temperature normalized by the square of the magnitude of mean slip velocity for different mean-slip Reynolds number. The solid volume fraction is 0.1 and the solid to fluid density ratio of the suspension is 100. (b) Evolution of the Reynolds number based on granular temperature ( $\mathit{Re}_{T}$ ) for the cases shown in figure 4(a). Error bars are the 95 % confidence intervals obtained from five independent simulations.

From figure 4 we can also see that at a given solid-to-fluid density ratio and solid volume fraction, the normalized granular temperature decreases with increasing $\mathit{Re}_{m}$ (see figure 4 a, whereas $\mathit{Re}_{T}$ increases according to figure 4 b). This result implies that $T/\!\left|\left\langle \boldsymbol{W}\right\rangle \right|^{2}\sim 1/{\mathit{Re}_{m}}^{p}$ , such that $p<2$ . This result also implies that the granular temperature increases with mean slip velocity, and that is expected because particles pick up more energy from the fluid as $\mathit{Re}_{m}$ increases.

A pertinent question that needs to be answered is whether the steady granular temperature attained by the suspension for a given $\left({\it\phi},\mathit{Re}_{m},{\it\rho}_{p}/{\it\rho}_{f}\right)$ depends on the choice of the initial particle granular temperature. To answer this question, we selected a gas–solid suspension with ${\it\phi}=0.1$ , $\mathit{Re}_{m}=20,{\it\rho}_{p}/{\it\rho}_{f}=100$ , and initialized the particles with a Maxwellian velocity distribution corresponding to $\mathit{Re}_{T}=2$ , 4, and 6. The temporal evolution of $\mathit{Re}_{T}$ is shown in figure 5(a). The evolution of $\mathit{Re}_{T}$ for a suspension initialized with a zero particle granular is also shown in the figure for reference. It is clear from the figure that the statistical steady state attained by the suspension is independent of the initial granular temperature. It is useful to examine the behaviour of source $S$ of (2.12) and dissipation ${\it\Gamma}$ of (2.13) with granular temperature for these suspensions to gain a better understanding of the flow physics. Figure 5(b) shows a phase-space plot $({\dot{T}},T)$ of (2.6) where the time derivative of temperature ${\dot{T}}$ is partitioned into source and dissipation terms. Figure 5(b) shows how all the gas–solid systems initialized with different granular temperatures arrive at the same stable steady state. The source and dissipation shown in this figure are computed from PR-DNS data using the estimators given in (B 1) and (B 2). In the suspension initialized with $\mathit{Re}_{T}=0$ , the particles pick up energy from the fluid and hence the source term (open triangles) is greater than dissipation (filled triangles) at initial time. Both source and dissipation grow from initially low values to their steady-state values. For the suspensions initialized with a granular temperature higher than the steady value, the particles have higher energy and thus they lose their energy to fluid. For these suspensions the dissipation (filled diamonds and filled circles) is initially greater than the source term (open diamonds and open circles).

Figure 5.

(a) Temporal evolution of scaled granular temperature for suspensions initialized with different particle granular temperatures. The solid volume fraction is 0.1, $\mathit{Re}_{m}$ is 20 and ${\it\rho}_{p}/{\it\rho}_{f}$ is 100. The error bars represent the 95 % confidence intervals obtained from five independent realizations for each case. (b) Variation of source and dissipation with granular temperature for the suspensions shown in figure 5(a). The suspension initialized with $\mathit{Re}_{T}=2$ is not shown here. Filled symbols denote dissipation and open symbols denote the source.

Having established that the steady state attained by a gas–solid suspension is unique for a given $\left({\it\phi},\mathit{Re}_{m},{\it\rho}_{p}/{\it\rho}_{f}\right)$ , we now determine how the steady granular temperature $T_{SS}$ depends on these physical parameters. To the authors’ knowledge, this is the first report of the dependence of granular temperature on $\left({\it\phi},\mathit{Re}_{m},{\it\rho}_{p}/{\it\rho}_{f}\right)$ . Figure 6 shows the variation of $T_{SS}/\left|\left\langle \boldsymbol{W}\right\rangle \right|^{2}$ with solid volume fraction and mean-slip Reynolds number for a solid-to-fluid density ratio of 100. The steady scaled granular temperature decreases with both volume fraction and mean-flow Reynolds number. Figure 7 shows the behaviour of the steady scaled granular temperature with solid-to-fluid density ratio for a solid volume fraction of 0.2 and mean-flow Reynolds number of 20. As expected, the steady scaled granular temperature decreases with increasing inertia of the particles.

Figure 6.

Variation of steady granular temperature with volume fraction and mean-flow Reynolds number for a density ratio of 100. Symbols are data obtained from PUReIBM PR-DNS and the solid lines indicate the fit given by (6.2). The error bars represent the 95 % confidence intervals obtained from five independent realizations for each case.

Figure 7.

Variation of steady granular temperature with solid to fluid density ratio for a solid volume fraction of 0.1 and mean-flow Reynolds number 20. Symbols denote the data obtained from PUReIBM PR-DNS and the solid line indicates the function $\left({\it\rho}_{p}/{\it\rho}_{f}\right)^{-1}$ . The error bars represent the 95 % confidence intervals obtained from five independent realizations for each case.

It is interesting to note that the behaviour of the steady scaled granular temperature with $\mathit{Re}_{m}$ is similar to the behaviour of the scaled fluid kinetic energy with $\mathit{Re}_{m}$ reported by Mehrabadi et al. (Reference Mehrabadi, Tenneti, Garg and Subramaniam2015). Hence the functional dependence of the scaled $T_{SS}$ on $\mathit{Re}_{m}$ is similar to that proposed for the fluid kinetic energy (Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015). As shown in figure 7, $\left({\it\rho}_{p}/{\it\rho}_{f}\right)^{-1}$ captures the dependence of $T_{SS}$ on the solid-to-fluid density ratio very well. Based on the data obtained from PUReIBM PR-DNS, the following functional form is proposed for the steady value of the scaled granular temperature:

which is valid in the range $0.1\leqslant {\it\phi}\leqslant 0.4$ , $10\leqslant \mathit{Re}_{m}\leqslant 100$ and $100\leqslant {\it\rho}_{p}/{\it\rho}_{f}\leqslant 2000$ .

7. Lagrangian model for hydrodynamic acceleration of particles

As discussed earlier, the closure model for the conditional particle acceleration affects both the mean momentum and the particle granular temperature, as these are derived from the NDF equation (2.1). A closure at the one-particle distribution level automatically implies closure of the mean momentum and particle velocity second-moment equations. Furthermore, closures at the one-particle distribution level are guaranteed to satisfy realizability criteria, whereas special care is needed to ensure the same in the case of moment closures. These considerations motivate the development of models for the unclosed terms in the transport equation for the one-particle distribution function corresponding to gas–solid flow.

The model for the instantaneous particle acceleration can be deduced by considering the factors that contribute to the hydrodynamic forces acting on a particle. In a typical gas–solid flow system, the presence of numerous particles moving in different directions will lead to continuously changing hydrodynamic interactions between particles (i.e. fluctuations in the fluid velocity and pressure fields). We propose the following stochastic model for the increment in the particle velocity:

The above equation is an isotropic form of the generalized Langevin model (Haworth & Pope Reference Haworth and Pope1986). In the above equation $\text{d}v_{i}^{}$ is the increment in the particle velocity, $v_{i}^{\prime \prime }$ is the fluctuation in the particle velocity and $\text{d}\mathscr{W}_{i}$ is a Wiener process increment. Fluctuations in the particle velocity are defined about the mean particle velocity, i.e. $v_{j}^{^{\prime \prime }}=v_{j}-\langle v_{j}\rangle$ . The model coefficients ${\it\beta}$ , ${\it\gamma}$ and $B$ will be specified in § 7.2. The first term on the right-hand side of (7.1) represents the mean acceleration experienced by the solid phase in the suspension, and accounts for the effect of the mean slip velocity. This term is related to the normalized force $F$ as defined by Hill et al. (Reference Hill, Koch and Ladd2001b ) and Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011) as follows:

The coefficient ${\it\beta}$ can thus be inferred by computing $F$ using any of the drag laws developed for gas–solid suspensions (Hill et al. Reference Hill, Koch and Ladd2001a ,Reference Hill, Koch and Ladd b ; van der Hoef et al. Reference van der Hoef, Beetstra and Kuipers2005; Beetstra et al. Reference Beetstra, van der Hoef and Kuipers2007; Yin & Sundaresan Reference Yin and Sundaresan2009a ,Reference Yin and Sundaresan b ; Holloway et al. Reference Holloway, Yin and Sundaresan2010; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011). The second term accounts for the fluctuation in the particle velocity and the last term models the effect of the hydrodynamic interaction of the neighbouring particles. The coefficient ${\it\gamma}$ is the inverse of the Lagrangian particle velocity autocorrelation time. It quantifies how long a particle retains memory of its initial velocity. These coefficients are functions of volume fraction ${\it\phi}$ , mean-flow Reynolds number $\mathit{Re}_{m}$ and particle-to-fluid density ratio ${\it\rho}_{p}/{\it\rho}_{f}$ . It is interesting to note that in the different problems of dilute, turbulent particle-laden Couette (Goswami & Kumaran Reference Goswami and Kumaran2010a ,Reference Goswami and Kumaran b ) and channel flows (Goswami & Kumaran Reference Goswami and Kumaran2011, Reference Goswami and Kumaran2010c ), Goswami & Kumaran found that the effect of fluid velocity fluctuations on particle acceleration can be described using a Langevin equation. The Langevin model for the particle acceleration we propose here is verified by computing the particle velocity autocorrelation function after the granular temperature reaches a steady state. Following a standard approach in PDF modelling of turbulent reacting flows (Pope Reference Pope1985, Reference Pope2000), this stochastic model based on the Langevin equation is developed for a statistically stationary problem with model coefficients that depend on moments of the PDF (e.g. mean velocity or second moments). The model is intended to be then applied to the non-stationary problem with local values of the temporally varying mean and second moments that appear in the model coefficients. This approach makes the assumption of local statistical stationarity.

8. Model assessment

In order to assess the performance of the proposed model, we implemented the particle acceleration model presented in the previous section in MFIX-DEM. In this code, the change in the velocity of the particles due to the hydrodynamic forces is given by (7.1). The change in velocity due to inter-particle collisions is modelled using the same soft-sphere collision model that was described in the PR-DNS numerical method. For a direct assessment of the performance of the model, the simulations are performed in the same periodic cubic domains that were used in the PR-DNS. Moreover, only a single cell is used so as to remove the effects of spatial discretization on the model assessment. Details about the numerical convergence of the implemented model with respect to time step and the number of particles are given in appendix C.

The first objective is to assess the impact of the particle acceleration model on the evolution of granular temperature. For this purpose, we selected a gas–solid suspension with $\mathit{Re}_{m}=10$ , ${\it\phi}=0.1$ , and ${\it\rho}_{p}/{\it\rho}_{f}=100$ , and initialized the particles with different initial granular temperatures. We performed these simulations with both the Lagrangian model for particle acceleration due to hydrodynamic forces that is based on simple extension of the mean drag law (see appendix A for details) and with the Langevin model proposed in this work. The evolution of granular temperature is shown in figure 12 (the results with the Lagrangian model based on simple extension of the mean drag law are denoted ‘without LM’ in the legend). We see that for different initial granular temperatures, the steady-state granular temperature obtained is the same. Moreover, the figure shows that, as expected from the analysis in appendix A, the granular temperature for the Lagrangian model based on simple extension of the mean drag law decays to zero as there is no source term in the evolution equation of the granular temperature. Only by using the stochastic Lagrangian model for hydrodynamic acceleration is a finite steady granular temperature obtained. Also, as a verification, the moment equation (2.11) was solved using Matlab and the comparison shows that the MFIX-DEM results are in accordance with the moment equation.

Figure 12.

Evolution of granular temperature with different initial granular temperatures.

Next we study how well the steady state proposed in § 6 is attained by gas–solid suspensions that are simulated using the acceleration model. For this we performed MFIX-DEM simulations at different solid volume fractions and mean-slip Reynolds numbers at a density ratio of 100. Figure 13 shows that the model predicts the steady-state granular temperature values very close to the fit proposed in (6.2). The maximum error in MFIX-DEM simulation results is around 5 % from the correlation value which is reasonable. As the correlation is obtained from the PUReIBM data corresponding to a wide range of volume fractions and mean-slip Reynolds numbers, the model could be applied to various regions in a riser flow in continuum simulations to compute the local granular temperature based on the local (cell) mean-slip Reynolds number and local volume fraction. This in turn would lead to a better prediction of particle velocities which is essential for understanding the mean flow structure in riser flows.

Figure 13.

Comparison of steady-state granular temperatures for various ${\it\phi}$ and $Re_{m}$ obtained from MFIX-DEM with the correlation proposed in (6.2).