2.1 Kinetic-theory-based two-fluid model

Table 1. Transport equations and constitutive relations of the GTSH KT-TFM. The solution variables are the solids concentration, $\unicode[STIX]{x1D719}$ , the solids and fluid velocity vectors, $\boldsymbol{U}_{j}=(u_{j},v_{j},w_{j})^{\text{T}}$ with $j=s,f$ , the fluid pressure, $p_{f}$ , and the granular temperature, $T$ . Additionally, $\unicode[STIX]{x1D70C}_{j}$ is the phasic density, $\boldsymbol{g}$ is the gravity vector, $e$ is the restitution coefficient, $d_{p}$ is the particle diameter and $m$ is the particle mass.

This study utilizes a two-fluid model (TFM) to solve gas–solid flow in a state of fluidization or sedimentation in a fully periodic domain. The solids-phase transport coefficients and constitutive relations are closed with a KT approach. We use the KT-TFM of Garzó et al. (Reference Garzó, Tenneti, Subramaniam and Hrenya2012) (hereafter referred to as GTSH). Unlike previous KT models for particulate flows, the instantaneous gas–solid interaction force appearing in the starting Enskog equation is decomposed into three components: (i) a mean drag force (proportional to $\unicode[STIX]{x1D6FD}$ in table 1) that models the interfacial momentum transfer due to a difference in mean relative velocities, (ii) a thermal drag force (proportional to $\unicode[STIX]{x1D6FE}$ in table 4) that models the dissipation of granular temperature (a measure of the fluctuating kinetic energy of the solids phase) due to the fluid phase, and (iii) a neighbour effect (proportional to $\unicode[STIX]{x1D709}$ in table 5), which is a source of granular temperature due to the fluid and the presence of nearby particles, i.e. it is necessarily both a multiphase and multi-particle contribution. For the sake of brevity, only the governing and constitutive equations of the KT-TFM are given in table 1. The additional relations needed to fully specify the GTSH model are provided in tables 3–5 in appendix A along with a discussion of the slight differences between closures used here and those originally proposed by Garzó et al. (Reference Garzó, Tenneti, Subramaniam and Hrenya2012). In §§ 2.1.1, 2.1.2 and 2.1.3 we provide a review of the relevant assumptions and approximations that impact the limitations of the GTSH theory.

2.2 Numerical solution of the KT-TFM

The GTSH KT-TFM is solved numerically using the open-source MFiX CFD code developed at the National Energy Technology Laboratory (https://mfix.netl.doe.gov/). MFiX solves the discretized transport equations with the finite volume method. The spatial discretization employs a staggered grid and an upwinded variable extrapolation scheme with a flux limiter. All simulations in this work were performed with the MUSCL flux limiter (van Leer Reference van Leer1979; Waterson & Deconinck Reference Waterson and Deconinck2007) on a uniform grid of $N_{x}=N_{z}=12$ and $N_{y}=50$ . Although the grid is uniform, it is not exactly cubic and the exact grid spacing $\unicode[STIX]{x1D6E5}_{i}=L_{i}/N_{i}$ for $i=x,y,z$ , changes slightly depending on the Archimedes number, see § 2.4 below. For both Archimedes numbers, the non-dimensionalized grid size, $\unicode[STIX]{x1D6E5}^{\ast }=(\unicode[STIX]{x1D6E5}_{x}\times \unicode[STIX]{x1D6E5}_{y}\times \unicode[STIX]{x1D6E5}_{z})^{1/3}/d_{p}$ , is approximately 0.7 which was determined to be necessary after carrying out a systematic grid dependence study of a test case. Such a fine grid is approximately an order of magnitude smaller than previously found to be necessary for grid-insensitive solutions (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Wang et al. Reference Wang, van der Hoef and Kuipers2009; Fullmer & Hrenya Reference Fullmer and Hrenya2016). It is believed that such a high level of numerical resolution is needed here due to the small system size; specifically, large spatial gradients at cluster interfaces account for a larger portion of the total domain in this study (Li et al. Reference Li, Wang, Rogers, Zhou and Ge2017). Time advancement is semi-implicit with a SIMPLE-type algorithm (Patankar Reference Patankar1980) for pressure–velocity coupling. The temporal discretization is first-order Euler with adaptive time stepping. A maximum time step of is specified, however in order to meet the tolerance limits, the code-limited time step is typically one to two orders of magnitude smaller than the maximum time step.

2.3 DNS method

The SUSP3D code developed by Ladd and co-workers (Ladd & Verberg Reference Ladd and Verberg2001; Nguyen & Ladd Reference Nguyen and Ladd2002) is used in this study to simulate the dynamics of fully periodic sedimentation/fluidization. This code uses a D3Q19 lattice Boltzmann velocity model and a two-relaxation-time collision operator to solve fluid flow on a uniform, space filling, cubic lattice. Incompressible flow is recovered with second-order accuracy with respect to the lattice Mach number, $Ma=u/c_{s}$ , where $u$ is the fluid velocity and $c_{s}=1/\sqrt{3}$ is the lattice speed of sound. When a fluid population crosses a fluid–solid interface, the linked-bounce-back rule is applied to return the fluid population back to the fluid domain, recovering the no-slip boundary condition. The momentum exchanged between the fluid and a particle through bounce-back is accumulated over the surface of the particle to give the fluid–particle force. This force is then used to update the particles velocity according to Newton’s law of motion.

For close pairs of particles, analytical expressions of lubrication interactions are imposed onto lattice Boltzmann solutions to ensure that lubrication interactions between particles are correctly captured (Nguyen & Ladd Reference Nguyen and Ladd2002). As analytical expressions contain singularities when the gap between particles approaches zero, a lubrication cutoff, $\unicode[STIX]{x1D716}$ is specified. This cutoff sets maximum values on the lubrication force and torque, and reduces the stiffness of the velocity update. As in the assumptions used to derive the GTSH KT-TFM, collisions between particles are treated as hard sphere collisions with a normal restitution coefficient. Tangential restitution due to friction was not considered in the simulations.

Table 2. Parameters used in the DNS.

Table 2 lists the properties of the lattice fluid, diameters of particles and lubrication cutoffs that were employed in the simulations. Note that the diameter $d_{p}$ is the effective hydrodynamic diameter, which is different from the geometric diameter, $d_{p}^{0}$ , used to distinguish fluid and solid nodes on the computational lattice. Generally $d_{p}>d_{p}^{0}$ , because the collection of solid nodes mapped by $d_{p}^{0}$ has a rough interface that comes from the staircase representation of a sphere. The difference between $d_{p}$ and $d_{p}^{0}$ , the hydrodynamic correction $\unicode[STIX]{x1D6E5}_{H}$ , is typically established through drag calibration (Ladd & Verberg Reference Ladd and Verberg2001), and is always less than one lattice spacing. As one increases $d_{p}$ and $d_{p}^{0}$ , the need for a hydrodynamic correction decreases, as the sphere is now represented by more lattice units (i.e. higher resolution). The choice of $d_{p}$ is generally based on the balance between computational cost and accuracy. In this study, we tested a higher resolution $d_{p}=10$ for selected cases, and the results were identical to those obtained using the lower $d_{p}$ listed in table 2. It was also tested and confirmed that the DNS results are essentially not affected by the selection of $\unicode[STIX]{x1D716}$ .

2.4 System description

The system considered is a uniform suspension in a fully periodic domain undergoing infinite sedimentation or fluidization. A constant body force, i.e. gravity, acts only in the negative vertical $y$ -direction. In the KT-TFM simulations, the fluid pressure is assumed to be decomposed into a fluctuating (local) component superimposed on a mean linear gradient. The mean fluid pressure gradient opposes the body force and supports the weight of the suspension,

(2.11) $$\begin{eqnarray}-\!\frac{\unicode[STIX]{x2202}P_{f}}{\unicode[STIX]{x2202}y}=(\langle \unicode[STIX]{x1D719}\rangle \unicode[STIX]{x1D70C}_{s}+(1-\langle \unicode[STIX]{x1D719}\rangle )\unicode[STIX]{x1D70C}_{f})|\boldsymbol{g}|\end{eqnarray}$$

so that no net force acts on the mixture in the $y$ -direction. The angle brackets in (2.11) indicate an averaging over the entire domain. The pressure gradient is imposed on the system numerically by specifying a pressure-drop boundary condition in the $y$ -direction. In DNS, a body force equivalent to that specified by (2.11) is applied to the fluid to ensure that the net force on the mixture in the $y$ -direction is zero. No forces are applied in the transverse $x$ - and $z$ -directions and standard periodic boundary conditions are used. Since the system is fully periodic, the global (system) concentration, $\langle \unicode[STIX]{x1D719}\rangle$ , does not change in time.

With the specification of (2.11) one can calculate the base state, i.e. the homogeneous (stable) steady state, which will be an important metric for comparison with the numerical solutions of the complete (transient and three-dimensional) GTSH KT-TFM. The mixture momentum equation, (2.3) $+$ (2.4), reduces to (2.11). The difference momentum equation, (2.3)–(2.4), can be rearranged to give the mean velocity difference in the vertical direction:

(2.12) $$\begin{eqnarray}|v_{s}^{({\mathcal{H}})}-v_{f}^{({\mathcal{H}})}|=\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}|\boldsymbol{g}|d_{p}^{2}}{18\unicode[STIX]{x1D707}_{f}F^{\ast }}=\frac{u_{\infty }}{F^{\ast }},\end{eqnarray}$$

