2. Gas-phase velocity variance

In the Eulerian two-fluid theory the fluid-phase Reynolds stress is defined as a phasic average, which is an average conditional on the presence of the fluid phase (Drew Reference Drew1983; Drew & Passman Reference Drew and Passman1998; Pai & Subramaniam Reference Pai and Subramaniam2009). If $Q\left(\boldsymbol{x},t\right)$ is any field, then its phasic average $\langle Q^{(f)}\rangle \left(\boldsymbol{x},t\right)$ , referred to as its fluid-phase mean, is defined as

where $I_{f}$ is the fluid-phase indicator function defined in § 1.

Using this definition, the ensemble-averaged kinetic energy in the fluid phase $\langle E^{(f)}\rangle$ is defined as

where $\boldsymbol{u}$ is the fluid velocity. It is easy to see that the average kinetic energy in the fluid phase is the sum of the kinetic energy in the mean fluid motion $\bar{E}^{(f)}$ and the average kinetic energy in the fluctuating motions $k^{(f)}$ . The average kinetic energy in the mean fluid motion is given by $\bar{E}^{(f)}=(\langle u_{i}^{(f)}\rangle \langle u_{i}^{(f)}\rangle )/2$ , where the quantity $\langle u_{i}^{(\,f)}\rangle$ is the phase-averaged fluid velocity. The average kinetic energy in the fluctuating motion of the fluid is given by

where fluctuations in the fluid velocity field are defined with respect to the phase-averaged fluid velocity, i.e. $u_{i}^{\prime \prime (f)}=u_{i}-\langle u_{i}^{(\,f)}\rangle$ . We now describe the procedure to compute $k^{(f)}$ from the solution of the flow past statistically homogeneous suspensions using PR-DNS.

3. Numerical method

Several numerical methods have been developed for PR-DNS of gas–solid suspensions (Johnson & Tezduyar Reference Johnson and Tezduyar1997; Prosperetti & Oguz Reference Prosperetti and Oguz2001; Peskin Reference Peskin2002). The PR-DNS methodology employed in this work is called the particle-resolved uncontaminated-fluid reconcilable immersed boundary method (PUReIBM).

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 PUReIBM are the continuity equation

and the Navier–Stokes equations

Here, $\boldsymbol{u}$ is the instantaneous velocity field, $\boldsymbol{S}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(\boldsymbol{u}\boldsymbol{u}\right)$ is the convective term in conservative form, $\boldsymbol{g}_{IBM}$ is the pressure gradient and $\boldsymbol{A}_{f}$ accounts for the acceleration of the frame of reference (Tenneti et al. Reference Tenneti, Garg, Hrenya, Fox and Subramaniam2010). The simulation of a freely evolving suspension in PUReIBM is carried out in an accelerating frame of reference that moves with the mean velocity of the particles. The formulation of PUReIBM in an accelerating frame of reference enables the simulation of not just sedimenting suspensions, but arbitrary mean slip values, while maintaining average particle motion at rest in the accelerating frame. In other words, the use of an accelerating frame of reference facilitates examination of the influence of each flow parameter while maintaining other parameters constant.

The no-slip and no-penetration conditions on the fluid velocity at the particle surfaces are imposed via the additional immersed boundary (IB) force term $\boldsymbol{f}$ in (3.2). The full details of the computation of the IB force can be found in the work of Tenneti et al. (Reference Tenneti, Garg and Subramaniam2011). To compute the IB force, the surface of the sphere is represented by a discrete number of points called the boundary points that are parametrized in spherical coordinates. Two additional sets of points, termed as 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}^{(k,d)}$ is obtained at the $k\text{th}$ interior point. Following the direct forcing method proposed by Mohd-Yusof (Reference Mohd-Yusof1996) the IB force $\boldsymbol{f}^{(k)}$ at the $k\text{th}$ interior point is specified to cancel the remaining terms in the momentum conservation, and to force the velocity to its desired value $\boldsymbol{u}^{(k,d)}$ :

