2.1. Spatial decomposition of the particle velocity field

In the context of high-inertia particles with response times that are long compared to the characteristic time scale of the turbulence, individual particle trajectories will retain information from previous collisions and interactions with distant turbulent eddies, causing them to deviate from fluid pathlines (Maxey Reference Maxey1987). Velocities of neighbouring particles may therefore be uncorrelated, while ensembles of particles collectively respond to large-scale motions of the flow. The instantaneous velocity of any particle  $i$ , $\boldsymbol{v}_{p}^{(i)}$ , can be decomposed into a spatially correlated contribution $\boldsymbol{u}_{p}$ and an uncorrelated (residual) component ${\bf\delta}\boldsymbol{v}_{p}^{(i)}$ (Février, Simonin & Squires Reference Février, Simonin and Squires2005), i.e.

where $\boldsymbol{x}_{p}^{(i)}$ is the centre position of particle  $i$ . For a single realization of the flow, the total particle-phase fluctuating energy is given by

where $\boldsymbol{v}_{p}^{\prime }=\boldsymbol{v}_{p}-\langle \boldsymbol{v}_{p}\rangle$ is the total fluctuation in particle velocity with the property $\langle \boldsymbol{v}_{p}^{\prime }\rangle =0$ . Note that the averaging operator $\langle (\cdot )\rangle$ is used throughout this work to denote a particle average when applied to a Lagrangian quantity, and a volume average when applied to an Eulerian quantity. Owing to the statistical stationarity and homogeneity of the flows considered in this work, $\langle (\cdot )\rangle$ is a function neither of the spatial coordinate $\boldsymbol{x}$ nor of time $t$ at steady state.

Analogous to Favre averaging in variable-density flows, the PA with respect to the particle phase, denoted by $\langle (\cdot )\rangle _{p}=\langle {\it\alpha}_{p}(\cdot )\rangle /\langle {\it\alpha}_{p}\rangle$ , is useful in multiphase modelling (Fox Reference Fox2014). Note that PA Eulerian quantities are identical to particle-average Lagrangian quantities, e.g.  $\langle \boldsymbol{u}_{p}\rangle _{p}=\langle \boldsymbol{v}_{p}\rangle$ . Fluctuations about the PA velocity are expressed as $\boldsymbol{u}_{p}^{\prime \prime }(\boldsymbol{x},t)=\boldsymbol{u}_{p}(\boldsymbol{x},t)-\langle \boldsymbol{u}_{p}\rangle _{p}$ , with $\langle \boldsymbol{u}_{p}^{\prime \prime }\rangle _{p}=0$ . Using this definition, the PA particle-phase turbulent kinetic energy is defined as

It is important to note that $\boldsymbol{u}_{p}^{\prime \prime }\neq \boldsymbol{u}_{p}^{\prime }$ , and therefore $\langle \boldsymbol{u}_{p}^{\prime \prime }\rangle \neq 0$ in general.

A quantitative measure of the local uncorrelated particle agitation is given by the granular temperature ${\it\Theta}$ , which is defined using the residual component of the instantaneous particle velocity ${\bf\delta}\boldsymbol{v}_{p}^{(i)}$ , such that the PA granular temperature is given by

With these definitions, the total particle-phase fluctuation energy ${\it\kappa}_{p}$ corresponds to the sum of the correlated and uncorrelated granular energy, i.e.

The distinction between $k_{p}$ and $\langle {\it\Theta}\rangle _{p}$ is crucial in turbulence modelling. For instance, in the context of moderately dense particulate flows, $\langle {\it\Theta}\rangle _{p}$ is needed to evaluate the particle-phase viscosity and pressure, which arise due to collisions. Thus, failure to separate these two contributions will lead to a gross overprediction of the collision rate (Fox Reference Fox2014). Moreover, previous works (Hrenya & Sinclair Reference Hrenya and Sinclair1997; Février et al. Reference Février, Simonin and Squires2005) have shown that the dissipation of $k_{p}$ enters as a source term for $\langle {\it\Theta}\rangle _{p}$ . This is analogous to single-phase flow, where dissipation of turbulent kinetic energy leads to viscous heating.

2.2. Mesoscale description of fluid–particle systems

A complete description of each particle and the interstitial fluid phase is provided by the microscale model, which solves the Navier–Stokes equation for the fluid with no-slip boundary conditions applied at the surface of each particle, along with Newton’s equations for the momentum of each individual particle (see e.g. Tenneti & Subramaniam Reference Tenneti and Subramaniam2014). However, its use is limited to relatively small ensembles of particles, typically too small to observe the macroscale turbulence that we wish to capture in this study. Consequently, a mesoscale model can be formulated from the microscopic description by introducing a set of variables that adequately describes the physics of the disperse phase (Subramaniam Reference Subramaniam2000; Fox Reference Fox2007). For systems containing a large number of particles, it is often convenient to define a number-density function $f(t,\boldsymbol{x},\boldsymbol{v}_{p})$ , which depends on time $t$ , position $\boldsymbol{x}$ and particle velocity $\boldsymbol{v}_{p}$ . The stochastic nature of the particle phase can be described by a kinetic equation given by

where $\mathbf{\mathscr{A}}_{p}$ is the acceleration due to fluid–particle interactions, $\boldsymbol{g}=[-g,0,0]^{\text{T}}$ is the acceleration due to gravity, and $\mathbb{C}$ is the collision operator. Owing to the large number of independent variables associated with $f$ in a three-dimensional flow, a direct solution of the kinetic equation is still expensive in most practical flows, and various mathematical approximations have been made to arrive at tractable mesoscale models. In a recent review paper, Fox (Reference Fox2012) summarizes such modelling strategies, in addition to providing a systematic approach for developing LES tools starting from microscale physics.

2.3. Eulerian description of the particle phase