where $u_{\infty }$ is the terminal Stokes velocity of a single particle in an infinite medium and $F^{\ast }$ is the mean drag law (A 1) and the $({\mathcal{H}})$ superscript denotes the homogeneous state. In §§ 3 and 4, it will be important to distinguish between values calculated analytically from the KT-TFM assuming homogeneity and those calculated from the results of CFD simulations of the KT-TFM which may be far from homogeneous. The transverse velocities in the homogeneous state are all zero: $u_{s}^{({\mathcal{H}})}=u_{f}^{({\mathcal{H}})}=w_{s}^{({\mathcal{H}})}=w_{f}^{({\mathcal{H}})}=0$ . Note that (2.12) only prescribes the magnitude of the velocity difference, not the phasic velocity values. The mixture continuity equation, (2.1) $+$ (2.2), states only that the total volumetric flux should be a constant,

(2.13) $$\begin{eqnarray}j_{y}^{({\mathcal{H}})}=[\unicode[STIX]{x1D719}v_{s}+(1-\unicode[STIX]{x1D719})v_{f}]^{({\mathcal{H}})}=\text{const.}\end{eqnarray}$$

In theory, the absolute values of the phasic velocities are arbitrary: hence the lack of distinction between sedimentation or fluidization in a fully periodic domain. However, to satisfy the $Ma$ criterion in the DNS, the velocities should be minimized and thus the choice $j_{y}^{({\mathcal{H}})}=0$ was selected and also used in the KT-TFM simulations for consistency. With $j_{y}^{({\mathcal{H}})}=0$ , equation (2.12) yields $v_{s}^{({\mathcal{H}})}=-(1-\langle \unicode[STIX]{x1D719}\rangle )u_{\infty }/F^{\ast }$ and $v_{f}^{({\mathcal{H}})}=\langle \unicode[STIX]{x1D719}\rangle u_{\infty }/F^{\ast }$ (note that $\unicode[STIX]{x1D719}^{({\mathcal{H}})}=\langle \unicode[STIX]{x1D719}\rangle$ ). In the homogeneous state, the granular temperature equation (2.5) reduces to the algebraic relation $S_{0}=0$ or,

(2.14) $$\begin{eqnarray}\left[\unicode[STIX]{x1D709}=\frac{2\unicode[STIX]{x1D6FE}}{m}T+\unicode[STIX]{x1D701}_{0}T\right]^{({\mathcal{H}})},\end{eqnarray}$$

which describes how granular temperature generated by the neighbour effect is balanced the by the dissipation due to thermal drag and inelasticity (if $e<1$ ).

For both numerical methods, the fluid and particles are initially at rest and accelerated by body forces to achieve the quasi-steady state of sedimentation. Though by different means, both methods enforce $\langle j_{y}\rangle =0$ dynamically throughout calculations to avoid drift caused by round-off error. To provide an initial perturbation in the KT-TFM simulations consistent with DNS randomly distributed particles are generated at initialization and filtered into solids concentration on the CFD grid with a simple Gaussian filter. Furthermore, in the KT-TFM simulations a negligibly small initial granular temperature is prescribed to avoid singularities in the constitutive equations. In this manner, the KT-TFM and DNS methods were initialized as similarly as possible.

Seven non-dimensional parameters characterize the system: the Archimedes number, $Ar=\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}|\boldsymbol{g}|d_{p}^{3}/\unicode[STIX]{x1D707}_{f}^{2}$ , characterizing the ratio of gravity-to-viscous forces, the density ratio, $\unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{s}/\unicode[STIX]{x1D70C}_{f}$ , the mean concentration, $\langle \unicode[STIX]{x1D719}\rangle$ , the coefficient of restitution, $e$ , and the system size in each direction, $L_{i}^{\ast }=L_{i}/d_{p}$ where $i=x,y,z$ . Two Archimedes numbers considered in this study, $Ar=71$ and 1432, lead to approximately one order of magnitude of variation in mean-flow Reynolds numbers, equation (A 26). (Note that (2.12) can be easily non-dimensionalized into the form $Re_{m}^{({\mathcal{H}})}=(1-\langle \unicode[STIX]{x1D719}\rangle )Ar/18F^{\ast }$ ). The density ratio serves as the control parameter and was varied from 10 to 1000. Four mean concentrations are studied, $\langle \unicode[STIX]{x1D719}\rangle =0.10$ , 0.15, 0.25 and 0.40 and two values of the restitution coefficient are used, $e=0.9$ and 1.0. Here, we primarily focus on the results obtained in the elastic case ( $e=1.0$ ) as inelasticity does not significantly affect the results as discussed in § 4.6. The system size varies slightly between the two Archimedes number cases. For the case $Ar=71$ : $L_{x}^{\ast }=L_{z}^{\ast }=8.56$ and $L_{y}^{\ast }=34.2$ . For the $Ar=1432$ case: $L_{x}^{\ast }=L_{z}^{\ast }=8.656$ and $L_{y}^{\ast }=34.624$ . In both instances, the domain size is very small – a consequence of the large computational cost associated with DNS. Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2015) recommend that the domain size in the flow direction be at least $L_{y}^{\ast }\approx 0.2\unicode[STIX]{x1D70C}^{\ast }Ar$ (assuming $\unicode[STIX]{x1D70C}^{\ast }\gg 1$ ). In this work, $0.2\unicode[STIX]{x1D70C}^{\ast }Ar$ ranges from 142 to 286 400 – orders of magnitude larger than $L_{y}^{\ast }$ considered here – so that the periodic boundaries influence the development of the instabilities. Therefore, when clustering occurs it forms a significantly different structure than observed in larger, truly unbounded domains, e.g. Capecelatro et al. (Reference Capecelatro, Desjardins and Fox2015), Capecelatro, Desjardins & Fox (Reference Capecelatro, Desjardins and Fox2016). However, the variety of structures observed as a result of the small system provide an equally, if not more, stringent test of the KT-TFM than traditional clustering alone, as shown in § 3.

2.5 Averaging operations

A few averaging operations that will be used in the following sections are formally defined here. Consider a general, arbitrary variable $f_{j}$ associated with phase- $j$ which is a function of space and time. The volume average,

(2.15) $$\begin{eqnarray}\langle f_{j}\rangle (t^{\ast })=\frac{1}{L_{x}^{\ast }L_{y}^{\ast }L_{z}^{\ast }}\int _{0}^{L_{z}^{\ast }}\int _{0}^{L_{y}^{\ast }}\int _{0}^{L_{x}^{\ast }}f_{j}(\boldsymbol{x}^{\ast },t^{\ast })\,\text{d}x^{\ast }\,\text{d}y^{\ast }\,\text{d}z^{\ast },\end{eqnarray}$$

defines the global or system-averaged quantity. When comparing with the results of DNS, continuum flow variables, e.g. solids velocity, must be Favre or concentration weighted,

(2.16a,b ) $$\begin{eqnarray}\langle \langle f_{s}\rangle \rangle (t^{\ast })=\frac{\langle \unicode[STIX]{x1D719}f_{s}\rangle (t^{\ast })}{\langle \unicode[STIX]{x1D719}\rangle (t^{\ast })}\quad \text{and}\quad \langle \langle f_{f}\rangle \rangle (t^{\ast })=\frac{\langle (1-\unicode[STIX]{x1D719})f_{f}\rangle (t^{\ast })}{1-\langle \unicode[STIX]{x1D719}\rangle (t^{\ast })},\end{eqnarray}$$

