Verification of Eulerian-Eulerian and Eulerian-Lagrangian simulations for particle-laden vertical channel flow

Particle-laden flows in a vertical channel were simulated using an Eulerian–Eulerian, Anisotropic-Gaussian (EE-AG) model. Two sets of cases varying the overall mass loading were done using particle sizes corresponding to either a large or small Stokes number. Primary and turbulent statistics were extracted from these results and compared with counterparts collected from Eulerian–Lagrangian (EL) simulations. The statistics collected from the small Stokes number particle cases correspond well between the two models, with the EE-AG model replicating the transition observed using the EL model from shear-induced turbulence to relaminarization to cluster-induced turbulence as the mass loading increased. The EE-AG model was able to capture the behavior of the EL simulations only at the largest particle concentrations using the large Stokes particles. This is due to the limitations involved with employing a particle-phase Eulerian model (as opposed to a Lagrangian representation) for a spatially intermittent system that has a low particle number concentration.


Introduction
Particle-laden flows are pervasive throughout the natural and industrial world. Modeling particle-laden flows is of major interest in developing a better understanding of these systems and, in an industrial context, improving process optimization, troubleshooting, and scale-up. There have been a variety of approaches to modeling these systems including direct numerical simulation (DNS), Eulerian-Lagrangian (EL) simulation, and Eulerian-Eulerian (EE) simulation (Pan et al., 2002;Feng & Musong, 2014;Capecelatro & Desjardins, 2013;Igci & Sundaresan, 2011;Chen et al., 2012;Hrenya & Sinclair, 1997;Schneiderbauer, 2018). One major challenge in doing so involves not only successfully resolving the features present, but also doing so in a way that is computationally feasible and efficient (Fox, 2012;van der Hoef et al., 2006). With that in mind, the purpose of this study is to evaluate the ability of an EE-anisotropic Gaussian (EE-AG) model (Levermore & Morokoff, 1996;Vié et al., 2015;Kong et al., 2017) to reproduce the results of particleladen channel flow modeled using an EL model (Capecelatro & Desjardins, 2013). Turbulent statistics and terms for both the particle and fluid phases corresponding to the Reynolds-averaged Navier-Stokes equations derived in (RO, 2014) are extracted from the results for different mass loadings and compared. Accurately obtaining these turbulent statistics is necessary for the future development of turbulent transport parameters (e.g., turbulent viscosity) in turbulence models based on Reynolds averaging. In prior work (Kong et al., 2017;Kasbaoui et al., 2019), the applicability of the EE-AG model to reproduce EL turbulence statistics for statistically homogeneous particle-laden flows with cluster-induced turbulence (CIT) has been demonstrated. In comparison, the turbulent channel flow investigated in this work is considerably more challenging due to the walls (Capecelatro et al., 2018), resulting in a spatial dependence of the turbulence statistics in the wall-normal direction.
EL simulations involve the tracking of individual particles while treating the fluid phase surrounding the particles as a continuum (Capecelatro & Desjardins, 2013). Special modeling considerations must be taken when using an EL model accounting for the volume of the particles in the fluid, inter-particle collisions, and interphase interactions between the particles and fluid (Capecelatro & Desjardins, 2013;Patankar & Joseph, 2001;Capecelatro et al., 2014). EL methods have successfully modeled a wide variety of particleladen system including fluidized beds, vertical channel risers, sprays, and impinging jets (Capecelatro & Desjardins, 2013). One practical limitation to the use of EL models is the large computational cost involved in systems containing high concentrations of particles due to the shear number of individual particles that must be tracked. Furthermore, parallelization issues can arise in systems which feature a significant amount of maldistribution in the particle phase such as in CIT (Kong et al., 2017). In those cases the bulk of the computational load is imbalanced towards the particular processors responsible for solving the governing equations in the portions of the system geometry that include the clusters. Most importantly, tracking and obtaining accurate turbulent statistics from EL data requires strenuous post-processing (Capecelatro et al., 2015;Capecelatro et al., 2018). This is especially true for the particle phase, which inherently has problems in areas such as in near-wall regions where there are consistently not enough particles present over a reasonable amount of modeling time to obtain accurate statistics.
By treating the particle phase as another continuous phase, EE simulations are able to alleviate some of the drawbacks present in EL models (Drew, 1983). The primary trade-offs are the modeling considerations that must be addressed to handle the effects of physical phenomena such as particle collisions and interphase drag. Like EL simulations, EE simulations have been successfully implemented to model a variety of particleladen flows. A key assumption in previous models with an Eulerian particle phase is that the particle-phase velocity distribution is nearly Maxwellian (Kong et al., 2017). This assumption is most valid for near-zero Knudsen number systems such as fluidized beds where particles are densely packed and highly collisional.
The assumption makes it so that physical phenomena such as significant particle-phase velocity anisotropy and particle trajectory crossing cannot be accounted for (Desjardins et al., 2008). Failing to capture this results in the spontaneous development of clusters even in highly dilute systems that are not replicated in EL modeling (Kong et al., 2017;Kasbaoui et al., 2019). One method that has been developed to remedy this issue that has been demonstrated in previous work involves considering the particle-phase velocity fluctuations as an anisotropic Gaussian distribution (Kong et al., 2017;Kasbaoui et al., 2019). This EE-AG model is used in this study to compare with EL simulation results.
In this work, we consider a vertical particle-laden channel flow spanning from dilute to moderately dense in terms of the number concentration. Our previous work (Capecelatro et al., 2018) demonstrated the transition from pure-fluid turbulent channel flow to a state of CIT as more particles are added to increase the mass loading. CIT occurs when the turbulence in the fluid phase is governed primarily by the presence of clusters in the particle phase (Agrawal et al., 2001). These clusters first originate from interphase momentum coupling resulting from any significant mean velocity difference between the fluid and particle phases and are sustained by the CIT in the fluid phase (Capecelatro et al., 2015;Kasbaoui et al., 2019). Prior to this work, most studies on wall-bounded particle-laden flows involve only dilute systems where CIT is absent due to weak interphase momentum coupling. In those studies, the fluid-phase turbulence is governed primarily though mean-shear production as what occurs in pure-fluid flow in a channel at high Reynolds number (Kim et al., 1987). In our prior work on particle-laden channel flow at significant mass loading (Capecelatro et al., 2016;Capecelatro et al., 2016), the channel Reynolds number was too low to observe transitions between CIT and fully developed turbulence. Here, our primary interest is whether the EE-AG model can reproduce the turbulence statistics found in the EL simulations for the same high-Reynolds-number channel flow investigated in (Capecelatro et al., 2018). For such cases, the computational cost is significant due to the grid requirements needed to resolved the turbulent boundary layers at the walls (Kim et al., 1987).
The remainder of this work is organized as follows. In Sec. , we briefly review the EL and EE-AG models developed in our prior work. More details can be found in the literature (Capecelatro & Desjardins, 2013;Kong et al., 2017). Of particular importance are the near-wall boundary conditions for the particle-phase moments described in Sec. . In Sec. , the EE-AG simulation algorithm implemented in OpenFOAM, an open-sourced CFD package (Weller et al., 1998), is described for the vertical channel flow case investigated in (Capecelatro et al., 2018) where details on the EL simulations can be found. The detailed comparison of the turbulence statistics found from the two different modeling approaches is presented in Sec. . Conclusions are drawn in Sec. .