It should be noted that all the terms in the above equation are evaluated at the $k\text{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 to 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 obtain $\boldsymbol{f}$ .

The governing equations in PUReIBM are solved in terms of fluctuating variables, and periodicity conditions are imposed on the appropriate fluctuating fields. The velocity field is decomposed into a spatially uniform mean flow $\langle \boldsymbol{u}\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 $\boldsymbol{S}$ , pressure gradient $\boldsymbol{g}$ and IB forcing $\boldsymbol{f}$ terms. Substitution of 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 the case of fixed assemblies, the no-slip and no-penetration boundary conditions on $\boldsymbol{u}$ are imposed at the fixed particle surfaces. For a freely evolving suspension, the particles in PUReIBM move under the influence of hydrodynamic and collisional forces, and are represented in a Lagrangian frame of reference at time $t$ by $\{\boldsymbol{X}^{(i)}(t),\boldsymbol{V}^{(i)}(t)\;i=1\dots \,N_{p}\}$ , with $\boldsymbol{X}^{(i)}(t)$ and $\boldsymbol{V}^{(i)}(t)$ being the position and velocity of the $i$ th particle respectively and $N_{p}$ being the total number of particles. The position and translational velocity of the $i$ th particle evolve according to Newton’s second law 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 the $j$ th particle. Particle–particle interactions are treated using a soft-sphere model originally proposed by Cundall & Strack (Reference Cundall and Strack1979). In the soft-sphere approach, the contact mechanics between two overlapping particles is modelled by a system of springs and dashpots in both normal and tangential directions. The spring causes 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 Pannala2010a ).

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 which 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 IB method approaches (including the original implementation of Mohd-Yusof Reference Mohd-Yusof1996) is that the IB forcing is restricted to those grid points that lie in 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. Since a flow field that is contaminated by additional forcing affects the calculation of $k^{(f)}$ and its dissipation rate, the uncontaminated flow field of PUReIBM is better suited for quantification of $k^{(f)}$ and its dissipation rate. 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 PR-DNS methodology has been extensively validated (Garg et al. Reference Garg, Tenneti, Mohd-Yusof, Subramaniam, Pannala, Syamlal and O’Brien2010b ; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011) by comparing the drag force obtained from PUReIBM with available experimental and simulation data in the literature in the following comprehensive suite of test cases.

1. (a) The drag acting on a single sphere (Garg Reference Garg2009; Garg et al. Reference Garg, Tenneti, Mohd-Yusof, Subramaniam, Pannala, Syamlal and O’Brien2010b ) with the experimental correlation of Schiller & Naumann (Reference Schiller and Naumann1935).

2. (b) The drag acting on simple cubic (SC) and face-centred cubic (FCC) arrangements (Tenneti et al. Reference Tenneti, Garg and Subramaniam2011) of particles in the Stokes flow regime with the results reported by Zick & Homsy (Reference Zick and Homsy1982) using the boundary integral method (semi-analytic solution).

3. (c) The drag acting on SC and FCC arrangements (Tenneti et al. Reference Tenneti, Garg and Subramaniam2011) of particles at moderate Reynolds numbers with the results published by Hill et al. (Reference Hill, Koch and Ladd2001b ) using the lattice Boltzmann method (LBM).

4. (d) The mean drag acting on a random arrangement (Tenneti et al. Reference Tenneti, Garg and Subramaniam2011) of particles in the Stokes flow regime with the results published by Hill et al. (Reference Hill, Koch and Ladd2001a ) and van der Hoef et al. (Reference van der Hoef, Beetstra and Kuipers2005) using the LBM.

5. (e) High-Reynolds-number flow past random arrays of monodisperse spheres with the ANSYS-FLUENT CFD package.

In addition to the comprehensive validation of the PUReIBM method (Garg et al. Reference Garg, Tenneti, Mohd-Yusof, Subramaniam, Pannala, Syamlal and O’Brien2010b ; Tenneti et al. Reference Tenneti, Garg and Subramaniam2011), we present selected additional validation tests to establish the numerical convergence and accuracy of PUReIBM near solid boundaries in figure 2. Figure 2(a) shows a comparison of the pressure coefficient along the surface of a sphere obtained from our PR-DNS with that reported by Clift, Grace & Weber (Reference Clift, Grace and Weber1978) (CGW) for an isolated sphere at a Reynolds number of 10. Figure 2(a) shows an excellent agreement of the coefficient of pressure profile on the surface of the sphere with the data reported in CGW. Figure 2(b) shows a comparison of the velocity field in a square duct at a Reynolds number of 20 with the analytical solution given by Cornish (Reference Cornish1928). We can see that the velocity profile obtained from PUReIBM is numerically converged and accurate. These plots show that in addition to computing the total drag accurately, our method also computes the correct contributions of pressure and viscous drag forces. In the following section we describe the simulation set-up used to compute the level of gas-phase velocity fluctuations and also discuss the choice of the numerical parameters needed to ensure numerically converged results.

Figure 2.

(a) Variation of the coefficient of pressure $C_{p}$ along the surface of the sphere. The symbols are the data obtained from PUReIBM simulations for a Reynolds number of 10 while the solid line is that reported in the book authored by Clift, Grace and Weber (CGW) (Clift et al. Reference Clift, Grace and Weber1978). (b) Comparison of the velocity profile in a square duct obtained from PUReIBM simulations at a Reynolds number of 20 with the analytical solution (Cornish Reference Cornish1928). It should be noted that the walls are generated using the IB method.

5. Simulation set-up

In our simulation set-up the particles are held stationary and a steady flow is established by imposing a pressure gradient that corresponds to the desired mean slip Reynolds number. A typical simulation of flow past a random arrangement of particles with contours of local kinetic energy $({u_{i}^{\prime \prime }}^{(f)}{u_{i}^{\prime \prime }}^{(f)})/2$ normalized by the mean energy is shown in figure 3. In all the simulations the mean flow is directed along the positive $x$ -axis.