denoted here by double brackets. Finally, the system is simulated in a statistical steady state and it is frequently desired to report time-averaged values,

(2.17) $$\begin{eqnarray}\overline{f_{j}}(\boldsymbol{x}^{\ast })=\frac{1}{t_{2}^{\ast }-t_{1}^{\ast }}\int _{t_{1}^{\ast }}^{t_{2}^{\ast }}f_{j}(\boldsymbol{x}^{\ast },t^{\ast })\,\text{d}t^{\ast }\end{eqnarray}$$

denoted with an overbar. The non-dimensional time is defined here as $t^{\ast }=t\unicode[STIX]{x1D708}_{f}/d_{p}^{2}$ . Since DNS is extremely computationally expensive, the time averaging windows are minimized by beginning time averaging as soon as the quasi-steady state is reached and terminating each simulation when a sufficient amount of statistics have been collected, typically controlled by the higher-order statistics. Five replicates (otherwise identical simulations with a different randomization of initial particle distribution) were run for each condition. The DNS results provided in § 4 represent the average of the five replicates with error bars signifying the standard deviation between the replicates. The KT-TFM simulations are carried out to $t^{\ast }=2000$ for $Ar=71$ and to $t^{\ast }=500$ for $Ar=1432$ . In both cases, the data are split into ten equal length segments. The first segment is neglected to avoid any transient data associated with the transition from the initially static but perturbed state to the statistically steady state. The remaining nine segments are time averaged and used to calculate the mean and standard deviation, the latter of which is represented by error bars throughout this work except in figure 11 where specifically noted otherwise.

4.1 Analysis of the Stokes number and the limit of GTSH theory

In the following §§ 4.2–4.4, the KT-TFM and DNS results are directly compared quantitatively for two Archimedes numbers, several mean concentrations and a wide range of solid-to-fluid density ratios. In general, the best agreement is observed at the highest density ratios and a critical analysis of the assumptions overviewed in §§ 2.1.1–2.1.3 was carried out to determine why some regions of the phase space show better comparisons than others. A key benefit of comparing to ideal DNS numerical data is that the particle-phase properties and most of the fluid-phase properties are identical between KT-TFM and DNS. Moreover, complex physics commonly encountered in experiments (e.g. cohesion, humidity, etc.) are entirely neglected in both KT-TFM and the validation DNS data, providing as much of an ideal comparison as possible. The only physical condition that differs between KT-TFM and DNS is related to particle-induced turbulence; the fine grid used in the KT-TFM simulations actively capture cluster-induced turbulence, but small-scale fluctuations generated by the flow around individual particles are not specifically modelled. The impact of this discrepancy is gauged with Sato’s simple mixing length model for particle-induced turbulence, $\unicode[STIX]{x1D708}_{PIT}=C_{PIT}\unicode[STIX]{x1D719}d_{p}|\unicode[STIX]{x0394}U|$ , (Sato, Sadatomi & Sekoguchi Reference Sato, Sadatomi and Sekoguchi1981). Using the suggested value of $C_{PIT}=0.6$ , an empirical fitting coefficient originally developed for bubbly flows, and the homogeneous (stable) solution (2.12), the ratio $\unicode[STIX]{x1D708}_{PIT}/\unicode[STIX]{x1D708}_{f}$ ranges from 0.09 to 0.11 for $Ar=71$ and from 1.24 to 1.85 for $Ar=1432$ . Therefore, assumptions in the KT derivation, § 2.1.3, are the chief suspects for observed discrepancies. We focus here on the assumed particle distribution function (high thermal Stokes number assumption) and reserve a discussion of the small Knudsen number assumption for § 4.5.

As discussed in § 2.1.3, the GTSH theory allows for a non-Maxwellian velocity distribution function, $f$ . However, since $f$ is expressed as an expansion of Sonine polynomials (to first order), $f$ is expected to be near Maxwellian, its deviation from which is measured via (A 22) and relaxed predominantly by particle–particle collisions. In the context of a fluid–particle suspension, this assumption amounts to a requirement that the particle–particle (collisional) time scale should be much shorter than the fluid–particle (viscous relaxation) time scale (Koch Reference Koch1990) – i.e. a high thermal Stokes number (as defined below). We take the mean collision time to be the ratio of the mean free path due to the mean (fluctuating) speed. In a dense gas, the mean free path is given by $\ell =d_{p}/6\sqrt{2}\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}$ (Garzó Reference Garzó2005). Ignoring the Maxwellian deviation given by (A 22), the mean speed is given by $2c/\sqrt{\unicode[STIX]{x03C0}}$ where $c=\sqrt{2T}$ is the thermal (most probable) speed. Combining, the mean collision time is approximated by

(4.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{coll}=\frac{d_{p}}{24\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}}\sqrt{\frac{\unicode[STIX]{x03C0}}{T}}.\end{eqnarray}$$

The viscous relaxation time is defined as the time constant in the exponential velocity decay of a test particle in a Lagrangian frame of reference (Wallis Reference Wallis1969). In multi-particle suspension at finite $Re_{m}$ , the viscous relaxation time is

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{visc}=\frac{m}{3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{f}d_{p}F^{\ast }},\end{eqnarray}$$

where $m$ is the particle mass and $F^{\ast }$ is the Stokes deviation of the mean drag law (A 1). Combining (4.1) and (4.2) the ratio for the viscous-to-collisional time scales becomes

(4.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70F}_{visc}}{\unicode[STIX]{x1D70F}_{coll}}=\frac{4}{3\sqrt{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}}{F^{\ast }}\unicode[STIX]{x1D70C}^{\ast }Re_{T}\gg 1,\end{eqnarray}$$

where the product of the density ratio and the thermal Reynolds number has been combined into the thermal Stokes number, $St_{T}=\unicode[STIX]{x1D70C}^{\ast }Re_{T}/9$ (Tenneti & Subramaniam Reference Tenneti and Subramaniam2014), which is more commonly used to characterize the particle response time. The validity condition of (4.3) is vague due to the ‘much greater than’ inequality. For demonstration purposes, we assume that half an order of magnitude ( $10^{1/2}$ ) is sufficient, as its square represents one order of magnitude change. Also note that (4.3) is a local instant quantity which should be averaged in some fashion. Accordingly, it is expected that the assumptions of the GTSH theory (nearly Maxwellian, or equivalently, high $St_{T}$ ) are approximately upheld if the space- and time-averaged viscous-to-collision time scale ratio, ${\mathcal{T}}$ , is at least $10^{1/2}$ , or

(4.4) $$\begin{eqnarray}{\mathcal{T}}=\overline{\left\langle \frac{12}{\sqrt{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}}{F^{\ast }}St_{T}\right\rangle }\geqslant 10^{1/2}.\end{eqnarray}$$

Figure 5 summarizes the calculation of ${\mathcal{T}}$ for the elastic conditions studied. The intersections of the ${\mathcal{T}}$ -curves with $10^{1/2}$ , i.e. the equality of (4.4), demarks a critical density ratio which is indicated by the vertical red lines in figures 6–9. As discussed below (§§ 4.2 and 4.3), the results indicate that this high- $St_{T}$ assumption (4.4) of the continuum theory is needed for accurate predictions.

Figure 5. Maximum normalized concentration difference in the KT-TFM simulations as a function of density ratio for $e=1$ and (a) $Ar=71$ and (b) $Ar=1432$ at mean concentrations of $\langle \unicode[STIX]{x1D719}\rangle =0.1$ (black), 0.15 (red), 0.25 (blue) and 0.40 (green). The thin dashed lines of similar colour show ${\mathcal{T}}$ in the homogeneous state.

4.2 Mean solids velocity

In this section we focus on mean solids velocity and make direct quantitative comparisons between the KT-TFM predictions and the DNS data. For the comparison, the continuum variable that is the solids velocity is Favre averaged. In the limit $\unicode[STIX]{x1D6E5}^{\ast }\rightarrow 0$ , $\unicode[STIX]{x1D719}$ becomes a phase indicator function and the Favre average becomes equivalent to a simple (discrete) particle average, consistent with the DNS results. Through a combination of volume, Favre and time averaging, equations (2.15)–(2.17), the space- and time-dependent data for an entire simulation can be reduced to a set of mean, scalar quantities for a given set of conditions ( $Ar$ , $\unicode[STIX]{x1D70C}^{\ast }$ , $\langle \unicode[STIX]{x1D719}\rangle$ , $e$ ). We present the mean sedimenting Reynolds number,