Governing equations and methods
This section briefly summarizes the equations used in the EL and EE-AG simulations applied in this study. The EL model is solved in the NGA environment (Capecelatro & Desjardins, 2013;Desjardins et al., 2008), while the EE-AG model is solved in the OpenFOAM environment (Kong et al., 2017).

Eulerian-Lagrangian model
The fluid-phase equations for the EL model are based on the conservation of mass and momentum for a continuous fluid phase derived from a volume-filtering approach originally developed by Anderson and Jackson (Anderson & Jackson, 1967). With the application of a volume-filtering with a kernel H(|x|) over a volume of fluid V f (1) any point variable represented as a such as phase velocity, pressure, or force can be replaced by a smoothed, locally-filtered volume field (Capecelatro & Desjardins, 2013). The fluid-phase continuity equation using these filtered variables is where α f is the fluid-phase volume fraction and U f is the filtered fluid-phase velocity. The fluid phase is assumed to be incompressible with a constant density of ρ f . The conservation of momentum equation is where τ f is the total stress tensor, g is the gravitational acceleration, F U is the subgrid stress, F µ is the residual viscous stress, and F interphase is the interphase momentum exchange. The viscous component of the total stress tensor is modeled using a gradient-viscosity model, where P f is the fluid-phase pressure and ν f,ef f ective is the effective fluid-phase viscosity (Kong et al., 2017). The subgrid stress term in the CIT regime is often neglected (Patel et al., 2017). The effect of the residual viscous stress term is rolled up into the effective viscosity using the fluidized bed model from Gibilaro et al., , where ν f is the bulk viscosity for the fluid phase (Gibilaro et al., 2007).
The interphase momentum exchange term is developed through putting the total force of the fluid acting upon each particle through the volume-filtering kernel function, interphase is the force of the fluid acting upon a single particle i and n p is the total number of particles (Patel et al., 2017). This is approximated by where V p is the particle volume and f (i) drag is the drag acting upon a single particle i. The latter is modeled using a Stokes drag law p are the mass and velocity of a particle i respectively and τ p is the Stokes drag timescale. The latter is defined as τ p = d 2 p ρ p /(18ν f ρ f ) where d p is the mean particle diameter (Passalacqua et al., 2010).
For the particle phase the displacement of each particle is calculated using a balance based on Newton's second law of motion with interphase, collisional, and gravitational forces considered where f (i) collisions is the collisional source term (Capecelatro & Desjardins, 2013). Two types of particleparticle interactions are applied in this model, normal is the normal collision force and f (i,j) f riction is the tangential collision force involving interparticle friction between a particle i and another particle j. The Cundall and Strack soft-sphere approach is applied for the normal force term where k n is the normal spring stiffness parameter, δ (i,j) n is the overlap between particles i and j, n (i,j) n is the normal unit vector between the two particles, η n is the normal dampening parameter, U (i,j) p,n is the relative normal velocity between the two particles, d (i,j) is the distance between the centers of the two particles, r (i) is the radius of a particle, and λ is the force range parameter (Cundall & Strack, 1979). The van der Hoef et al. model (van der Hoef et al., 2006) is applied for the frictional model where k t is the tangential spring stiffness parameter, δ t is the tangential displacement, η t is the tangential dampening parameter, U (i,j) p,t is the relative tangential velocity between the two particles, and µ t is the tangential friction coefficient.