For flow past homogeneous particle assemblies, a Reynolds number based on the magnitude of the mean slip velocity between the two phases is defined as

where $|\langle \boldsymbol{W}\rangle |$ is the magnitude of the mean slip velocity, $d_{p}$ is the particle diameter and ${\it\phi}$ is the solid volume fraction. The mean slip velocity $\langle \boldsymbol{W}\rangle =\langle \boldsymbol{u}^{(p)}\rangle -\langle \boldsymbol{u}^{(f)}\rangle$ is defined as the difference between the average solid- and gas-phase velocities. In the simulations, the mean slip Reynolds number (or the desired flow rate) is specified as an input and since for fixed assemblies $\langle \boldsymbol{u}^{(\,p)}\rangle =0$ , the desired fluid-phase mean velocity $\langle \boldsymbol{u}^{(\,f)}\rangle$ is known in terms of the input Reynolds number and other physical properties. The mean pressure gradient evolves in time until it attains the value required to drive the fluid at the desired flow rate.

Figure 3.

Contours of local kinetic energy $({u_{i}^{\prime \prime }}^{(f)}{u_{i}^{\prime \prime }}^{(f)})/2$ in the gas phase normalized by the mean energy for steady flow past a random assembly of spheres at a solid volume fraction of 0.1 and mean slip Reynolds number of 50.

Particles are initialized corresponding to a specified mean solid volume fraction ${\it\phi}$ . The particles are fixed in a random equilibrium configuration they attain following elastic collisions (in the absence of ambient fluid), starting from a lattice arrangement with a Maxwellian velocity distribution. The elastic collisions are simulated using a soft-sphere discrete element model (Cundall & Strack Reference Cundall and Strack1979; Garg et al. Reference Garg, Galvin, Li and Pannala2010a ).

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 1.

The convergence of $k^{(f)}/\bar{E}^{(f)}$ with respect to the grid resolution $D_{m}$ for a solid volume fraction of 0.3 and mean slip Reynolds number of 20 is shown in figure 4. The grid resolutions used in the PUReIBM simulations (see table 1) have been chosen such that they yield numerically converged solutions. The values provided in table 1 guarantee that the particle mean drag force is also converged (Tenneti et al. Reference Tenneti, Garg and Subramaniam2011). The level of $k^{(f)}/\bar{E}^{(f)}$ for a specific random particle assembly clearly shows numerical convergence as the grid resolution is increased.

Figure 4.

Convergence characteristics of $k^{(f)}/\bar{E}^{(f)}$ with grid resolution $d_{p}/{\rm\Delta}x$ for flow past a random particle configuration at ${\it\phi}=0.3$ and $\mathit{Re}_{m}=20$ . The random configuration is the same for all grid resolutions.

Table 1.

Numerical parameters (number of particles $N_{p}$ , number of MIS $\mathscr{M}$ , particle diameter in grid units $D_{m}$ and the ratio of the length of the box to the particle diameter $L/d_{p}$ ) used for random arrays in PUReIBM simulations. Different numerical parameters are used for $\mathit{Re}_{m}\leqslant 100$ and $\mathit{Re}_{m}>100$ . These are separated by a ‘/’. Numbers before the ‘/’ correspond to $\mathit{Re}_{m}\leqslant 100$ while numbers after the ‘/’ correspond to $\mathit{Re}_{m}>100$ . At volume fraction 0.5 PUReIBM simulations are performed only up to a Reynolds number of 100.

The convergence of the ensemble-averaged $k^{(f)}$ with 65 realizations for the case ${\it\phi}=0.4$ and $Re_{m}=100$ is shown in figure 5. This figure indicates that the ensemble-averaged mean obtained with 65 realizations lies within the 95 % confidence intervals obtained with five realizations (highlighted as the filled symbol in figure 5). Moreover, the 95 % confidence intervals for five realizations represent less than 5 % deviation from the evaluated average. Therefore, five independent realizations were simulated for all the cases shown in table 1.

Figure 5.

Convergence characteristics of $k^{(f)}/\bar{E}^{(f)}$ with number of realizations for flow past random arrays of spheres at ${\it\phi}=0.4$ and $\mathit{Re}_{m}=100$ . The error bars denote 95 % confidence intervals in the estimation of the average $k^{(f)}$ . The filled symbol represents the level of $k^{(f)}$ corresponding to five realizations.

Besides convergence with the grid resolution and number of realizations, it is also important to check whether the box size is adequate or not. The box size is deemed adequate if the two-point correlation functions in the fluid phase decay to zero within the box length. To check this, the two-point velocity correlation function for the streamwise component of the velocities has been computed for the highest Reynolds number simulated. The fluid-phase velocity autocorrelation ${\it\rho}_{u_{\Vert }}\left(r\right)$ is defined as