(4.5) $$\begin{eqnarray}Re_{S}=\frac{\unicode[STIX]{x1D70C}_{f}d_{p}}{\unicode[STIX]{x1D707}_{f}}|\overline{\langle \langle v_{s}\rangle \rangle }|,\end{eqnarray}$$

which measures the mean streamwise (vertical) velocity of the solids. Non-dimensionalization is required to compare the results of KT-TFM simulations performed in physical units with DNS simulations performed in lattice units. Here, $Re_{S}$ is preferred to $Re_{m}$ because the averaging is simpler; however, the two measurements are roughly similar. In fact, in the homogeneous state it can be shown that $Re_{S}^{({\mathcal{H}})}=Re_{m}^{({\mathcal{H}})}$ due to the specification $j_{y}^{({\mathcal{H}})}=0$ .

Figure 6 compares KT-TFM predictions of $Re_{S}$ to the DNS data as a function of $\unicode[STIX]{x1D70C}^{\ast }$ for $Ar=71$ , $e=1.0$ and the four mean concentrations: $\langle \unicode[STIX]{x1D719}\rangle =0.10$ , 0.15, 0.25 and 0.40. Each DNS data point represents an average of five replicate simulations with different initial conditions, although each simulation is run for a shorter time than the continuum simulations. The standard deviation of the five DNS replicates are shown as error bars in figure 6. In most cases, the error bars are too tight to extend above/below the data points representing the mean. The horizontal dashed lines show $Re_{S}^{({\mathcal{H}})}$ calculated from (2.12), which is independent of $\unicode[STIX]{x1D70C}^{\ast }$ for the drag law, equations (A 1)–(A 3), of Beetstra, van der Hoef & Kuipers (Reference Beetstra, van der Hoef and Kuipers2007) used here. The vertical red lines indicate the critical density ratio needed to uphold the high- $St_{T}$ (nearly Maxwellian) assumption of (4.4). The KT-TFM theory is expected to be valid to the right of the red line, i.e. at higher $\unicode[STIX]{x1D70C}^{\ast }$ . Generally, good agreement is observed in these regions, which appears to be a conservative estimate in some instances, especially at lower concentrations. Although the exact shape of the four KT-TFM $Re_{S}$ -curves in figure 6 are different, a similar trend is observed: $Re_{S}$ is largest at an intermediate density ratio ( $\unicode[STIX]{x1D70C}^{\ast }\approx 100$ ) and appears to decay at low and high density ratios. At lower density ratios, $Re_{S}$ decreases because the KT-TFM solution approaches a homogeneous state ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{max}\rightarrow 0$ in figure 1) which is most evident at higher concentrations as $Re_{S}$ coincides with $Re_{S}^{({\mathcal{H}})}$ over a range of $\unicode[STIX]{x1D70C}^{\ast }$ . The discrepancy between the DNS results and $Re_{S}^{({\mathcal{H}})}$ (determined from a DNS-based drag law) at low density ratio may be due to the specification of a fixed-bed-type drag law. Recently, Rubinstein et al. (Reference Rubinstein, Derksen and Sundaresan2016) postulated that a Wen & Yu (Reference Wen and Yu1966) type drag law is more appropriate for such low mean-flow Stokes number flows. Although the incorporation of such a mean-flow-Stokes-number-dependent drag law may also improve the KT-TFM predictions in this region, the most significant difference would be expected to occur outside of the intended range of validity, i.e. to the left of the red lines.

Figure 6. Comparison of KT-TFM predictions (thick lines) to DNS data (symbols) of the sedimenting Reynolds number for $Ar=71$ and (a) $\langle \unicode[STIX]{x1D719}\rangle =0.10$ , (b) $\langle \unicode[STIX]{x1D719}\rangle =0.15$ , (c) $\langle \unicode[STIX]{x1D719}\rangle =0.25$ and (d) $\langle \unicode[STIX]{x1D719}\rangle =0.40$ . Black dashed lines show homogeneous (stable) solutions, $Re_{S}^{({\mathcal{H}})}$ . Right-hand side of vertical red lines demarcates high- $St_{T}$ region of (4.4).

Conversely at large $\unicode[STIX]{x1D70C}^{\ast }$ , the suspension remains inhomogeneous ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{max}>0$ in figure 1) yet $Re_{S}$ decays to the point of $Re_{S}^{({\mathcal{H}})}$ and below, in some instances. A hindered sedimenting velocity, i.e. $Re_{S}<Re_{S}^{({\mathcal{H}})}$ , is a consequence of an increase in mean drag – a counterintuitive result since the clustering instability of gas–solid particle flow is most often associated with a mean drag reduction occasionally referred to as flow bypass, e.g. Agrawal et al. (Reference Agrawal, Loezos, Syamlal and Sundaresan2001), Li & Kwauk (Reference Li and Kwauk2001), Heynderickx et al. (Reference Heynderickx, Das, De Wilde and Marin2004), Beetstra, van der Hoef & Kuipers (Reference Beetstra, van der Hoef and Kuipers2006), Wang et al. (Reference Wang, Lu, Liu, Sheng, Xu and Gidaspow2011), Shah et al. (Reference Shah, Utikar, Tade, Evans and Pareek2013), Zhou et al. (Reference Zhou, Xiong, Wang, Wang, Ren and Ge2014) among others. Note that $Re_{S}<Re_{S}^{({\mathcal{H}})}$ at $\unicode[STIX]{x1D70C}^{\ast }\approx 1000$ is observed in the DNS data as well. The observed decay in $Re_{S}$ at high $\unicode[STIX]{x1D70C}^{\ast }$ is due to the change in flow structure discussed previously, namely the transition from the unsteady chaotic regime to the steady plug regime. For horizontal structure formed in the chaotic regime, the fluid has enough inertia to penetrate the plug and flow bypass is possible, although intermittent (see figure 2). In the steady (nonlinearly stable) plug regime, the plug structure is persistent and the fluid is forced to flow through the width-spanning, dense cluster (see figure 3). As the density ratio increases, the plug becomes more one-dimensional which causes more of the solids concentration to move from the funnel region into the horizontal plug region, further decreasing the sedimenting velocity. While the clustering structures observed here are a result of the small domain size, the accurate prediction of the counter-intuitive results displayed in figure 6 bodes well for the predictive capability of the GTSH KT-TFM, at least where it is expected to be valid (to the right of the vertical red line).

Figure 7. Comparison of KT-TFM predictions (thick lines) to DNS data (symbols) of the sedimenting Reynolds number for $Ar=1432$ and (a) $\langle \unicode[STIX]{x1D719}\rangle =0.10$ , (b) $\langle \unicode[STIX]{x1D719}\rangle =0.15$ , (c) $\langle \unicode[STIX]{x1D719}\rangle =0.25$ and (d) $\langle \unicode[STIX]{x1D719}\rangle =0.40$ . Black dashed lines show homogeneous (stable) solutions, $Re_{S}^{({\mathcal{H}})}$ . Right-hand side of vertical red lines demarcates high- $St_{T}$ region of (4.4).