The Eulerian description of the particle phase is formulated in terms of moments of the mesoscopic kinetic model (2.6). The governing equations take the form of Eulerian transport equations with coupling terms between the phases due to mass and momentum exchange. The KT model of granular flows (Jenkins & Savage Reference Jenkins and Savage1983) is now generally accepted as a valid macroscale description of moderately dense particle flows. In this regime, the particle volume fraction is large enough such that particle–particle collisions drive the velocity distribution function towards the equilibrium (Maxwellian) distribution, though not too large that particles remain in sustained contact. Assuming that the particles are frictionless hard spheres of equal density and diameter (i.e. monodisperse) and that collisions are nearly elastic, the conservation of mass and momentum in the presence of a constant-density fluid are given by

and

From (2.16), the transport equation for the particle-phase velocity tensor product can be obtained as

It is typically assumed that the particle-phase Knudsen number is very small so that a Chapman–Enskog expansion can be used to provide constitutive relations for the pressure tensor $\unicode[STIX]{x1D64B}$ (Lun et al. Reference Lun, Savage, Jeffrey and Chepuriny1984). However, for non-equilibrium flows, a transport equation for the pressure tensor is necessary. Defining the second-order moment of the velocity distribution function as $\unicode[STIX]{x1D614}^{2}=\int \boldsymbol{v}_{p}\otimes \boldsymbol{v}_{p}f\,\text{d}\boldsymbol{v}_{p}$ , from (2.6) and (2.17) the transport equation for the central moment ${\it\alpha}_{p}\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D614}^{2}-{\it\alpha}_{p}\boldsymbol{u}_{p}\otimes \boldsymbol{u}_{p}$ can be derived:

In this equation, $\unicode[STIX]{x1D64C}$ is a heat-flux tensor that contains the third-order central moments of the velocity distribution function, and the last term on the right-hand side is the particle–particle collision term that has been closed using the Bhatnagar–Gross–Krook (BGK) approximation (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954) extended to inelastic collisions (Passalacqua et al. Reference Passalacqua, Galvin, Vedula, Hrenya and Fox2011), where ${\it\Theta}=\text{tr}(\unicode[STIX]{x1D64B})/3$ is the granular temperature, $0\leqslant e\leqslant 1$ is the coefficient of restitution and ${\it\Delta}^{\star }$ is the second-order moment of the collisional equilibrium distribution, given by

In the limit of elastic collisions ( $e=1$ ), the collision term on the right-hand side of (2.18) contains only contributions from the deviatoric pressure tensor.

The particle-phase pressure tensor can be decomposed into an isotropic contribution ${\it\Theta}\unicode[STIX]{x1D644}$ and the particle-phase viscous stress tensor defined by ${\bf\sigma}_{p}={\it\Theta}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D64B}$ . By taking one-third of the trace of (2.18), the equation for the granular temperature is found to be

where $\boldsymbol{q}$ is the granular energy flux. The second term on the right-hand side of (2.20) represents the loss of granular energy due to drag, and the last term is dissipation due to inelastic collisions. The first term on the right-hand side of (2.20) can be rewritten using the decomposition $\unicode[STIX]{x1D64B}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{p}={\it\Theta}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{p}-{\bf\sigma}_{p}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{p}$ . The contribution $-{\it\alpha}_{p}{\it\Theta}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{p}$ represents compressive heating of the particle phase, while the second contribution represents production due to viscous dissipation of mean kinetic energy. As shown in figure 2, taken from Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2014a ), compressive heating upstream of falling clusters makes an important contribution to the granular temperature balance in CIT.

Figure 2.

Instantaneous field of granular temperature generated using data from Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2014a ) corresponding to $\mathit{Re}_{p}=1$ in this work (refer to table 1). Full domain showing iso-contours of ${\it\alpha}_{p}=3\langle {\it\alpha}_{p}\rangle$ (a) and zoomed-in snapshot showing particle velocity vectors in the vicinity of a cluster (b). High values of the granular temperature correspond to regions in front of falling clusters where the particle-phase velocity $\boldsymbol{u}_{p}$ is strongly compressed in the vertical direction. Particle size and mesh are shown to scale.

2.4. Reynolds-average equations of motion

The wide range of length and time scales associated with turbulent multiphase flows makes a direct solution of the transport equations presented thus far intractable for many environmental and industrial applications. It is therefore necessary to obtain from first principles a set of macroscopic equations that can correctly capture the statistical properties of the flow. Here we present the exact RA transport equations recently derived by Fox (Reference Fox2014), and extend the analysis to include the full Reynolds stress and granular stress equations. The relevant unclosed terms that are retained in the context of homogeneous gravity-driven flows are then presented and discussed.

3.1. Configuration and simulation parameters

The volume-filtered fluid-phase equations presented in § 2.2.2 and the Lagrangian description of the particle phase given in § 2.2.1 are implemented in the framework of NGA (Desjardins et al. Reference Desjardins, Blanquart, Balarac and Pitsch2008), a second-order, fully conservative computational fluid dynamics code tailored for turbulent flow computations. To assess the unclosed terms that appear in the RA equations from § 2.4, a range of simulations are conducted in the context of fully developed gravity-driven CIT. The independent dimensionless parameters that characterize homogeneous fluid–particle flows include the density ratio ${\it\rho}_{p}/{\it\rho}_{f}$ , the mean volume fraction $\langle {\it\alpha}_{p}\rangle$ and the mean particle Reynolds number defined as