where $r=|\boldsymbol{r}|$ and $u_{\Vert }^{\prime \prime (f)}$ is the streamwise component of the gas-phase velocity fluctuation. Figure 6(a) shows the convergence of the fluid-phase velocity autocorrelation function with the grid resolution as well as the box size for a random configuration of particles at a solid volume fraction of 0.2 and a Reynolds number of 20. The autocorrelation function has also been computed for the highest Reynolds number that we simulated and is shown in figure 6(b). These results clearly indicate that the numerical parameters used in our simulation are adequate to obtain numerically converged results for $k^{(f)}$ . We now present the results obtained from PUReIBM simulations of flow past monodisperse fixed particle assemblies.

Figure 6.

(a) Convergence of the fluid-phase velocity autocorrelation function with the grid resolution as well as the box size for a random configuration of particles at a solid volume fraction of 0.2 and a Reynolds number of 20. (b) Decay of the fluid velocity autocorrelation function obtained from PUReIBM simulation of steady flow past a random configuration of spheres at a solid volume fraction of 0.2 and mean slip Reynolds numbers of 20 (squares) and 300 (triangles). In these simulations $L/d_{p}$ ratios of 6 and 4.5 are used for Reynolds numbers of 20 and 300 respectively.

6. Gas-phase velocity fluctuations in steady flow past fixed particle assemblies

While the evolution of the gas-phase velocity fluctuations due to particles has been studied for zero mean slip (Cate et al. Reference Cate, Derksen, Portela and van den Akker2004; Zhang & Prosperetti Reference Zhang and Prosperetti2005; Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2010), the mean slip velocity between the phases is non-zero in fluidized beds. Xu & Subramaniam (Reference Xu and Subramaniam2010) studied the interaction of a turbulent upstream flow with a fixed particle assembly in an attempt to reproduce the experimental findings of Moran & Glicksman (Reference Moran and Glicksman2003). In their study the upstream flow is initialized with isotropic turbulence and drawn through a uniform configuration of fixed particles by imposing a mean constant pressure gradient. This flow is inhomogeneous in the mean flow direction, but gas-phase velocity fluctuations reach a constant value within three particle diameters of entering the bed. The significant result of this study is the enhancement of turbulence as it interacts with solid spheres of the fixed bed.

In the current study we simulate initially turbulent flow in fixed particle assemblies to quantify the relative magnitude of gas-phase velocity fluctuations arising from turbulent and pseudo-turbulent sources. It should be noted that the initial isotropic turbulence is not sustained by artificial forcing and decays in time. The particle configurations correspond exactly to the uniformly distributed particle configurations of Xu & Subramaniam (Reference Xu and Subramaniam2010). These configurations are uniform distributions of non-overlapping spheres generated using the Matérn hard-core point process (Stoyan, Kendall & Mecke Reference Stoyan, Kendall and Mecke1995). This process is essentially a Poisson point process for particle centres from which overlapping spheres have been removed using an approach called dependent thinning. To address the level of turbulent and pseudo-turbulent gas-phase velocity fluctuations we consider three types of simulations:

1. (a) Case I initialized with a non-turbulent uniform laminar flow.

2. (b) Case II initialized with an isotropic turbulent flow.

3. (c) Case III initialized with the steady solution of Case I to which the initial isotropic turbulence of Case II is added.

The simulations start by imposing a mean pressure gradient along the mean flow such that the Reynolds number reaches the desired value. The mean slip Reynolds number is 50 with four realizations for all cases. Detailed information on the flow parameters is provided in table 2.

Table 2.

Numerical and physical parameters of the turbulent/pseudo-turbulent simulations: solid-phase volume fraction ${\it\phi}$ , mean slip Reynolds number $Re_{m}$ , Taylor microscale turbulent Reynolds number $Re_{{\it\lambda}}$ , particle diameter to dissipative length scale ratio $d_{p}/{\it\eta}$ , particle diameter in grid units $d_{p}/{\rm\Delta}x$ , computational box length to particle diameter ratio $\mathscr{L}/d_{p}$ and number of particles $N_{p}$ .

Case I is simply PR-DNS of a uniformly distributed fixed particle assembly initialized with non-turbulent uniform laminar flow, and the corresponding steady value of $k^{(f)}$ is shown as circles in figure 7. The initial isotropic turbulence in Case II is generated by the method described by Rogallo (Reference Rogallo1981) for a Taylor microscale turbulent Reynolds number $Re_{{\it\lambda}}$ of 12 with the energy spectrum function given by Pope (Reference Pope2000). The corresponding Kolmogorov length scale is selected such that the ratio $d_{p}/{\it\eta}$ is 5. Although the $Re_{{\it\lambda}}$ obtained from scaling analysis in experimental studies could be relatively high (the turbulent Reynolds number $Re_{{\it\lambda}}$ is approximately 50 in the work of Moran & Glicksman Reference Moran and Glicksman2003), we are naturally limited by computational resources to lower initial $Re_{{\it\lambda}}$ because particle-resolved simulations for $d_{p}/{\it\eta}>1$ impose additional computational overhead. This is because the mean slip velocity generates flow structures that need to be resolved in PR-DNS (for estimates of computational cost, see for instance Xu & Subramaniam Reference Xu and Subramaniam2010).