Eulerian-Eulerian anisotropic Gaussian model
As in the EL model, the fluid-phase equations for the EE-AG model uses the conservation of mass to solve for the fluid-phase volume fraction and velocity profiles respectively where U f and U p are the fluid and particle velocity fields respectively. The total stress tensor, τ f , uses an identical gradient viscosity model as that applied in the EL model, where The Gibiliaro et al. model (Gibilaro et al., 2007) is likewise used for the effective viscosity. The main distinction between the EL and EE-AG in this equation is in the direct application of a Stokes drag law term to handle interphase interactions.
The equations for the particle phase solved in the EE-AG model are derived from the first three loworder moments of the particle velocity normal distribution function: et al., 2017). The zeroth-order moment transport equation is simply the equation for the conservation of particle mass where α p is the particle-phase volume fraction and ρ p is the bulk density of the particle phase.
The first-order moment transport equation is the conservation of particle momentum where P p is the particle-pressure tensor and G p is the collisional flux tensor. The latter is derived using Enskog-Boltzmann kinetic theory G p = 4 5 η c α p g 0 (3Θ p I + 2P p ) where η c = 1 2 (1 + e c ) is derived from e c , the collisional restitution coefficient, g 0 is the particle radial distribution function, and Θ p is the granular temperature (Marchisio & Fox, 2013). The particle-pressure tensor can be found inside the definition of the second-order particle-phase velocity moments.
The transport equation for the second-order moment can be rearranged into a transport equation for the second-order particle-pressure tensor: where k p,cond is the granular conductivity, τ c is the collisional timescale, and ∆ * is the collisional energy source term. The granular conductivity is modeled as k p,cond = ν p /P r p , where ν p = 2.0/(c c (1 + e c ) 2 ), is the particle-phase turbulent viscosity, P r p = (16 − 11e c )/(15 − 5e c ) is the particle-phase Prandtl number, and c c is the interparticle collision frequency coefficent (Capecelatro et al., 2016). Both the collisional timescale and collisional energy source term are developed from a linearized Bhatnagar-Gross-Krook (BGK) inelastic collision model (Bhatnagar et al., 1954). The BGK collisional energy timescale is τ c = √ πd p /(6α p Θ p ) and the BGK collisional energy source term is ∆ * = η 2 c Θ p I + (1 − η c ) 2 P p .