where the characteristic velocity $\mathscr{V}={\it\tau}_{p}g$ is the particle settling velocity in an undisturbed flow field when ${\it\rho}_{p}\gg {\it\rho}_{f}$ . Note that combining these non-dimensional numbers yields a Froude number $\mathit{Fr}={\it\rho}_{p}\mathit{Re}_{p}/(18{\it\rho}_{f})$ and the mean mass loading  ${\it\varphi}$ . To obtain an a priori measure of the mesoscale features that arise due to the coupling between the phases, previous studies have introduced a characteristic length scale $\mathscr{L}={\it\tau}_{p}^{2}g$ in order to assess an appropriate domain size such that the results are unaffected by the periodic boundary conditions (see e.g. Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Igci et al. Reference Igci, Andrews, Sundaresan, Pannala and O’Brien2008; Ozel et al. Reference Ozel, Fede and Simonin2013). This length scale can be used to formulate a ‘mesoscale’ Reynolds number to characterize the fluid-phase inertia, i.e.

The simulations presented in this work are solved on a triply periodic domain with gravity directed in the negative $x_{1}$ direction. The flows are made statistically stationary by a source term added to the pressure gradient in the fluid-phase momentum equation (2.10), which is computed dynamically in order to prevent the development of a net mass flow rate (i.e.  $\langle \boldsymbol{u}_{f}\rangle _{f}=0$ ). The particles are monodisperse and initially uniformly distributed throughout the domain of dimensions $3584d_{p}\times 896d_{p}\times 896d_{p}$ with a mesh size of $2048\times 512\times 512$ , corresponding to a uniform grid spacing of ${\rm\Delta}x=1.75d_{p}$ . The unclosed terms that arise in the RA transport equations are investigated for a range of $\mathit{Re}_{p}$ by keeping all parameters constant and modifying $g$ . A list of each simulation case and their corresponding parameters is given in table 1. Note that the Archimedes number, $\mathit{Ar}=({\it\rho}_{p}/{\it\rho}_{f}-1)d_{p}^{3}g/{\it\nu}^{2}$ , is commonly used in the literature to characterize gravity-driven fluid–particle flows (e.g. Noymer & Glicksman Reference Noymer and Glicksman2000; Ozel et al. Reference Ozel, Fede and Simonin2013), and is related to our definition of the particle Reynolds number by $\mathit{Ar}=18\mathit{Re}_{p}$ , or equivalently to the square of the Galileo number defined in the introduction.

Table 1.

Fluid–particle parameters used in CIT simulations.

3.2. Implementation of the mesoscale equations

In this work, the Navier–Stokes equations are solved on a staggered grid with second-order spatial accuracy for both the convective and viscous terms. Time advancement is accomplished using the second-order accurate semi-implicit Crank–Nicolson scheme of Pierce (Reference Pierce2001). Momentum exchange between the phases is accounted for by the drag term $\mathbf{\mathscr{A}}$ given in (2.13). To account for unresolved fluid velocity fluctuations generated by the particles, a model derived by Gibilaro et al. (Reference Gibilaro, Gallucci, Di Felice and Pagliai2007) is used for the effective viscosity, given by

We briefly note that contributions from ${\it\nu}_{f}^{\star }$ were found to be negligible compared to ${\it\nu}_{f}$ in the simulation cases considered in this work.

The particles are distributed among the processors based on the underlying domain decomposition of the fluid phase. For each particle, its position and velocity are solved using a second-order Runge–Kutta scheme. The collision force $\boldsymbol{F}_{c}$ that appears in (2.7) is modelled using a soft-sphere approach originally proposed by Cundall & Strack (Reference Cundall and Strack1979). When two particles come into contact, a repulsive force of particle $j$ on particle $i$ is created as

where $k$ and ${\it\eta}$ are the spring stiffness and the damping parameter, respectively, $\boldsymbol{n}_{ij}$ is the unit normal vector from particle $i$ to particle $j$ , ${\it\delta}_{ij}$ is the overlap between the particles, and the normal relative velocity between particles $i$ and $j$ is given by

The collision force is computed when the distance between the centres of the particles $d_{ij}$ is less than the particle diameter $d_{p}$ and a force range ${\it\lambda}$ . The force range ${\it\lambda}$ is adjusted dynamically in order to recover the correct random close-packing limit when particles are in sustained contact, while providing a robust model for high-speed collisions (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013a ). A model for the damping parameter is given by

where $m_{p}={\rm\pi}{\it\rho}_{p}d_{p}^{3}/6$ is the particle mass. In this work, the restitution coefficient is taken to be $e=0.9$ . In order to limit inter-particle overlap during collisions, the simulation time step, ${\rm\Delta}t$ , is chosen such that particles cannot move more than one-tenth of their diameter per time step. The spring stiffness is related to the collision time, ${\it\tau}_{c}$ , according to

To properly resolve the collisions without requiring an excessively small time step, ${\it\tau}_{c}$ is chosen to be an order of magnitude larger than ${\rm\Delta}t$ . Once each individual collision force is computed, the full collision force that particle $i$ experiences can be expressed as a sum of collisions with all other particles $j$ , i.e.

Details on the parallel implementation of the Lagrangian particle tracking framework can be found in Capecelatro & Desjardins (Reference Capecelatro and Desjardins2013a ).

3.3. Two-way coupling

Coupling between the fluid phase and particles appears in the form of the volume fraction ${\it\alpha}_{p}$ and interphase exchange term $\mathbf{\mathscr{A}}$ . These terms are first computed at the location of each particle using information from the fluid, and are then transferred to the Eulerian mesh. To interpolate the fluid variables to the particle location, a second-order trilinear interpolation scheme is used. To extrapolate the particle data back to the Eulerian mesh, the volume-filtering approach used in the mathematical formulation of the fluid-phase mesoscale equations is applied. We begin by defining a filtering kernel $G$ with a characteristic length ${\it\delta}_{f}$ , such that $G(r)>0$ decreases monotonically with increasing $r$ , and is normalized such that it integrates to unity. Given a quantity $A^{(i)}(t)$ located at the centre of particle $i$ at time $t$ , and assuming $G$ does not vary significantly over the volume of the particle (i.e.  ${\it\delta}_{f}\gg d_{p}$ ), its Eulerian projection is given by