Figure 7.

Evolution of $k^{(f)}$ normalized by the mean flow energy $\bar{E}^{(f)}$ versus time normalized by the characteristic flow time scale $d_{p}/|\langle \boldsymbol{W}\rangle |$ . The symbols ▿ correspond to the simulation initialized with an isotropic turbulence. The symbols ▵ represent simulation of the perturbed velocity field. The symbol ○ shows the steady value of the gas-phase pseudo-turbulence.

Figure 7 also indicates that for Case II, the level of $k^{(f)}$ starts from the specified initial isotropic turbulence level and increases to the steady-state pseudo-turbulent value. The data reveal that the velocity fluctuations arising from the presence of particles in the flow with non-zero mean slip velocity are much higher than the turbulent velocity fluctuations for $Re_{{\it\lambda}}=12$ , and the principal contribution to the gas-phase velocity fluctuations is the pseudo-turbulent part.

Since in this study the configurations and simulation parameters are the same as those chosen by Xu & Subramaniam (Reference Xu and Subramaniam2010), we can compare the $k^{(f)}$ values between these homogeneous PUReIBM simulations and the inflow/outflow simulations of Xu & Subramaniam (Reference Xu and Subramaniam2010). In their study, the flow is inhomogeneous along the flow direction and $k^{(f)}$ is reported along the mean flow direction. The level of $k^{(f)}$ increases from the beginning of the bed; it represents a 100 % increase from the entrance of the particle bed and reaches a constant value after an entrance length of $3d_{p}$ . Our homogeneous simulation results are similar for Case II, where the normalized $k^{(f)}$ starts from an initial value of 0.1 and reaches approximately twice its initial value (0.19) at steady state. This indicates that the enhancement observed by Xu & Subramaniam (Reference Xu and Subramaniam2010) is mainly due to the pseudo-turbulent velocity fluctuations, an observation corroborated by experiment (Moran & Glicksman Reference Moran and Glicksman2003).

Case III is basically used to examine the relaxation of turbulence in particle-laden flows. The evolution of $k^{(f)}$ for Case III in figure 7 shows that the excess amount of $k^{(f)}$ decays and the pseudo-turbulent steady-state level of $k^{(f)}$ is recovered. These results, showing both attenuation and enhancement of turbulence depending on the initial level of velocity fluctuations, indicate that for this choice of ${\it\phi}$ and $Re_{m}$ , initial turbulent motions in the gas phase do not influence the steady value of $k^{(f)}$ , which corresponds to the pseudo-turbulent fluctuations arising from the interaction of particles with the mean flow. Since the mechanisms for generation and dissipation of these contributions to $k^{(f)}$ are different, this provides useful information for the development of predictive multiphase turbulence models in this regime of gas–solid flow. Therefore, in the following sections we focus on quantifying the level of pseudo-turbulence and the corresponding dissipation rate over a wide range of flow parameters.

Figure 8.

Variation of $k^{(f)}$ normalized by $\bar{E}^{(f)}=(1/2)\langle \boldsymbol{u}^{(f)}\rangle \boldsymbol{\cdot }\langle \boldsymbol{u}^{(f)}\rangle$ with $\mathit{Re}_{m}$ and ${\it\phi}$ . (a) The behaviour of $k^{(f)}/\bar{E}^{(f)}$ with ${\it\phi}$ for different mean slip Reynolds numbers.(b) The behaviour of $k^{(f)}/\bar{E}^{(f)}$ with $\mathit{Re}_{m}$ for different solid volume fractions.

7. Quantification of PTKE from PR-DNS

We performed PUReIBM DNS of flow past fixed assemblies of monodisperse spheres over a wide range of solid volume fractions ( $0.1\leqslant {\it\phi}\leqslant 0.5$ ) and mean slip Reynolds numbers ( $0.01\leqslant \mathit{Re}_{m}\leqslant 300$ ) to quantify the level of PTKE in terms of ${\it\phi}$ and $\mathit{Re}_{m}$ . All simulations start with the initial condition of uniform fluid velocity. We have shown in the previous section that initializing the simulations with a homogeneous isotropic turbulent velocity field does not affect the steady value of $k^{(f)}$ attained by the system. Therefore, the steady-state value of $k^{(f)}$ obtained in a fixed particle assembly depends only on the solid volume fraction and the mean slip Reynolds number. Figure 8(a) shows the variation of $k^{(f)}/\bar{E}^{(f)}$ with the solid volume fraction for different mean slip Reynolds numbers, while figure 8(b) shows the variation of $k^{(f)}/\bar{E}^{(f)}$ with mean slip Reynolds number for different solid volume fractions. As is evident from figure 8(a), the PTKE normalized by the mean energy in the gas phase increases dramatically with the volume fraction. As shown in figure 8(b), at a given volume fraction $k^{(f)}/\bar{E}^{(f)}$ decreases rapidly with increasing mean slip Reynolds number up to $\mathit{Re}_{m}=50$ , and beyond $\mathit{Re}_{m}=50$ it has a weak power law dependence on $\mathit{Re}_{m}$ . This behaviour is a result of the normalization of $k^{(f)}$ by $\bar{E}^{(f)}$ . It implies that the variance of gas velocity increases approximately as the square of the mean slip Reynolds number. Since the total kinetic energy of the gas increases with increasing mean flow Reynolds number, we expect the strength of gas-phase velocity fluctuations to also increase. Using the data obtained from PUReIBM DNS we found that the following function fits the data for PTKE with an average deviation of 5 %:

As shown in (7.1), the correlation is proposed from simulations in the range $0.1\leqslant {\it\phi}\leqslant 0.5$ and $0.01\leqslant \mathit{Re}_{m}\leqslant 300$ . The value of $k^{(f)}/\bar{E}^{(f)}$ from (7.1) tends to appropriate values in the limit of infinite dilution and creeping flow. In the limiting case of infinite dilution, i.e. ${\it\phi}\rightarrow 0$ , the value of $k^{(f)}/\bar{E}^{(f)}$ is zero. This limiting value is consistent with the fact that in the absence of particles the flow field is uniform. In the Stokes flow regime ( $\mathit{Re}_{m}\rightarrow 0$ ) the value of $k^{(f)}/\bar{E}^{(f)}$ reaches an asymptote and depends only on the solid volume fraction. This behaviour is consistent with the fact that the mean drag (which is shown to be the source of $k^{(f)}$ in the next section) acting on the particles is linear in the Stokes flow regime and thus the normalized quantity $k^{(f)}/\bar{E}^{(f)}$ is independent of the Reynolds number.

Although our study is for homogeneous gas–solid suspensions, the dependence of $k^{(f)}$ on the solid volume fraction and mean slip Reynolds number has implications for transport of the gas-phase Reynolds stress in inhomogeneous flows also. The strong dependence of $k^{(f)}$ on ${\it\phi}$ suggests that the transport of $k^{(f)}$ could be significant in statistically inhomogeneous flows with spatial variation of ${\it\phi}$ . The quantification of the dissipation of PTKE in the gas phase is discussed in the following section.

8. Reynolds stress tensor of pseudo-turbulent gas-phase velocity fluctuations

The gas-phase Reynolds stress is defined as

and its trace represents the level of gas-phase velocity fluctuations $k^{(f)}$ . The presence of particles introduces anisotropy in the gas-phase Reynolds stress that in turn modifies the structure of the carrier phase by extending them along the mean flow (Uhlmann Reference Uhlmann2008; Xu & Subramaniam Reference Xu and Subramaniam2010). The gas-phase anisotropy tensor

is computed along the streamwise (parallel to mean slip) and spanwise (perpendicular to mean slip) directions and presented in figure 9(a,b) respectively. The results indicate that the cross-correlation of velocity fluctuations $b_{\Vert ,\bot }^{(f)}$ is negligible (not shown here), and the normal component in the parallel direction $b_{\Vert ,\Vert }^{(f)}$ is dominant compared with the component in the perpendicular direction $b_{\bot ,\bot }^{(f)}$ . Hence, the state of pseudo-turbulent gas-phase velocity fluctuations is axisymmetric (Pope Reference Pope2000).

Figure 9 indicates that the level of anisotropy increases with $Re_{m}$ from Stokes flow to moderate Reynolds numbers (ranging from 10 to 40) and then smoothly decreases. We hypothesize that this variation in anisotropy is related to the ratio of the length scale associated with gas-phase velocity fluctuations to the local average interparticle spacing. These length scales are quantified using our simulations.

Figure 9.

Deviatoric tensor components are shown for different ${\it\phi}$ and $Re_{m}$ obtained from fixed-bed simulations. The symbols show our PR-DNS data while the solid lines represent the model proposed in (8.6) and (8.7). The shaded region indicates where the maximum anisotropy occurs. (a) Parallel to the mean slip. (b) Perpendicular to the mean slip.

Increase of Reynolds number from low- $Re_{m}$ flow initiates the separation of boundary layer around particles and a standing eddy forms behind spheres with a characteristic length (Pruppacher, Le Clair & Hamielec Reference Pruppacher, Le Clair and Hamielec1970), which is a source of anisotropy in gas-phase Reynolds stress (Mohd-Yusof Reference Mohd-Yusof1996; Uhlmann Reference Uhlmann2008). However, definition of an appropriate length scale for gas-phase velocity fluctuations in a random particle assembly is not straightforward since the wakes are affected by the presence of neighbouring particles. We use the Eulerian two-point correlation of gas-phase velocity fluctuations to define a characteristic length for gas-phase velocity fluctuations $L_{\Vert ,\Vert }$ as follows:

where $R_{\Vert ,\Vert }^{(f)}(r)$ is

the two-point correlation in the parallel direction. This length scale of gas-phase velocity fluctuations, shown in figure 10, decreases with $Re_{m}$ due to the fact that the flow structures become finer and less correlated to other surrounding structures. In addition, increase of the solid volume fraction makes the gas-phase velocity fluctuations become less spatially correlated due to the influence of nearby particles, and hence the length scale decreases.

We use the radial distribution function $g(r)$ to define a characteristic local interparticle spacing within a neighbourhood of a test particle. It should be noted that other studies used a similar approach to compute the nearest neighbour distance for point particles (Hertz Reference Hertz1909) and finite size particles (Torquato, Lu & Rubinstein Reference Torquato, Lu and Rubinstein1990) in random particle arrangements. The radial distribution function is shown in figure 11 for different volume fractions. The radial distribution is the probability of finding a particle at separation $r$ given that there is a particle at the coordinate origin. By assuming a spherical shell of volume $4{\rm\pi}r^{2}\;{\it\delta}r$ at separation $r$ , the number of particles in the shell is $2{\rm\pi}N_{p}\;n\;g(r)r^{2}\;{\it\delta}r$ . Thus, we compute the local interparticle spacing of the suspension using the expression

which is a weighted average of distances among particles in a spherical shell between the minimum separation $d_{p}$ and $R$ . In (8.5), $l$ is the surface to surface distance, defined as $r-d_{p}$ , and $f(l)$ is a weight function. We choose $R$ as the second peak of $g(r)$ , which includes all particles in the local neighbourhood of the test particle. The weight function is assumed to have the form $1/l^{p}$ with $p\geqslant 0$ . It should be noted that the choice of $p=0$ causes $L_{int}$ to be equally weighted by neighbouring particles, while choosing higher values for $p$ preferentially weights the most proximate particles. We select $p=1$ , a value that takes into account the importance of all particles in the region between $d_{p}$ and $R$ in a manner where proximate particles contribute more to the local interparticle spacing.

Figure 10.

The symbols show the length scale of gas-phase velocity fluctuations $L_{\Vert ,\Vert }$ , while the straight lines show the corresponding local interparticle spacing $L_{int}$ . The shaded region shows where the maximum anisotropy occurs. In (a)  ${\it\phi}=0.1$ (▫), ${\it\phi}=0.2$ (▵). In (b)  ${\it\phi}=0.3$ (♢), ${\it\phi}=0.4$ (▿).

Figure 11.

Radial distribution function estimated from 100 realizations for different  ${\it\phi}$ .

The results indicate that the local interparticle spacings (shown as straight lines in figure 10) intersect the length scales of gas-phase velocity fluctuations at the Reynolds numbers where the anisotropy starts to decrease (figure 9). At low Reynolds numbers, flow structures formed behind particles elongate with $Re_{m}$ , which gives rise to the increase of anisotropy. At moderate Reynolds numbers ( $10\leqslant Re_{m}\leqslant 40$ ), the wakes become as large as the gaps among particles and are broken up due to interaction with neighbouring spheres. Thus, the magnitude of anisotropy is determined by the ratio of the gas-phase velocity fluctuation length scale to the local interparticle spacing. The breakup of elongated structures redistributes the energy of the fluctuating velocities among Reynolds stress components and decreases the anisotropy.

It is also observed that increase of the solid volume fraction causes an attenuation in the level of anisotropy. Increase of the solid volume fraction lessens the local interparticle gaps and does not allow the formation of distinct wake structures (Reddy et al. Reference Reddy, Sathe, Joshi, Nandakumar and Evans2013), leading to a decrease in anisotropy. This is confirmed in figure 10, since the length scale of gas-phase velocity fluctuations decreases with the volume fraction. Breault, Guenther & Shadle (Reference Breault, Guenther and Shadle2008) also reported that the axial solid dispersion coefficient of particles in a particle-laden flow decreases with the volume fraction, which is interpreted as decrease of $R_{\Vert ,\Vert }$ while $R_{\bot ,\bot }$ increases, and our data are in agreement with this trend observed in experiments.

Using the gas-phase anisotropy data obtained from our PR-DNS and considering the aforementioned physics explaining the behaviour of gas-phase anisotropy, we found that the following function fits the data for the gas-phase anisotropy:

which is plotted in figure 9 with solid lines. The model constants in (8.6) are given as $a=0.523$ , $b=0.305$ , $c=0.114$ , $d=3.511$ , $e=1.801$ and $f=0.005$ . This model asymptotically reaches a constant value for a fixed ${\it\phi}$ (not shown here) that guarantees that the normal components of the gas-phase anisotropy remain non-zero. This model and the model proposed for $k^{(f)}$ in (7.1) can be used together to propose an ASM for the pseudo-turbulent gas-phase Reynolds stress in the form of