EE-AG wall boundary conditions
A no-slip boundary condition is applied at the walls in solving the fluid-phase momentum transport equation in (12). The inlet and outlet boundaries are cyclically periodic for all fluid and particle quantities. The pressure gradient in the fluid-phase momentum transport equation is adjusted to maintain an average fluidphase velocity over the inlet and outlet boundary at the channel velocity parameter specified in Table 1.
The wall boundaries for the particle-phase volume fraction, momentum, and pressure tensor in (13), (14), and (??) use a flux reflection model at the wall (Kong & Fox, 2017). In this model a flux is developed from the wall in opposition to the flux going to the wall taking into account the particular physical properties of a wall-particle collision present in a given system. The particle velocity moment fluxes of order n coming from this wall boundary, F w (M n p ), are calculated as where f r (M n p ) is the moment generation function for the particle velocity reflection distribution coming from the wall surface, S w , and n w is the wall-normal direction vector. The moment generating function is represented by contributions from the specular-reflection component, f r,s (M n p ), and the diffusive-reflection component, f r,d (M n p ), of particle velocity: φ s is the tunable specularity factor that can be adjusted to balance the contribution of both components of a particle hitting the wall. The specular-reflection component of particle-wall collisions is defined by where f i (M n p ) is the incoming velocity distribution at the wall and e w is the wall restitution coefficient. In this work, only the specular-reflection component of particle-wall collisions is considered and thus the specularity factor is set at zero and the diffusive-reflection component is neglected. For the zeroth-order particle velocity moment applicable to the particle continuity equation in (13) the flux reflection model results in a zero-gradient boundary condition.
TheIn addition to the flux refection model at the wall boundary, sub-particle-diameter flux blending is applied near the wall boundary for the particle volume fraction, momentum, and pressure. This is done through amalgamating the fluxes of a specified number of highly-refined cells close to the wall boundary (which are small relative to the particle diameter as seen in Fig. 1), effectively creating a single larger cell out of smaller ones. This method contrasts with that found in (Patel et al., 2019) such that the geometric integrity of the system is maintained rather than shifting the wall fluxes outside of highly refined cells. The impetus for this is the large gradient in particle pressure present in the near-wall region that is exacerbated by the small size of the near-wall cells. This results in a large gradient for the granular temperature in the near-wall region as seen in Fig. 1. While maintaining the original refinement resulted in better correspondence with EL data in the near-wall region in regard to the particle velocity, the large granular temperature gradient significantly limits the size of the timestep, especially for the large Stokes number particle cases. With the application of near-wall blending, the timestep increased by an order of magnitude for all mass loadings and particle sizes.