where $N_{p}$ is the total number of particles in a single realization of the flow and $V_{p}={\rm\pi}d_{p}^{3}/6$ is the particle volume. This expression replaces the discontinuous Lagrangian data with an Eulerian field that is a smooth function of the spatial coordinate  $\boldsymbol{x}$ .

Note that (3.9) will only yield useful information if the spatial variations in the particle field can be decomposed into contributions on a scale comparable with the particle spacing and on a much larger scale corresponding to mesoscopic features in the flow (e.g. clusters), provided the filter size ${\it\delta}_{f}$ is within these scales. For ratios of ${\rm\Delta}x/d_{p}\approx 1$ , a brute-force implementation of (3.9) would require looping through a large number of cells for each particle, making this operation prohibitively expensive. Therefore, the filtering procedure is solved in two steps (Capecelatro & Desjardins Reference Capecelatro and Desjardins2013a ). First, the particle data are transferred to the nearest-neighbouring cells via trilinear extrapolation. The data are then diffused such that the final width of the filtering kernel is independent of the mesh size. In this work, $G$ is taken to be Gaussian with a characteristic length scale ${\it\delta}_{f}=8d_{p}$ , defined as the full width at half the height of the kernel. This value of ${\it\delta}_{f}$ was found to provide accurate solutions of the spatial distribution of the particle field in our previous work (Capecelatro et al. Reference Capecelatro, Desjardins and Fox2014a ). To keep the cost low and ensure unconditional stability, the diffusion process is solved in a single implicit step by utilizing the approximate factorization scheme of Briley & McDonald (Reference Briley and McDonald1977).

3.4. Numerical approach for extracting multiphase statistics

We postulate that evaluating the multiphase statistics presented in § 2.4, in particular $k_{p}$ and $\langle {\it\Theta}\rangle _{p}$ , can be performed by introducing a separation of length scales when sampling the particle data. However, this requires validation, which is the purpose of this section. To accomplish this separation of length scales, we employ the spatial filtering procedure given by (3.9). The spatial filter is used in the first step of a two-step process for obtaining relevant statistics. The first step computes instantaneous Eulerian fields from Lagrangian particle data. The second step computes RA statistics from these Eulerian fields. In this section, we describe how the filter is defined such that the Eulerian turbulence statistics are computed accurately.

4.1. Visualizations of instantaneous fully developed CIT

The simulations are run until the initial transient is complete and the flows reach a statistically stationary state. Owing to the non-trivial coupling between the phases, the transients were observed to persist longer than $80{\it\tau}_{p}$ . Once a statistical steady state has been reached, results are measured at each computational time step ( ${\rm\Delta}t=3.25\times 10^{-4}\,{\it\tau}_{p}$ , $5.2\times 10^{-4}\,{\it\tau}_{p}$ and $8.0\times 10^{-4}\,{\it\tau}_{p}$ for $\mathit{Re}_{p}=0.25$ , 0.5 and 1.0, respectively), over a duration of approximately $10{-}15\,{\it\tau}_{p}$ . As shown in figure 9, the resulting particle field is highly non-uniform, with significant segregation in volume fraction. The corresponding fields of fluid velocity are shown in figure 10. The fluid is observed to be entrained by the clusters, leading to strong vertical velocities in dilute regions of the flow, indicative of jet bypassing. Furthermore, the magnitude of the local fluid-phase velocity, normalized by their respective terminal velocities, is seen to decrease with increasing  $\mathit{Re}_{p}$ .

Figure 9.

Instantaneous fields of particle volume fraction in a two-dimensional plane of the fully developed CIT with (a)  $\mathit{Re}_{p}=0.25$ , (b)  $\mathit{Re}_{p}=0.5$ and (c)  $\mathit{Re}_{p}=1$ .

Figure 10.

Instantaneous fields of fluid velocity in a two-dimensional plane of the fully developed CIT with (a)  $\mathit{Re}_{p}=0.25$ , (b)  $\mathit{Re}_{p}=0.5$ and (c)  $\mathit{Re}_{p}=1$ . Iso-countours of ${\it\alpha}_{p}=3\langle {\it\alpha}_{p}\rangle$ are shown in black.

Table 2.

Relative contributions to kinetic energy production and dissipation.

4.2. Overall kinetic energy balance

The unclosed terms appearing in the RA equations are evaluated via the adaptive spatial filter introduced in § 3.4, and each equation is analysed in detail in § 4.4. Based on the magnitude of the terms in the RA equations, we can identify the overall kinetic energy balance for CIT shown in figure 11. In this balance, the drag production (DP) is proportional to $\langle u_{f,1}^{\prime \prime \prime }\rangle _{p}$ , which is zero in the absence of clusters. As in single-phase turbulence, the viscous dissipation (VD) results from the resolved small-scale velocity fluctuations in the fluid phase. In contrast, the drag exchange terms represent two separate phenomena: (i) viscous dissipation of unresolved fluid velocity fluctuations ( $\text{DE}_{2}$ ), e.g. in the viscous boundary layers around individual particles; and (ii) energy transferred to the particles at the fluid–particle interface ( $\text{DE}_{1}$ ). As described by Sundaram & Collins (Reference Sundaram and Collins1999), the unresolved dissipation arises in the point-particle model due to the difference between the fluid velocity at the particle location and the particle velocity. Denoting the sum of the two contributions as $\text{DE}_{t}$ , in table 2 we observe that the relative contribution of the interphase energy exchange ( $\text{DE}_{1}/\text{DE}_{t}$ ) is approximately 22 % for all $\mathit{Re}_{p}$ . Furthermore, the ratio of the resolved to unresolved viscous dissipation is less than 6 % and decreases with increasing $\mathit{Re}_{p}$ . It is interesting to note in table 2 that the fraction of the total kinetic energy production going to DP is largest for the smallest $\mathit{Re}_{p}$ . Presumably this fraction would approach unity as $\mathit{Re}_{p}$ approaches zero if the clusters were fully developed. Numerically, this conjecture would be difficult to verify starting from uniformly distributed particles because the transient to fully developed clusters increases with decreasing  $\mathit{Re}_{p}$ .