Figure 7 shows the same comparison as in figure 6 but at the higher $Ar=1432$ condition, which results in a wholesale increase of $Re_{S}$ by approximately one order of magnitude. As in figure 6, the red lines approximate where the high- $St_{T}$ condition (4.4) is valid, but are shifted to significantly lower density ratios compared to the $Ar=71$ case. Again, good agreement except for a slight over-prediction is generally seen where the theory is expected to be valid, and this agreement deteriorates to the left of the red line. The over-prediction of the DNS data by KT-TFM – even in the valid range – occurs at all concentrations, hinting that there may be a discrepancy between the mean drag law and the current set of DNS data. This discrepancy may be due, in part, to the increased uncertainty in the mean drag law at this higher Reynolds number, e.g. see comparisons in Tenneti, Garg & Subramaniam (Reference Tenneti, Garg and Subramaniam2011), Tang et al. (Reference Tang, Peters, Kuipers, Kriebitzsch and van der Hoef2015). Similar to the $Ar=71$ case, figure 7 shows that $Re_{S}$ decays with increasing $\unicode[STIX]{x1D70C}^{\ast }$ and a hindered mean sedimenting velocity, $Re_{S}<Re_{S}^{({\mathcal{H}})}$ , is observed in the KT-TFM predictions, consistent with the DNS data. However, two significant distinctions relative to the lower $Ar$ results (figure 6) can be seen: (i) the lowest density ratios are no longer small enough to observe the near-homogeneous state with the KT-TFM and (ii) the decay of $Re_{S}$ with $\unicode[STIX]{x1D70C}^{\ast }$ is not unbounded; at a given density ratio (in the vicinity of $\unicode[STIX]{x1D70C}^{\ast }\approx 200$ ) the trend changes, rather abruptly in some instances, and $Re_{S}$ begins to increase. The latter is a new feature of higher $Ar$ and one that is captured remarkably well by the KT-TFM. The local minima in the sedimenting velocity curves of figure 7 coincide with the loss of the 2-D plug structure, i.e. transition in flow structure from that of figure 3 to that of figure 4. As the density ratio continues to increase, the peak concentration of the plug is reduced as the magnitude of the 1-D wave relaxes, leading to a reduction in the maximum drag experienced. Although the dilute regions are also less dilute, the strong nonlinearity in the concentration dependence of the drag law (A 1)–(A 3) causes an overall decrease in drag, causing the increase in $Re_{S}$ observed in figure 7. A few exploratory studies were carried out on a coarser grid to investigate the increasing $Re_{S}$ behaviour beyond $\unicode[STIX]{x1D70C}^{\ast }=1000$ . Up to $\unicode[STIX]{x1D70C}^{\ast }=100\,000$ , $Re_{S}$ does not cross the homogeneous solution again but asymptotes to it, i.e. $Re_{S}<Re_{S}^{({\mathcal{H}})}$ as $\unicode[STIX]{x1D70C}^{\ast }\rightarrow \infty$ .

4.3 Fluctuating solids velocity

In this section we move from mean to fluctuating statistics. Before jumping directly to a comparison with the DNS data, the KT-TFM variables will be reviewed, specifically how the continuum variables are defined from the probabilistic KT variables. The three variables used in the governing equations in table 1, concentration, mean velocity and granular temperature, are defined by

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}={\mathcal{V}}_{p}\displaystyle \int f(\boldsymbol{x},\boldsymbol{v}_{p},t)\,\text{d}\boldsymbol{v}_{p}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{U}_{s}=\frac{{\mathcal{V}}_{p}}{\unicode[STIX]{x1D719}}\int \boldsymbol{v}_{p}f(\boldsymbol{x},\boldsymbol{v}_{p},t)\,\text{d}\boldsymbol{v}_{p}, & \displaystyle\end{eqnarray}$$

and

(4.8) $$\begin{eqnarray}T=\frac{1}{3}\frac{{\mathcal{V}}_{p}}{\unicode[STIX]{x1D719}}\int (\boldsymbol{v}_{p}-\boldsymbol{U}_{s})^{2}f(\boldsymbol{x},\boldsymbol{v}_{p},t)\,\text{d}\boldsymbol{v}_{p},\end{eqnarray}$$

respectively (Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012). In (4.6)–(4.8), ${\mathcal{V}}_{p}$ is the particle volume, $\boldsymbol{v}_{p}=(u_{p},v_{p},w_{p})^{\text{T}}$ is the particle velocity and $f$ is the single-particle velocity distribution function. Similarly, the total (solids) kinetic energy can be defined by

(4.9) $$\begin{eqnarray}E_{s}=\frac{1}{2}\frac{{\mathcal{V}}_{p}}{\unicode[STIX]{x1D719}}\int \boldsymbol{v}_{p}^{2}f(\boldsymbol{x},\boldsymbol{v}_{p},t)\,\text{d}\boldsymbol{v}_{p}.\end{eqnarray}$$

The total kinetic energy is infrequently encountered in the literature because it is redundant, i.e. it can be equivalently determined by

(4.10) $$\begin{eqnarray}E_{s}={\textstyle \frac{1}{2}}(\boldsymbol{U}_{s}^{2}+3T).\end{eqnarray}$$

Although (4.10) could be directly compared with discrete particle data, it is more informative to introduce an approximation to decompose the kinetic energy into its directional components in order to study any anisotropy present in the system. By assuming that the (microscopic) fluctuating velocities are isotropic, i.e. $u_{p}u_{s}\approx v_{p}v_{s}\approx w_{p}w_{s}\approx T/3$ , equation (4.8) can be manipulated to give

(4.11a-c ) $$\begin{eqnarray}E_{s,x}={\textstyle \frac{1}{2}}(u_{s}^{2}+T),\quad E_{s,y}={\textstyle \frac{1}{2}}(v_{s}^{2}+T),\quad \text{and}\quad E_{s,z}={\textstyle \frac{1}{2}}(w_{s}^{2}+T).\end{eqnarray}$$

Now, considering again the analogy between a simple discrete particle average and a volume average with infinite resolution. Previously it was mentioned that the appropriate quantity to compare to the particle mean velocity (discrete particle, DNS) is the Favre-averaged velocity (continuum, KT-TFM). Likewise, the appropriate quantity to compare to the particle velocity variance is the difference between the Favre-averaged kinetic energy and the Favre-averaged velocity squared, a measure of the fluctuating kinetic energy. As before, the fluctuating kinetic energy is time averaged and non-dimensionalized into a Reynolds number:

(4.12a,b ) $$\begin{eqnarray}Re_{\unicode[STIX]{x1D70E},x}=\frac{\unicode[STIX]{x1D70C}_{f}d_{p}}{\unicode[STIX]{x1D707}_{f}}\overline{[\langle \langle u_{s}^{2}\rangle \rangle +\langle \langle T\rangle \rangle -\langle \langle u_{s}\rangle \rangle ^{2}]^{1/2}},\quad Re_{\unicode[STIX]{x1D70E},y}=\frac{\unicode[STIX]{x1D70C}_{f}d_{p}}{\unicode[STIX]{x1D707}_{f}}\overline{[\langle \langle v_{s}^{2}\rangle \rangle +\langle \langle T\rangle \rangle -\langle \langle v_{s}\rangle \rangle ^{2}]^{1/2}}.\end{eqnarray}$$

The total fluctuating kinetic energy is the sum of two contributions: (i) a modelled ‘microscopic’ component arising from the granular temperature, $\langle \langle T\rangle \rangle$ , and (ii) a resolved ‘macroscopic’ component arising from variations in the mean velocity, $\langle \langle u_{s}^{2}\rangle \rangle -\langle \langle u_{s}\rangle \rangle ^{2}$ . Moving forward, it is worthwhile to note that (4.12) is an approximate comparison due to the directional decomposition, although the assumption of isotropy only extends to the microscopic contribution.

Figure 8. Comparison of KT-TFM predictions (thick lines) to DNS data (symbols) of the Reynolds number based on the standard deviation of the solids velocity in the $x$ - (black, squares) and $y$ -directions (red, circles) for $Ar=71$ and (a) $\langle \unicode[STIX]{x1D719}\rangle =0.10$ , (b) $\langle \unicode[STIX]{x1D719}\rangle =0.15$ , (c) $\langle \unicode[STIX]{x1D719}\rangle =0.25$ and (d) $\langle \unicode[STIX]{x1D719}\rangle =0.40$ . Black dashed lines show homogeneous (stable) solutions, $Re_{T}^{({\mathcal{H}})}$ . Right-hand side of vertical red lines demarcates high- $St_{T}$ region of (4.4).