Simulation algorithm and parameters
The geometry used in the simulations is the rectangular vertical channel from (Capecelatro et al., 2018). The dimensions of the channel and their corresponding labels are detailed in Table 1. Both EE-AG and EL cases use an identical block grid as the solution mesh including near-wall refinement based on the work of Kim et al. (Kim et al., 1987). As discussed in (Capecelatro et al., 2018) and shown in Fig. 1, the near-wall grid spacing required to resolve the turbulent channel flow is significantly smaller than the particle diameter.
The solution algorithm for the EE-AG model is summarized in the following steps: 1. Initialization of all moments and variables.
where Ω is the computational domain, C CF L is the CFL number, andũ α ,ṽ α , andw α are the particlephase velocity abscissas.
3. Particle-phase moment fluxes are computed at all surfaces and are used to update the corresponding variables, α p , U p , and P p .
4. The updated U p and P p variables are corrected by the application of the collisional terms ∇ · ρ p α p G p and 2ρpαp τc (∆ * − P p ), respectively. 5. The updated U p and P p variables are further corrected by the particle-phase interphase drag and gravitational accelerations terms ρ p α p g + ρpαp τp (U f − U p ) and − 2ρpαp τp P p , respectively. 6. All particle-phase moments are updated with the corrected values of U p and P p . 7. The fluid-phase momentum transport equation from (12), U f , is constructed as a semi-discretized equation separating diagonal and off-diagonal elements, where λ f = (A + K drag ) −1 is the inverse of the sum of the diagonal coefficients A and the overall drag coefficient K drag = ρpαp τp and H is the off-diagonal contributions to the fluid-phase equation.
8. The fluid-phase pressure-gradient equation constructed using the fluid-phase velocity fluxes on each cell surface S, develops the fluid-phase pressure field and its result is applied to update the fluid-phase velocity field.
9. Iterate steps 7 and 8 until the fluid-phase pressure converges.
10. Advance in time and repeat from step 2 until the solution is complete. Table 1 contains the parameters for the fluid and particle phases used in the simulations (Capecelatro et al., 2018). The two particle diameters correspond to flow with a Stokes number of 1.74 for the smaller particle diameter, and 17.86 for the larger particle diameter. The overall concentration of particles for each case is adjusted to achieve a specific mass loading, defined by ϕ = ρ p α p /ρ f α f . The cases for the small Stokes number particles include mass loadings of 0.2, 1, 2, and 4. For the large Stokes number particles, the cases are 0.2, 1, 2, 4, and 10. The lower mass loading ceiling for the small Stokes number particle cases is a result of the computational cost for the EL simulation becoming too costly with the number of particles that must be tracked at larger mass loadings. The 0.2 mass loading case corresponds to a particle count of 1.47 million using small Stokes number particles versus 44.8 thousand using large Stokes particles. An additional particle-free pure-fluid case was also done. All wall-normal profiles displayed in this work are averaged spatially in both the spanwise and streamwise directions and temporally over 1 s of real time.