Figure 11.

Schematic of the kinetic energy balance in CIT. Kinetic energy is produced by particles falling due to gravity. In the absence of clustering, mean kinetic energy production is balanced by mean drag ( ${\it\phi}u_{p}^{2}/{\it\tau}_{p}$ with $u_{p}=\mathscr{V}$ ), which generates viscous heating of the fluid (i.e. increases the thermodynamic temperature  $T_{f}$ ). In fully developed CIT, drag production (DP) generates fluid-phase turbulent kinetic energy ( $k_{f}$ ), which is dissipated to heat the fluid directly by viscous dissipation (VD) and indirectly by drag exchange ( $\text{DE}_{2}$ ). The remaining fraction of the fluid-phase turbulent kinetic energy is transferred to the particle phase by drag exchange ( $\text{DE}_{1}$ ), and is subsequently dissipated to heat the fluid through the fluctuating drag ( ${\it\kappa}_{p}/{\it\tau}_{p}$ ).

As will be shown in the following section, owing to the formation of clusters, in fully developed CIT the average falling velocity of particles ( $|\langle u_{p,1}\rangle _{p}|$ ) increases significantly relative to the particle falling velocity without clustering ( $\mathscr{V}$ ). Some of this additional particle-phase kinetic energy is dissipated by mean drag, while the remainder is transferred to fluid-phase turbulent kinetic energy through drag production. In the fluid phase, nearly 80 % of the turbulent kinetic energy is dissipated by unresolved viscous dissipation (i.e.  $\text{DE}_{2}$ ) and the rest is transferred to particle-phase fluctuation energy. As we shall see in § 4.4, only a very small fraction of the particle-phase fluctuation energy is dissipated due to inelastic particle–particle collisions.

4.3. One-point statistics

Having established the overall kinetic energy balance for fully developed CIT, we now turn our attention to specific one-point statistics for the fluid and particle phases. In table 3, we present statistics for the fluid phase for each value of $\mathit{Re}_{p}$ . First, note that, although the scaled value of $k_{f}$ decreases with $\mathit{Re}_{p}$ , the anisotropy between the vertical and horizontal components is nearly constant. Second, the fluid velocity statistics exhibit non-Gaussian behaviour. The vertical velocity component has significant positive skewness, while its kurtosis is close to the Gaussian value of 3. The horizontal component has negligible skewness (as would be expected due to symmetry), but with a large kurtosis. The non-Gaussian behaviour of the one-point fluid velocity statistics is confirmed from the one-point p.d.f.s shown in figure 12. Details on the formulation of the one-point p.d.f.s can be found in appendix A. Finally, note that these p.d.f.s exhibit reasonably good collapse for all  $\mathit{Re}_{p}$ .

Figure 12.

The p.d.f.s of the (a) streamwise and (b) spanwise components of PA fluid-phase velocity from fully developed CIT with $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ) compared to a Gaussian distribution (▵).

Table 3.

Fluid-phase velocity statistics.

The one-point statistics for the particle-phase volume fraction are given in table 4. The corresponding one-point p.d.f.s are shown in figure 13. For a homogeneous distribution of particles, devoid of any clustering, the p.d.f. is given by the discrete Poisson distribution (Squires & Eaton Reference Squires and Eaton1991; Eaton & Fessler Reference Eaton and Fessler1994), which takes the form

where $N$ is the observed number of particles in a given sample, and $\overline{N}$ is the average particle number. In figure 13, the volume-fraction p.d.f.s display a higher probability of local regions empty of particles, as well as local regions of higher solid fraction, compared to (4.1). Although the characteristic cluster size appears to increase with $\mathit{Re}_{p}$ in figure 9, the distribution of volume fraction is shown to be unaffected by the Reynolds number. Using the mean and variance particle volume fraction extracted from $\mathit{Re}_{p}=1$ , the forms of the p.d.f.s are seen to closely resemble a lognormal distribution, indicating a potential opportunity for future modelling efforts (see appendix B).

Figure 13.

The p.d.f. of particle volume fraction from single realizations of the fully developed CIT with $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ), corresponding to lognormal (○) and Poisson (▵) distributions.

In table 5, we present statistics for the particle phase for all values of $\mathit{Re}_{p}$ . First, we observe that the fall velocity $|\langle u_{p,1}\rangle _{p}|$ is more than twice as large as the reference velocity $\mathscr{V}$ due to the formation of clusters. Second, note that, as for the fluid, the scaled value of ${\it\kappa}_{p}$ decreases with $\mathit{Re}_{p}$ , but the anisotropy between the vertical and horizontal components is nearly constant. The same is true for $k_{p}$ and $\langle {\it\Theta}\rangle _{p}$ . Note also that, while the anisotropy in the correlated components (e.g.  $\langle u_{p,1}^{\prime \prime 2}\rangle _{p}$ ) is larger than in the uncorrelated components (e.g.  $\langle \unicode[STIX]{x1D617}_{11}\rangle _{p}$ ), the latter is too large to be approximated as isotropic. Next, the one-point p.d.f.s of the correlated components of the particle velocity statistics are shown in figure 14, displaying a reasonably good collapse for all $\mathit{Re}_{p}$ . While the vertical component is observed to be positively skewed, the overall behaviour of the correlated particle-phase velocity statistics is close to Gaussian. In contrast, the velocity statistics for the uncorrelated components of the particle velocity are significantly non-Gaussian, as seen by the large third-order moments in table 5, as well as the one-point p.d.f.s shown in figure 15. Sundaram & Collins (Reference Sundaram and Collins1997) report a non-Gaussian p.d.f. for the relative particle velocity in isotropic, particle-laden turbulence. This observation is consistent with figure 15 because the relative particle velocity will be mainly due to the spatially uncorrelated particle velocity field represented by ${\bf\delta}\boldsymbol{v}_{p}$ . It is evident from the exponential tails of the p.d.f.s that the uncorrelated components retain a strong imprint of their production terms from the correlated components that is not erased by particle–particle collisions.