Comparisons of $Re_{\unicode[STIX]{x1D70E}}$ between the DNS data and the KT-TFM simulations are provided in figures 8 and 9 for the elastic cases $Ar=71$ and 1432, respectively. Error bars of the DNS data again show the standard deviation of the five replicates, which are considerably larger here than in figures 6 and 7 since $Re_{\unicode[STIX]{x1D70E}}$ is based on second-order statistics (particle velocity variance). Both the streamwise, $Re_{\unicode[STIX]{x1D70E},y}$ (red circles), and transverse, $Re_{\unicode[STIX]{x1D70E},x}$ (black squares), directions are shown for comparison. In all cases, DNS and KT-TFM results are insensitive to the which transverse direction is considered, i.e. $Re_{\unicode[STIX]{x1D70E},x}\approx Re_{\unicode[STIX]{x1D70E},z}$ , and therefore only one is provided in figures 8 and 9 for clarity. The similarity of the statistics in the $x$ - and $z$ -directions is not surprising, given that the system is identical in both transverse directions. The homogeneous solution is also provided for comparison, which, unlike $Re_{S}$ , has a functional dependence on $\unicode[STIX]{x1D70C}^{\ast }$ . Since the homogeneous (stable) state is constant and uniform in space and time, $\langle \langle u_{s}^{2}\rangle \rangle =\langle \langle u_{s}\rangle \rangle ^{2}$ and $\langle \langle v_{s}^{2}\rangle \rangle =\langle \langle v_{s}\rangle \rangle ^{2}$ and, hence, (4.12) reduces to the thermal Reynolds number, i.e. $Re_{\unicode[STIX]{x1D70E}}^{({\mathcal{H}})}=Re_{T}^{({\mathcal{H}})}$ . As before, the vertical red lines represent the high- $St_{T}$ criterion of (4.4). Again, good agreement is generally observed where theory is expected to be valid and the agreement deteriorates moving to the left of the red line. Two notable exceptions are the comparisons at $Ar=71$ and $\langle \unicode[STIX]{x1D719}\rangle =0.10$ and 0.15 where good agreement is also observed to the left of the red line except at very low $\unicode[STIX]{x1D70C}^{\ast }$ where the KT-TFM becomes homogeneous; such agreement beyond the expected region of validity in these instances may be fortuitous. All cases show that the streamwise fluctuations ( $y$ -direction) are significantly larger than the transverse fluctuations ( $x$ - and  $z$ -directions), which vary only slightly from $Re_{T}^{({\mathcal{H}})}$ , at least where good agreement with the DNS data is observed. The good agreement of the KT-TFM predictions of $Re_{\unicode[STIX]{x1D70E},y}$ and $Re_{\unicode[STIX]{x1D70E},x}$ with the DNS data supports the assumption of isotropic microscale fluctuations, i.e. the primary cause for the difference in magnitudes between $Re_{\unicode[STIX]{x1D70E},y}$ and $Re_{\unicode[STIX]{x1D70E},x}$ is due to fluctuating kinetic energy generated by macroscopic (resolved) variations in the mean flow rather than microscopic (modelled) granular temperature. Although these comparisons can only be considered approximate, the fluctuating kinetic energy is the greatest level of detail compared between DNS and highly resolved KT-TFM to date and the good agreement observed speaks to the accuracy of KT-TFM in predicting clustering.

Figure 9. Comparison of KT-TFM predictions (thick lines) to DNS data (symbols) of the Reynolds number based on the standard deviation of the solids velocity in the $x$ - (black, squares) and $y$ -directions (red, circles) for $Ar=1432$ and (a) $\langle \unicode[STIX]{x1D719}\rangle =0.10$ , (b) $\langle \unicode[STIX]{x1D719}\rangle =0.15$ , (c) $\langle \unicode[STIX]{x1D719}\rangle =0.25$ and (d) $\langle \unicode[STIX]{x1D719}\rangle =0.40$ . Black dashed lines show homogeneous (stable) solutions, $Re_{T}^{({\mathcal{H}})}$ . Right-hand side of vertical red lines demarcates high- $St_{T}$ region of (4.4).

Although attention here has been focused on the high- $St_{T}$ region, where the KT-TFM predictions are expected to be valid, one consistent trend in the KT-TFM data in the lower- $St_{T}$ region is worth noting. Except for where the KT-TFM is at or approaching the homogeneous condition, i.e. $Re_{S}\rightarrow Re_{S}^{({\mathcal{H}})}$ and $Re_{\unicode[STIX]{x1D70E}}\rightarrow Re_{T}^{({\mathcal{H}})}$ , there is a direct correlation between the over-prediction of $Re_{S}$ and $Re_{\unicode[STIX]{x1D70E}}$ . Although fluctuating kinetic energy is the sum of mean flow and granular temperature components, deviations in $T$ from $T^{({\mathcal{H}})}$ are predominantly due to viscous heating, $\unicode[STIX]{x1D748}_{s}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{U}_{s}$ in (2.5), which is an inherently inhomogeneous term (if the system is near the homogeneous condition, $U_{s}$ is uniform and there is no viscous heating). Therefore, $Re_{\unicode[STIX]{x1D70E}}$ and in particular its deviation from $Re_{T}^{({\mathcal{H}})}$ can serve as an indication of the degree of inhomogeneity in the system. Consequently, the correlation between the over-prediction of $Re_{\unicode[STIX]{x1D70E}}$ and $Re_{S}$ suggests that the KT-TFM is over-predicting the degree of segregation (inhomogeneity), or ‘over-clustering’, which is producing larger fluctuating energy levels via larger variations in mean flow and consequently leading to an increased mean velocity via enhanced flow bypass. The enhanced segregation of the KT-TFM has been previously observed by Desjardins, Fox & Villedieu (Reference Desjardins, Fox and Villedieu2008) in the context of preferential concentration in frozen fluid-phase turbulence. Furthermore, it has been suggested that over-segregation of KT-TFM predictions is not only limited to the dilute, low mean-flow Stokes number regions (Passalacqua et al. Reference Passalacqua, Fox, Garg and Subramaniam2010). We conclude by noting that while the issue of over-segregation may contribute to disagreement between the KT-TFM predictions and the DNS data near the validity criteria (red line), good agreement is generally observed moving away from the estimated point of critical validity.

4.4 Regime boundary

As evidenced by three specific instances in § 3, the qualitative agreement between the KT-TFM and DNS is quite remarkable. Here, we attempt to compare the structure of the suspension in a more quantitative manner by calculating the transition between the chaotic and (2-D) plug regimes which was observed at both high and low $Ar$ numbers and all concentrations considered here. The transition will be characterized by a critical density ratio, $\unicode[STIX]{x1D70C}_{c}^{\ast }$ , for both KT-TFM and DNS.

As mentioned in § 3, the transition between the chaotic and plug regimes is gradual, which makes the identification of a single transition value, i.e. $\unicode[STIX]{x1D70C}_{c}^{\ast }$ , complicated. Figure 2 shows that plug-like structures form even in the chaotic region, at least temporarily. However, the two regimes may be distinguished by whether or not the fluid has enough inertia to break up the plugs. Therefore, $\unicode[STIX]{x1D70C}_{c}^{\ast }$ may be equivalently identified as the density ratio at which plugs persist (rather than simply exist).

In DNS, persistent plugs are identified using the structure factor. The structure factor is the Fourier transform of the particle pair distribution and can be calculated using

(4.13) $$\begin{eqnarray}S(\unicode[STIX]{x1D73F})=\frac{1}{N}\mathop{\sum }_{i}\mathop{\sum }_{j}\text{e}^{-\unicode[STIX]{x1D73F}\boldsymbol{\cdot }\boldsymbol{r}_{ij}\sqrt{-1}}.\end{eqnarray}$$

Here, $\boldsymbol{r}_{ij}$ is the vector from the centre of particle $i$ to particle $j$ , and $\unicode[STIX]{x1D73F}$ is the wave vector. Excess or deficit of $S(\unicode[STIX]{x1D73F})$ from its equilibrium value of unity is an indication that the spatial distribution of particles has a Fourier component in the direction of $\unicode[STIX]{x1D73F}$ with wavelength $2\unicode[STIX]{x03C0}|\unicode[STIX]{x1D73F}|^{-1}$ . In a periodic domain, the number of structures must be a positive integer. The first, second, and third structure factors along the direction of sedimentation are therefore given by

(4.14) $$\begin{eqnarray}S_{y}^{\{k\}}=\frac{1}{N}\mathop{\sum }_{i}\mathop{\sum }_{j}\text{e}^{-(2k\unicode[STIX]{x03C0}/L_{y}^{\ast })(y_{i}^{\ast }-y_{j}^{\ast })\sqrt{-1}},\end{eqnarray}$$