Pure-fluid case
For modeling pure-fluid channel flow, the mass conservation and Navier-Stokes equations in (2) for the EL model and (12) for the EE-AG model are identical. The only difference is that the EL model is solved using the NGA environment and the EE-AG is solved in OpenFOAM. Figure 2 compares the streamwise velocity of the EE-AG and EL models as well as the Reynolds-stress components in a pure-fluid channel. The pure-fluid mean velocity profiles of the EL and EE-AG models demonstrate good agreement. The Reynolds stresses, on the other hand, display minor deviations between the two simulation codes, especially towards the center of the channel. The EE-AG code shows more anisotropy towards the streamwise direction, with the wall-normal and spanwise turbulent components being more prominent in the EL code. The fluid kinetic energy of the EE-AG model also overpredicts that observed in the EL model. Nonetheless, the pure-fluid comparisons give us confidence that the two codes are capturing correctly the turbulent channel flow and, thus, can be employed for comparisons between the EE an EL approaches for particle-laden flows.  Fluid-phase turbulent kinetic energy is denoted as () and (), correlated particle-phase turbulent kinetic energy as (   Figures 3-6 show the mean velocities and Reynolds-stress components for all cases using particles corresponding to the smaller Stokes number. With the exception of the mass loading of 4 case in Fig. 6, there is good agreement between EL and EE-AG results for the fluid and particle velocities. The significant difference seen in the mass loading of 4 case is attributed to a particular parameter deviation in the execution of the EL model where a slower overall channel velocity of 4.24 (m/s) was applied rather than the intended value of 5.02 (m/s). The result was that the velocity profile of the EL is much more laminar in shape.

Small Stokes number particles
The lowest particle mass loading in Fig. 3 demonstrates the similar kind of deviations of the Reynolds stresses seen in the pure fluid case in Fig. 2. As more particles are applied, the Reynolds stresses become more anisotropic and display almost perfect agreement between the two models in Figs. 4 and 5. The agreement is maintained in Fig. 6 despite the differing channel velocity. The particle-phase Reynolds stresses in Figs. 3-6 show consistently good agreement between the EE-AG and EL models at all mass loadings. On the other hand, the EE-AG model for the lower mass loading cases underpredicts the uncorrelated kinetic energy of the particle phase near the wall. This may be an undesirable side effect of the near-wall blending operation used to increase the time step. Recall, however, that the EL prediction for the particle-phase kinetic energy is most accurate for the sum of the correlated and uncorrelated components because the individual components are found by a filtering operation (Capecelatro et al., 2014).  As the concentration of particles increases in Fig. 7, the production of turbulent kinetic energy shifts from classical velocity-gradient-driven production to production dominated by drag. This drag is significant only in the streamwise direction, resulting in the progression to anisotropy seen in Figs. 3-6. With this change, the dissipation of fluid-phase turbulent kinetic energy shifts in Fig. 8 from fluid-phase dissipation to interphase drag exchange to the particle phase. It is due to latter mechanism that CIT appears at higher mass loadings.
The effects of this transition can be observed directly in the EE-AG results through examining instantaneous profiles of the streamwise velocity and solid fraction in Fig. 9. In these plots, areas of relatively high particlephase concentration are superimposed onto the instantaneous fluid-phase velocity profiles. These clusters are identified for each case as cells where the particle volume fraction is more than 1.75 times the mean particle-phase volume fraction. In the low mass loading case in Fig. 9(A), even the areas with relatively large particle concentration do not have much of an effect on the turbulent behavior of the fluid phase. As more particles are applied as seen in Fig. 9(B) and (C), the fluid behavior visibly relaminarizes as more of the turbulent kinetic energy is transferred to the particle phase. At the largest mass loading for the small Stokes number particles in Fig. 9(D), the flow of the fluid phase can be seen carving itself around the areas where there are high concentrations of particles, showing the beginning of CIT behavior. Figure 10: Small Stokes particle-phase volume fraction profiles comparison.
The particle hold-up profiles in Fig. 10 also demonstrates a transition in both the EE-AG and EL models where the particles redistribute from near the wall at the lowest mass loading to the center of the channel as more particles are added. The particle-phase pressure tensor profiles in Fig. 3-6 show that both the EL and EE-AG model maintain the anisotropy of the particle-phase pressure tensor at each mass loading.  Figure 11 shows the instantaneous fluid-phase velocity profiles overlaid with the cells with 1.5 times the mean particle concentration for the large Stokes particles. The streamwise fluid velocities for the dilute cases in Fig. 11(A), (B), and (C) show no variation down the length of the channel and with only minor variation in (D). The lack of fluid-phase flow instabilities is in complete contrast to not only the small Stokes particles at equivalent mass loadings in Fig. 9, but also the pure fluid case. Additionally, particle clusters only appear at the highest mass loading case in Fig. 11(E), indicating a high degree of uniformity of the particle profile for the other cases. This is supported by the averaged wall-normal particle volume fraction profiles in Fig. 12 showing that the particle profiles are nearly completely uniform at the lowest mass loading cases. In contrast to the EE-AG model, the EL simulations remain turbulent at the lower mass loading for the large Stokes particles.  Fluid-phase turbulent kinetic energy is denoted as () and (), correlated particle-phase turbulent kinetic energy as (   Figure 15: ML= 10 case comparison of EE-AG and EL data using large Stokes number particles. Streamwise components (z or zz) for the EE-AG and EL cases are denoted respectively as () and (), wall-normal components (yy) as () and (), spanwise components (xx) as (]) and (), and shear components (yz) as () and (). Fluid-phase turbulent kinetic energy is denoted as () and (), correlated particle-phase turbulent kinetic energy as ( At these greater mass loadings shown in Figs. 13-15, the EE-AG fluid-velocity gradient matches its EL counterpart quite well. However, the EE-AG particle-phase velocities include a wall-normal dip in the nearwall region for each case that is absent from the EL results. At all mass loadings the EE-AG model fluidphase turbulence is highly anisotropic while the EL model is only the same at the greater mass loadings. The particle-phase turbulence is consistently anisotropic between both models. The particle pressures in Figs. 13-15 show consistent agreement of highly anisotropic particle pressure in both models. The behavior of the EE-AG small mass loading cases using the large Stokes particles at first appears to be a counter-intuitive result given both the small Stokes particle agreement and the interphase momentum coupling terms in fluid-phase and particle-phase momentum transport in (12) and (14), respectively. With large Stokes particles, the Stokes drag timescale is larger and allows for more decoupling between the particle and fluid phases because of the smaller drag momentum coupling term. Instead of the smaller term freeing up the fluid phase to behave more like the pure fluid case, even an extremely small volume of large Stokes particles is able to completely eliminate any significant wall-normal and spanwise flux and the corresponding turbulence in those directions.
The key element to explaining this divergence is the ability of the EL model to capture the specific local impact of each individual particle present in the system on the fluid phase that is absent in the EE-AG model. The EE-AG model instead involves redistributing the impact of those particles into a locally-averaged drag term. The difference between these two approaches is greatest where there is only a small particle number concentration. In the mass loading of 0.2 case using large Stokes particles, for example, only a tiny fraction of the 31.6 million cells contain one of the 44.8 thousand particles. Under these conditions, particle-particle collisions are rare, making the AG distribution inaccurate for representing the particle velocity distribution. When the local averaging in the EE-AG model spreads out the fluid-particle coupling too thin at lower large Stokes particle concentrations, this difference is observed in how the fluid-phase flow instabilities are dampened due to the smaller drag term. With the greater drag term present in the small Stokes particle cases at equivalent mass loadings, any given tiny flow instability that develops into the characteristic turbulent eddies and churning in the fluid phase is more readily transferred via drag to the particle phase. The particle phase ends up behaving much like the fluid phase in those cases and is itself turbulently redistributed throughout the channel. But with the larger particles, the tiny flow instabilities that initially occur in the fluid phase have less leverage to coax the particle phase to move along with it. At the smallest particle concentrations, this results in all the acceleration from the tiny instabilities getting wasted trying to move an uncooperative particle phase. This is why the distribution of particles in those cases is so uniform. As the particle concentration is increased, the fluid-phase drag term increases with it, as do particle-particle collisions. This is why beginning with the mass loading of 4 case there begins to be maldistribution in the particle holdup profile as the instabilities finally begin to propagate through to the particle phase.

Conclusions
This study modeled gas-particle flow in a vertical channel (Capecelatro et al., 2018) using both EL and EE-AG methodologies at a variety of dilute particle concentrations with particles that correspond to both large and small Stokes numbers. Primary and second-order turbulent statistics were extracted from all models and cases as a means of comparison. Using the small Stokes particles, the EE-AG model was able to successfully capture the behavior of the EL model. With the exception of the 4 mass-loading case where there was a parameter deviation in the EL model, the velocity, particle holdup, particle-pressure tensor, and Reynolds stresses extracted from both models showed reasonably good agreement. Minor differences may be attributed to the different statistics extraction methodologies between EE-AG and EL models as well as the contrasting solution methodology employed between the NGA and OpenFOAM environments. Importantly, the EE-AG model was able to capture the transition in the fluid-phase turbulence observed in the EL results from fully developed turbulent flow to a relaminarized flow regime to CIT as more and more particles were added to the channel (Capecelatro et al., 2018). From these results, the EE-AG model presents itself as a promising alternative to the EL model for dilute systems containing small Stokes particles, i.e., when the particle number concentration is large enough to warrant an Eulerian description.
The EE-AG model had far more mixed success in matching the EL results when large Stokes particles were used in combination with dilute conditions where particle-particle collision are infrequent. The low mass-loading cases for the EE-AG model in particular showed a nearly uniform particle holdup profile and a flow regime absent of significant flux in non-streamwise directions. Both of these observations are inconsistent with EL results. This difference can be attributed to the inability of the continuum drag model to describe the behavior of systems that are only sparsely populated by particles, resulting in the dampening of the formation of particle-phase instabilities. At larger mass loadings the system becomes spatially uniform enough to overcome this dampening effect and the primary flow statistics for the EE-AG model begin to agree much more with the EL results. Given that the cases where this is a problem involve dilute concentrations of large particles, this is not an issue in modeling these systems efficiently given that EL models in general excel in modeling systems under those conditions. These results indicate that the EE-AG model may still offer a good alternative to EL modeling with large Stoke particles at larger mass loadings when the number of particles begins to become too computationally costly for the EL model. A key point in future work for both small and large Stokes particles involves exploring these larger concentrations to check whether the EE-AG model is capable of accurately capturing CIT and dense particulate flows (e.g., fluidized beds), ideally with the aid of experimental data.