Table 4.

Particle-phase volume-fraction statistics.

Table 5.

Particle-phase velocity statistics.

Figure 14.

The p.d.f.s of the (a) streamwise and (b) spanwise components of PA particle-phase velocity from fully developed CIT with $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ) compared to a Gaussian distribution (▵).

Figure 15.

The p.d.f.s of the (a) streamwise and (b) spanwise components of velocity residual ${\bf\delta}\boldsymbol{v}_{p}$ used in defining the particle pressure tensor $\unicode[STIX]{x1D64B}$ and granular temperature ${\it\Theta}$ , with $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ) compared to a Gaussian distribution (▵).

As described in appendix A, it is straightforward to compute statistics for the fluid phase seen by the particles from the Eulerian fields. Thus, in table 6, selected fluid and fluid–particle statistics seen by the particles are presented. First, the normalized mean fluid velocity seen by the particles $|\langle u_{f,1}\rangle _{p}|/\mathscr{V}$ can be seen to decrease with $\mathit{Re}_{p}$ . As discussed earlier, a non-zero value for $\langle u_{f,1}\rangle _{p}$ is a direct result of particle clustering and the fact that clusters entrain the fluid in their vicinity due to momentum coupling. Second, note that the fluid-phase turbulence kinetic energy seen by the particles $k_{f@p}$ is nearly the same as $k_{f}$ ; however, its components are more anisotropic due to the drag production term. Third, the fluid–particle kinetic energy $k_{fp}$ is approximately equal to $0.87\sqrt{k_{f}k_{p}}$ for all $\mathit{Re}_{p}$ , and the anisotropy of its components is nearly the same as for $k_{f}$ and $k_{p}$ . Finally, the fluid velocity statistics seen by the particles are super-Gaussian, which is confirmed by the one-point p.d.f.s shown in figure 16.

Table 6.

Fluid-phase turbulence statistics seen by particles.

Figure 16.

The p.d.f.s of the (a) streamwise and (b) spanwise components of fluid velocity seen by the particles from fully developed CIT with $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ) compared to a Gaussian distribution (▵).

In summary, in this section we have presented one-point statistics for the principal quantities of interest in fully developed CIT. From these statistics, we see that the momentum coupling between the two phases leads to significant differences from the behaviour observed for very dilute systems with one-way coupling (Minier & Peirano Reference Minier and Peirano2001). Most importantly, owing to momentum coupling, the particle phase spontaneously organizes into clusters (or, equivalently, volume-fraction fluctuations with a near lognormal p.d.f.) that significantly increases the mean particle velocity relative to an isolated particle in an undisturbed flow. In turn, the increased mean particle velocity generates a drag production term for turbulent kinetic energy in the fluid phase that is highly anisotropic. Furthermore, the fluid-phase velocity p.d.f.s are non-Gaussian. Owing to the compressibility of the particle phase, the uncorrelated components of the particle-phase velocity statistics are highly non-Gaussian. This situation can be contrasted with the one-way coupling case where in the homogeneous limit all of the velocity statistics are nearly Gaussian (Minier & Peirano Reference Minier and Peirano2001). These differences would suggest that a predictive p.d.f. model for CIT must explicitly account for the compressibility of the particle phase. Finally, it is noteworthy that the uncorrelated particle velocity statistics (i.e.  $\unicode[STIX]{x1D64B}$ ) are significantly anisotropic in our simulations. At least for systems with moderately low mean particle volume fractions for which particle–particle collisions are too infrequent to ensure that $\unicode[STIX]{x1D64B}$ is nearly isotropic, this fact comprises the effectiveness of traditional kinetic-theory closures for $\unicode[STIX]{x1D64B}$ wherein a constitutive equation is used to represent  ${\bf\sigma}_{p}$ .

4.4. Turbulence budgets

Having provided an overall analysis of the kinetic energy budget for CIT in § 4.2, we now turn to a detailed analysis of the individual RA balance equations derived in § 2.4. In particular, we compute each term in the steady-state balance equations for homogeneous gravity-driven flow and use the notation introduced in § 2.4 to denote each term. From tables 7–11, it can be noticed that, after time averaging the various RA terms, some of the equations remain unbalanced. Defining an error as the residual normalized by the largest term appearing in the equation, it is seen that the RA balance equations for $\langle u_{p,1}\rangle _{p}$ and $k_{f}$ exhibit errors below 11 %, while the particle-phase energy budgets ( $k_{p}$ and $\langle {\it\Theta}\rangle _{p}$ ) are observed to have errors as large as 50 %. The errors associated with the particle-phase balance equations are attributed to inconsistencies with the Lagrangian reconstruction of the various Eulerian terms, in particular the terms involving gradients of the particle-phase Eulerian fields. Meanwhile, the relative importance of the key terms appearing in the RA equations are captured and provide useful insight on the energy exchange in CIT.

4.5. Conditional statistics