where $k=1$ , 2 and 3, respectively. Persistent plugs are identified by the first, second or third structure factors attaining a value greater than 25 for 50 % of the simulation time. Figure 10 shows the evolution of $S_{y}^{\{k\}}$ , obtained from a system with $Ar=1432$ , $\unicode[STIX]{x1D70C}^{\ast }=25$ , $\langle \unicode[STIX]{x1D719}\rangle =0.25$ and $e=1.0$ . At time A, the highest structure factor was approximately 8, and the system only showed a slight heterogeneity that cannot be clearly distinguished from statistical fluctuations. At time B, the second structure factor reached 17, and heterogeneity started to become noticeable. After time B, the heterogeneity continues to increase and at time C two horizontal clusters can be clearly identified when the second structure factor has a value of 36. As this simulation has strong clusters for more than 50 % of the time, it is identified as one with persistent plugs.

Figure 10. The graph on the far right shows the evolution of the first three structure factors in the $y$ -direction for a simulation with $Ar=1432$ , $\unicode[STIX]{x1D70C}^{\ast }=25$ , $\langle \unicode[STIX]{x1D719}\rangle =0.25$ and $e=1.0$ . In this case, the second structure factor (dotted line) increases significantly as the simulation progresses. Visualization of distributions of particle centres (left three graphs) at three different times (A, B and C) shows heterogeneities that grow with time, and specifically the two-cluster structure at time C.

The structure factor is calculated based on particle locations and is therefore not directly applicable to continuum (KT-TFM) results. A successful regime indicator was found by examining the temporal signatures, which provide a measure of the dynamics of the system rather than the spatial structure of the system. Namely, the 2-D and 1-D plug regimes manifest as travelling wave limit cycles which are cyclical and regular, while the chaotic regime is dynamic and irregular. The two types of system dynamics are differentiated with a type of auto-correlation. First, the near-homogeneous cases are removed by excluding cases with $\unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{max}<0.5$ . With the remaining inhomogeneous cases, the solids concentration field is decomposed into $N_{xyz}=N_{x}\times N_{y}\times N_{z}$ time signals denoted $\unicode[STIX]{x1D719}_{ijk}(t^{\ast })=\unicode[STIX]{x1D719}(x^{\ast }=x_{i}^{\ast },y^{\ast }=y_{j}^{\ast },z^{\ast }=z_{k}^{\ast },t^{\ast })$ , where $x_{i}^{\ast }$ , $y_{j}^{\ast }$ and $z_{k}^{\ast }$ are the $N_{xyz}$ discrete grid locations. Only the time in the averaging window, $[t_{1}^{\ast },t_{2}^{\ast }]$ , is considered. In the middle of each $\unicode[STIX]{x1D719}_{ijk}(t^{\ast })$ string, a time unit sample of length $\unicode[STIX]{x1D6FF}t\ast =200$ and 50 for $Ar=71$ and 1432, respectively, is selected. The sample is cross-correlated with the entire $\unicode[STIX]{x1D719}_{ijk}(t^{\ast })$ signal and by moving the sample in time by $\pm n\unicode[STIX]{x1D6E5}_{t}^{\ast }$ , where $n$ is a positive integer, The procedure is terminated when the sample reaches the beginning and end of the full signal. The $N_{xyz}$ temporal cross-correlations are then averaged over all spatial locations. When the sample is compared with itself, i.e. the auto-correlation at $n=0$ , the averaged cross-correlation is unity. As the sample moved to the left or right (forward or backward in time) the averaged cross-correlation decays. If the signal is perfectly cyclic – indicative of an exact travelling wave limit cycle – the cross-correlation will increase to unity again at secondary peaks indicating an exact correlation in the signals. The averaged cross-correlation will remain near zero if the two signals show no correlation (similarity) whatsoever. A simulation is considered to be in a plug-flow regime if a secondary peak (i.e. outside of the auto-correlation window) of the averaged cross-correlation attains a value of $2/3$ . Conversely, a simulation is considered to be in the chaotic regime if no secondary peak attains a value of $1/3$ (i.e. all secondary peaks ${<}1/3$ ). This treatment will yield a range of $\unicode[STIX]{x1D70C}_{c}^{\ast }$ rather than a single value, which is representative of the smooth transition between the chaotic and plug regimes that $\unicode[STIX]{x1D70C}_{c}^{\ast }$ is meant to represent. The two criteria signifying the transition region between chaotic and plug regimes are represented as upper and lower error bars centred about their geometric mean.

Figure 11. Comparison of continuum predictions (thick lines with small, open symbols) to DNS data (filled symbols) of the critical density ratio for the transition of dynamic to plug-flow patterns at $Ar=71$ (black squares) and $Ar=1432$ (blue circles). The KT-TFM points show the geometric mean of the criteria for chaotic and plug regimes which are indicated by the lower and upper error bars, respectively. Dashed red and green lines indicate high- $St_{T}$ criteria of (4.4) for $Ar=71$ and 1432, respectively.

A comparison of the chaotic-to-plug transition criteria determined from the DNS and KT-TFM simulations is provided in figure 11. Except at the highest concentration, $\langle \unicode[STIX]{x1D719}\rangle =0.40$ , the comparisons are very favourable. It should be pointed out that the seemingly good agreement for $\langle \unicode[STIX]{x1D719}\rangle =0.40$ at $Ar=71$ is apparently fortuitous as the error bar separates a single pair of simulations which, unlike all other cases, transition directly from homogeneous to plug, skipping the chaotic regime entirely. It should also be mentioned that most cases fall just below the required thermal Stokes number criterion which may be responsible for the discrepancy. However, the good agreement in figure 11 is further evidence that the qualitative similarities between KT-TFM and DNS extend beyond the three select cases observed previously in § 3.

4.5 Knudsen number analysis

Although it has been demonstrated that the breakdown of the high- $St_{T}$ condition can explain the poor agreement between KT-TFM predictions and the DNS data at low density ratios, it is worthwhile to analyse another assumption used in the derivation of the KT-TFM: the low Knudsen number assumption. Recall that this assumption stems from the Chapman–Enskog expansion about a small parameter, namely the Knudsen number defined as $Kn\equiv \ell /L$ , where $\ell$ is the mean free path defined previously in § 4.1 and $L$ is a macroscopic length scale characterizing the local gradient, i.e. $L$ is small for large gradients and vice versa. Accordingly, the low- $Kn$ assumption is also referred to as the ‘small gradient’ assumption. A priori, this assumption appears at odds with the presence of clustering instabilities, since a relatively large concentration gradient is present normal to the cluster interface. This question motivates the subsequent analysis.

Following a previous work (Mitrano et al. Reference Mitrano, Zenk, Benyahia, Galvin, Dahl and Hrenya2014), we take $L$ to be the length at which a variable of interest changes by 20 %, or $L_{f}=0.2|f|/|\unicode[STIX]{x1D735}f|$ , where $f$ is a variable of interest. Here, we consider the solids concentration, the vertical component of solids velocity and the granular temperature, i.e. $f=\unicode[STIX]{x1D719}$ , $v_{s}$ and $T$ . Special care is taken with the solids velocity. In the current system, the $v_{s}$ is only defined relative to the frame of reference. Hence, $Kn_{v}$ should be normalized with the thermal speed, $c$ , in order to eliminate the possibility of calculating a large value $Kn$ simply due to $v_{s}$ becoming (locally) null (Garzó & Santos Reference Garzó and Santos2003). The three Knudsen numbers of interest are then:

(4.15a-c ) $$\begin{eqnarray}Kn_{\unicode[STIX]{x1D719}}=\frac{5}{6\sqrt{2}}\frac{|\unicode[STIX]{x1D735}^{\ast }\unicode[STIX]{x1D719}|}{\unicode[STIX]{x1D719}^{2}\unicode[STIX]{x1D712}},\quad Kn_{v}=\frac{5}{12}\frac{|\unicode[STIX]{x1D735}^{\ast }v_{s}|}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}\sqrt{T}},\quad \text{and}\quad Kn_{T}=\frac{5}{6\sqrt{2}}\frac{|\unicode[STIX]{x1D735}^{\ast }T|}{\unicode[STIX]{x1D719}\unicode[STIX]{x1D712}T},\end{eqnarray}$$

where $\unicode[STIX]{x1D735}^{\ast }$ is the gradient operator non-dimensionalized by particle diameter.

Figure 12. Instantaneous Knudsen numbers in a KT-TFM simulation at $Ar=71$ , $\langle \unicode[STIX]{x1D719}\rangle =0.15$ , and $\unicode[STIX]{x1D70C}^{\ast }=100$ and $e=1.0$ . A slice in the $xy$ -plane shows contours of (a) $Kn_{\unicode[STIX]{x1D719}}$ , (b) $Kn_{v}$ and (c) $Kn_{T}$ . The cumulative distribution function of each is provided in (d) considering the entire domain.