Comparison of this ASM with PR-DNS data as shown in figure 12 indicates that the model appropriately predicts the non-zero components of the pseudo-turbulent gas-phase Reynolds stress. A noteworthy feature of this model is that in the limit of ${\it\phi}\rightarrow 0$ , the gas-phase Reynolds stress becomes zero and the system reduces to a laminar single-phase flow without any pseudo-turbulent velocity fluctuations.

9. Gas-phase velocity fluctuations in freely evolving suspensions

Although fixed beds are good approximations to particle-laden flows at high particle Stokes number, in reality each sphere moves with an acceleration arising from hydrodynamic and collisional forces. Simulation of freely evolving suspensions enables us to study the effect of physical parameters such as the solid-to-gas density ratio and coefficient of restitution, in addition to the solid volume fraction and Reynolds number that are used to characterize fixed-bed simulations. It should be noted that these freely evolving simulations complement the study of Yin et al. (Reference Yin, Zenk, Mitrano and Hrenya2013) by exploring the $Re_{m}$ axis of the parameter space, whereas their work explores the granular temperature dependence for zero $Re_{m}$ .

Figure 12.

The gas-phase Reynolds stress tensor is shown for different ${\it\phi}$ and $Re_{m}$ . The symbols show our PR-DNS data while the solid lines represent the model proposed in (8.8). (a) Parallel to the mean slip. (b) Perpendicular to the mean slip.

In our freely evolving suspension study, a mean pressure is imposed along the mean flow direction such that a desired mean slip Reynolds number is achieved. All our simulations are performed at a mean slip Reynolds number of $Re_{m}=20$ with five independent realizations for each case. Detailed information on particle configurations is provided in table 3. The influence of the solid volume fraction is studied by simulating two volume fractions ${\it\phi}=0.1$ and 0.2. The simulations are initialized with a uniform mean flow for the gas phase and zero granular temperature for the solid phase, and are carried out until a steady value of $k^{(f)}$ is attained.

Table 3.

The numerical and physical parameters of the freely evolving suspensions.

10. Quantification of the dissipation rate of PTKE

Viscous dissipation can be understood by analysing the evolution equation of the mixture energy in statistically homogeneous gas–solid flow given by (9.5). In this simple flow, the steady $e^{(m)}$ results from a balance of mixture interphase TKE transfer and viscous dissipation, bearing in mind that collisional dissipation is negligible. An expression for the mixture interphase TKE transfer ${\it\Pi}^{(m)}$ can be derived by its representation in terms of the average drag force acting per particle:

In the above expression, $F\left({\it\phi},\mathit{Re}_{m}\right)$ is the normalized average drag force per particle given by

where $\langle \boldsymbol{F}\rangle$ is the average hydrodynamic force per particle and $F_{Stokes}=3{\rm\pi}{\it\mu}^{(f)}d_{p}\left(1-{\it\phi}\right)|\langle \boldsymbol{W}\rangle |$ is the Stokes drag acting on an isolated sphere moving with a slip velocity of $\left(1-{\it\phi}\right)|\langle \boldsymbol{W}\rangle |$ . The expression in (10.1) for the source of $e^{(m)}$ due to interphase transfer of kinetic energy is similar to the one derived by Crowe (Reference Crowe2000). While Crowe (Reference Crowe2000) used the single-sphere drag correlation for $F\left({\it\phi},\mathit{Re}_{m}\right)$ , here we obtain this value directly from the PR-DNS. An accurate correlation for $F\left({\it\phi},\mathit{Re}_{m}\right)$ has been developed using the data obtained from PUReIBM simulations (Tenneti et al. Reference Tenneti, Garg and Subramaniam2011). The drag correlation is summarized below for the sake of completeness. The average normalized drag force acting per particle in flow past a random assembly of monodisperse spheres is given by

where $F_{\mathit{isol}}$ is the drag force acting on an isolated sphere moving in an unbounded medium. The drag on an isolated sphere is taken to be the correlation proposed by Schiller & Naumann (Reference Schiller and Naumann1935). The remaining two terms in (10.3) are given by

The dissipation rate of PTKE can be quantified by noting that the source and sink of $e^{(m)}$ must balance each other at steady state, i.e.

Since the source of $e^{(m)}$ depends only on the average interphase momentum transfer and the mean slip velocity, knowledge of the average drag force acting on the particles is enough to quantify the dissipation. At any solid volume fraction and mean slip Reynolds number, the dissipation rate can be readily inferred from (10.6) and (10.1) and by using a drag correlation (e.g. (10.3)).

It should be noted that these expressions for mixture interphase TKE transfer and viscous dissipation have been developed using a drag model obtained from PR-DNS of uniformly distributed fixed particle assemblies. However, in practical systems where particle clustering is observed and in unsteady and inhomogeneous flows where the particle granular temperature is far from its steady value, these expressions may not be applicable in two-fluid RANS calculations. Nevertheless, these expressions are still valuable for modelling the residual fluid stresses in LES calculations that are able to resolve mesoscale features of the flow such as particle clustering.