As shown earlier, the particle-phase volume fraction has a one-point p.d.f. close to lognormal and exhibits large fluctuations about its mean value. In order to account for these fluctuations, one possible approach is to consider one-point statistics conditioned on ${\it\alpha}_{p}$ . Formally, the one-point joint p.d.f. of ${\it\alpha}_{p}$ and any other flow quantity $Z$ can be written as $f(Z,{\it\alpha}_{p})=f_{1}(Z|{\it\alpha}_{p})f_{2}({\it\alpha}_{p})$ , where $f_{1}(Z|{\it\alpha}_{p})$ is the conditional p.d.f. of $Z$ given ${\it\alpha}_{p}$ and $f_{2}({\it\alpha}_{p})$ in the p.d.f. of ${\it\alpha}_{p}$ (i.e. nearly lognormal). The conditional expected value of $Z$ given ${\it\alpha}_{p}$ , denoted by $\langle Z|{\it\alpha}_{p}\rangle$ , is defined by

and is related to the PA statistic by

For CIT, some prime candidates for conditioning are (i) the vertical fluid velocity component $u_{f,1}$ , which appears in the drag production term, (ii) the vertical velocity difference $u_{p,1}-u_{f,1}$ , which appears in the mean momentum drag exchange term, and (iii) the components of the particle-phase fluctuation energy $k_{p}$ and ${\it\Theta}$ . (Recall that $\langle u_{f,1}\rangle _{f}=0$ and hence $\langle u_{f,1}\mid {\it\alpha}_{p}\rangle$ is equivalent to the conditional expected value of the drift velocity introduced in appendix A.)

In figure 17, $\langle u_{f,1}\mid {\it\alpha}_{p}\rangle /\mathscr{V}$ is plotted versus ${\it\alpha}_{p}/\langle {\it\alpha}_{p}\rangle$ for the three $\mathit{Re}_{p}$ . As was observed in figure 10, for small ${\it\alpha}_{p}$ the conditional drift velocity is positive and for large ${\it\alpha}_{p}$ it is negative, corresponding to stronger fluid entrainment by larger clusters. However, we can also note that the dependence on ${\it\alpha}_{p}$ decreases with increasing $\mathit{Re}_{p}$ , so that, presumably, for large enough $\mathit{Re}_{p}$ the drift velocity would become independent of ${\it\alpha}_{p}$ . As noted in table 2, the fraction of the kinetic energy production that goes to turbulent kinetic energy (i.e.  $\text{DP}/({\it\varphi}\langle u_{p,1}\rangle _{p}g)$ ) decreases with increasing $\mathit{Re}_{p}$ . Because the drift velocity determines DP, the fact that the conditional drift velocity becomes flatter for increasing $\mathit{Re}_{p}$ explains why $\text{DP}/({\it\varphi}\langle u_{p,1}\rangle _{p}g)$ decreases. Physically, the decorrelation of drift velocity with ${\it\alpha}_{p}$ is probably due to the less effective transfer of momentum from the particles to the fluid as $\mathit{Re}_{p}$ increases. Because our simulations are done with a constant ${\it\tau}_{p}$ (i.e. there is no explicit dependence on $\mathit{Re}_{p}$ ), one possible explanation would be that turbulent kinetic energy is dissipated locally in fully resolved wakes and boundary layers around the particles, and therefore it is not effective in producing the large-scale fluid turbulence seen in CIT. Experimental measurements by Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002) of heavy particles in decaying homogeneous isotropic turbulence suggested that clusters act as large pseudo-particles that add to the settling velocity of the particles inside them. Indeed, if clusters are approximated as pseudo-particles with length scale proportional to ${\it\alpha}^{1/3}$ , then the effective cluster Reynolds number will scale like $\mathit{Re}_{c}\approx \mathit{Re}_{p}(1+{\it\beta}{\it\alpha}_{p}^{1/3})$ , where ${\it\beta}$ is a scaling parameter. With this scaling, large clusters will have a significantly larger $\mathit{Re}_{c}$ , implying that they will generate significant fluid-phase vorticity in their wakes and boundary layers. Note that this interpretation is consistent with figure 17, where we observe that the curves flatten out faster for large ${\it\alpha}_{p}$ . In any case, based on the results in figure 17, it is clear that, in order to accurately model the drag production term DP, it will be necessary to account for loss of turbulent kinetic energy production with increasing cluster size.

Turning now to the normalized conditional slip velocity shown in figure 18, we observe that the three curves exhibit good collapse. Recall from § 4.4 that the unconditional slip velocity $|\langle u_{p,1}\rangle _{p}-\langle u_{f,1}\rangle _{p}|$ is very close to $\mathscr{V}$ . Thus, the result in figure 18 implies that

where the function $h>0$ with the property $\langle h\rangle _{p}=1$ is the collapsed curve in figure 18. For example, a simple function that produces the correct asymptotic behaviour for large ${\it\alpha}_{p}$ is $h(x)=0.35+0.65/x$ . Note that, because $\mathscr{V}={\it\tau}_{p}g$ , we can interpret $h(x)$ to be a cluster-size-dependent correction to the drag time scale, i.e.  ${\it\tau}_{c}={\it\tau}_{p}h({\it\alpha}_{p}/\langle {\it\alpha}_{p}\rangle )$ is the effective cluster drag time scale in CIT.

Figure 18.

Slip velocity conditioned on ${\it\alpha}_{p}$ from single realizations of fully developed CIT with $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –), $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ) and $h(x)=0.35+0.65/x$ (○).

In figure 19 the normalized conditional granular energy components are plotted versus ${\it\alpha}_{p}$ . First, note that the uncorrelated granular temperature represents a relatively small fraction of the total conditional granular energy. For the correlated component, it is clear from the curves for conditional $k_{p}$ that no collapse is observed for CIT. In fact, the behaviour for small ${\it\alpha}_{p}$ is quite different over the range of $\mathit{Re}_{p}$ investigated. From these observations, we can conclude that conditioning the particle turbulent kinetic energy on ${\it\alpha}_{p}$ is unlikely to simplify the modelling effort. In comparison, the conditional granular temperature curves in figure 19 show somewhat better collapse. Because the granular temperature is needed to estimate the particle–particle collision rate, it may be worthwhile to attempt to model the conditional granular temperature. However, in CIT we observed that energy dissipation due to inelastic collisions plays a minor role in the overall kinetic energy balance. Thus, given the relatively weak dependence of $\langle {\it\Theta}|{\it\alpha}_{p}\rangle$ on ${\it\alpha}_{p}$ in figure 19, using the unconditional mean $\langle {\it\Theta}\rangle _{p}$ may suffice for CIT.