A typical example of the three Knudsen numbers defined by (4.15) is shown in figure 12 at a particular instant in a fully developed, statistically steady state. The contour plots and the cumulative distribution functions (CDF) all use a logarithmic scale. It is observed that all three Knudsen numbers span several orders of magnitude from approximately $O(10^{-2})$ to $O(10^{2})$ . Hence, calculating a simple $Kn$ mean value would be misleading, tending to skew the results to very large $Kn$ values. To avoid such bias, the ‘reverse’ analysis is performed: rather than seeing if the mean values exceed a given criterion, we seek to find what percentage of the domain exceeds the specified criterion. Similar to (4.4), $Kn$ is assumed to be sufficiently small if it is less than a negative half order of magnitude, i.e. $Kn<10^{-1/2}$ . For example, at the time shown in figure 12 approximately 80 % of the domain is violating the small- $Kn_{\unicode[STIX]{x1D719}}$ assumption, almost 90 % of the domain is violating the small- $Kn_{T}$ assumption, and nearly the entire domain is violating the small- $Kn_{v}$ assumption.

Figure 13. Time-averaged values of the percentage of the KT-TFM simulation domain with $Kn\geqslant 10^{-1/2}$ for $Ar=71$ and $\langle \unicode[STIX]{x1D719}\rangle =$ (a) 0.10, (b) 0.15, (c) 0.25 and (d) 0.40. $Kn_{\unicode[STIX]{x1D719}}$ , $Kn_{v}$ and $Kn_{T}$ are indicated by black, red and blue lines with open square, circle and diamond symbols, respectively.

Figure 14. Time-averaged values of the percentage of the KT-TFM simulation domain with $Kn\geqslant 10^{-1/2}$ for $Ar=1432$ and $\langle \unicode[STIX]{x1D719}\rangle =$  (a) 0.10, (b) 0.15, (c) 0.25 and (d) 0.40. $Kn_{\unicode[STIX]{x1D719}}$ , $Kn_{v}$ and $Kn_{T}$ are indicated by black, red and blue lines with square, circle and diamond symbols, respectively.

The percentage of the domain that violates the small- $Kn$ assumption, i.e. fractional region where $Kn\geqslant 10^{-1/2}$ , is time averaged and presented in figures 13 and 14 for the elastic cases $Ar=71$ and 1432, respectively. Two trends are readily apparent: (i) large portions of the domain are violating the low- $Kn$ assumption for all three variables of interest at both $Ar$ conditions and all mean concentrations considered and (ii) the velocity Knudsen number is almost always the largest of the three. The latter appears to be a result of the concentration field being out of phase with the solids velocity gradient, unlike the concentration and temperature gradients. In other words, where the concentration field exhibits a local minima, the gradients of concentration and temperature are also small and the two effects tend to counteract each other. The velocity gradients, on the other, are not typically small near concentration minima producing larger $Kn_{v}$ relative to $Kn_{\unicode[STIX]{x1D719}}$ and $Kn_{T}$ . At low $\unicode[STIX]{x1D70C}^{\ast }$ for $Ar=71$ , where some of the largest discrepancy between the KT-TFM and DNS results was observed due to violation of the high- $St_{T}$ assumption, all three Knudsen numbers are very small due to the near-homogeneous condition of the KT-TFM simulations. Conversely, in some regions where good agreement with DNS data was noted (high- $St_{T}$ regions), wholesale violations of the small- $Kn$ assumption are observed. The results of this analysis indicate that the small- $Kn$ assumption is not the culprit for the breakdown in the agreement with the DNS data, unlike the high- $St_{T}$ assumption. It is worthwhile to note that a similar observation (significant violation of $Kn$ assumption despite good agreement) was also reported in a previous work for the onset of clustering in the granular HCS (Mitrano et al. Reference Mitrano, Zenk, Benyahia, Galvin, Dahl and Hrenya2014) where clustering was caused by particle inelasticity alone. The ability of the KT-TFM to perform well outside of its intended $Kn$ -range of validity is reminiscent of molecular gases, i.e. Navier–Stokes, and is very fortunate as higher-order theories such as linear-Burnett, Burnett, super-Burnett, etc. are unwieldy, particularly for engineering applications of particulate flows. Alternatively, moments of the starting kinetic equation can be solved directly, though unlike the Chapman Enskog method, such numerical solutions do not result in analytical expressions for the constitutive relations; see, for example, Desjardins et al. (Reference Desjardins, Fox and Villedieu2008).

4.6 Effect of inelasticity

All of the results presented thus far have been elastic, i.e. $e=1$ . Inelastic simulations, both DNS and KT-TFM, were also carried out for all cases with $e=0.9$ . For the sake of brevity, the results are not presented here since inelasticity did not have a significant impact on the systems examined. While $e$ has a widespread effect on the KT-TFM governing equations, its impact on the dissipation of granular temperature through the zeroth-order cooling rate, $\unicode[STIX]{x1D701}_{0}$ , see (A 21), may be the most important due to the role of granular temperature dissipation as a clustering mechanism (Fullmer & Hrenya Reference Fullmer and Hrenya2017b ). The contribution of inelasticity to the total dissipation (inelastic collisions $+$ thermal drag) is quantified in figure 15 for the homogeneous state. It is seen that inelasticity plays a relatively minor role in the total dissipation for the lower $Ar=71$ case throughout the $\unicode[STIX]{x1D70C}^{\ast }$ range. Only at $Ar=1432$ and the highest density ratios does inelasticity contribute at a similar level to the total dissipation as the thermal drag – which is also the only region where noticeable differences with the elastic simulations were observed.

Figure 15. (a) Ratio of thermal drag to total granular temperature dissipation (inelastic collisions $+$ thermal drag) for $\langle \unicode[STIX]{x1D719}\rangle =0.10$ (black), 0.15 (red), 0.25 (blue) and 0.40 (green) and $Ar=71$ (lower, solid lines) and $Ar=1432$ (upper, dashed lines). (b) Sedimenting Reynolds number at $Ar=1432$ , $\unicode[STIX]{x1D70C}^{\ast }=1000$ and $e=1.0$ (black) or $e=0.9$ (red) from DNS (symbols) and KT-TFM (thick lines) simulations. Dashed line show homogeneous (stable) solution, $Re_{S}^{({\mathcal{H}})}$ .

Figure 15 also shows the effect of inelasticity on the sedimenting Reynolds number for the case most sensitive to $e$ , i.e. $Ar=1432$ and $\unicode[STIX]{x1D70C}^{\ast }=1000$ . The DNS data show a decrease in $Re_{S}$ for all concentrations – a trend that is matched in the KT-TFM predictions. Quantitatively, the inelastic KT-TFM predictions agree with the DNS data slightly more favourably than the elastic results. The decrease of $Re_{S}$ has come about by a change in structure that can be best understood when viewed in terms of the most important consequence of particle inelasticity: a source of clustering. In general, clustering in gas–solid flows is induced by either a difference in mean relative motion and/or the dissipation of granular temperature (Fullmer & Hrenya Reference Fullmer and Hrenya2017b ). Due to thermal drag (which leads to a dissipation of granular temperature), both sources of clustering are present in the elastic cases, so the inclusion of particle inelasticity does not change the stability pattern of the system appreciably. However, the inelastic cases always have more dissipation and are therefore more unstable than the elastic cases which has the effect of moving the flow pattern regime boundaries to lower density ratios. One can imagine that the $Re_{S}$ curves of figure 7 have been shifted to the left (but a non-uniform shift due to the diminishing impact of $\unicode[STIX]{x1D701}_{0}$ with decreasing $\unicode[STIX]{x1D70C}^{\ast }$ shown in figure 15 a). Therefore, while the flow structure for all cases presented in figure 15(b) are in the 1-D plug regime, the amplitude of the plug is still larger in the inelastic cases. In closing, it is worth noting that the minor impact of inelasticity observed here should not be extrapolated to highly inelastic conditions or, perhaps more importantly, to higher mean-flow Stokes number conditions (i.e. higher $Ar$ and/or  $\unicode[STIX]{x1D70C}^{\ast }$ ).