Figure 19.

Components of the fluctuating granular energy conditioned on ${\it\alpha}_{p}$ from single realizations of fully developed CIT. PA particle-phase kinetic energy $k_{p}$ from $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\boldsymbol{\cdots }\!\!$ ). Granular temperature ${\it\Theta}$ from $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ).

Finally, it is noteworthy that the instantaneous granular temperature dynamics shown in figure 2 do not have a significant effect on $\langle {\it\Theta}|{\it\alpha}_{p}\rangle$ . On closer inspection, the reason for this behaviour is that the value of ${\it\alpha}_{p}$ is not strongly correlated with regions of strong compression (i.e.  $-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{p}\gg 0$ ). For example, the zoomed-in snapshot in figure 2 is a region where ${\it\alpha}_{p}$ changes rapidly and compression is strong. However, on the back side of the same cluster the same range of ${\it\alpha}_{p}$ is observed, but with expansion instead of compression. Overall, the two regions average out for the same value of ${\it\alpha}_{p}$ , leaving no net effect on $\langle {\it\Theta}|{\it\alpha}_{p}\rangle$ .

4.6. Two-point statistics

As already shown in § 3.4.2, two-point statistics contain detailed information of the spatial structure inherent to the turbulent flow field. Along with the two-point particle correlations defined in § 3.4.2, the two-point PA fluid velocity correlation can be expressed as

From this it is possible to define the vertical and spanwise fluid-phase integral length scales as

and

where the subscript $\Vert$ denotes a pair separation ( $r_{1}$ ) measured in the vertical direction and $\bot$ denotes a pair separation ( $r_{2}$ ) measured in the horizontal direction. Note that, owing to the statistical homogeneity of the flows considered in this work, the integral length scales are not a function of space or time. Similarly, particle integral length scales can be defined as

and

Here we have defined the particle-phase integral length scales using the Lagrangian two-point correlations, as they do not depend on the accuracy of the filter and represent the exact spatial distribution of the particles. Another integral length scale of interest can be computed from the RDF in order to define a characteristic cluster size, given by

The integral length scales computed from the two-point velocity correlations can be used to define the following integral time scales for the fluid phase,

and for the particle phase,

The integral length and time scales for the fluid and particle phases are reported in tables 12 and 13, respectively.

Table 12.

Integral length and time scales associated with the fluid phase.

Table 13.

Integral length and time scales associated with the particle phase.

Over the range of $\mathit{Re}_{p}$ considered, we observe that, in general, the integral length scales normalized by $\mathscr{L}$ decrease with increasing $\mathit{Re}_{p}$ . The same trend is observed for the integral time scales. Results for larger $\mathit{Re}_{p}$ would be needed to determine whether these length and time scales become independent of $\mathit{Re}_{f}$ in the high-Reynolds-number limit. In figures 20–22 the radial distribution functions and two-point velocity correlations are plotted using the corresponding integral length scale to normalize the separation distance. For the radial distribution functions in figure 20, an excellent collapse of the curves is found for all $\mathit{Re}_{p}$ . For the two-point velocity correlations in figures 21 and 22, the collapse is not as good, but still reasonable. Moreover, note that the velocity correlations exhibit good decay, confirming that the periodic domain used in our simulations is large enough to avoid finite-domain-size effects. Based on preliminary simulations on smaller domains, we find that, when the periodic domain is too small, particles tend to cluster in vertical columns, which has a strong negative effect on the accuracy of the one-point turbulence statistics. For example, the magnitude of the fluid velocity $\langle u_{f,1}\rangle _{p}$ is significantly underpredicted on smaller domains.

Figure 20.

Radial distribution functions from single realizations of the fully developed CIT computed in the (a) streamwise and (b) spanwise directions. $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ).

Figure 21.

Two-point particle velocity correlations from single realizations of the fully developed CIT simulations. $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ).

Figure 22.

Two-point fluid velocity correlations from single realizations of the fully developed CIT simulations. $\mathit{Re}_{p}=0.25$ (——), $\mathit{Re}_{p}=0.5$ (– – –) and $\mathit{Re}_{p}=1$ ( $\cdots \cdots$ ).

Finally, note that the integral length scales from the radial distribution functions are smaller than those for the fluid velocity correlations. Given that the fluid-phase turbulence is generated by momentum coupling with the clusters in the particle phase, the larger integral length scales in the fluid phase provide a mechanism for breaking up particle clusters on smaller scales, thereby sustaining the nearly lognormal p.d.f. shown in figure 13. On the contrary, if the fluid integral length scales were smaller than the average cluster size (e.g. due to pseudo-turbulence generated by individual particles), then clusters would probably continue to grow to the length scale of the periodic domain (which is observed when the simulation domain is too small). In addition, it can be seen that the length scale $\mathscr{L}$ yields a reasonable prediction for the characteristic cluster size ( $L_{\Vert }^{{\it\alpha}}$ and $L_{\bot }^{{\it\alpha}}$ ), but, owing to the formation of wakes, this length scale severely underpredicts the velocity length scales of each phase. Nonetheless, the exact mechanism that determines the integral length scales of the clusters and their one-point p.d.f. has yet to be elucidated, and thus the scaling behaviour of CIT as a function of $\mathit{Re}_{f}$ remains an open question.