A stochastic model for the hydrodynamic force in Euler–Lagrange simulations of particle-laden flows

Eulerian–Lagrangian (EL) methods, in which particles are tracked individually and the fluid is solved on an Eulerian grid, have gained considerable traction for modeling particle-laden flows due to a balance between speed and resolution [1, 2, 3, 4]. In recent years, emphasis has been placed on moderate to high mass loading where particles have a first-order effect on the underlying fluid flow disturbances, which feed back to the particle dynamics. Since EL methods do not resolve the fluid boundary layer at the surface of each particle, they enable grid spacings on the order of, or larger than, the particle diameter. The reduced resolution in EL methods requires a model for the hydrodynamic fluid-particle force. By contrast, particle-resolved direct numerical simulation (PR–DNS) resolves the fluid boundary layers around each particle, and thus, the interphase drag force is an output from such simulations. For strongly-coupled flows with inertial particles, increasing the quantitative agreement between EL methods and PR–DNS requires critical assessment of the drag force model. Generally speaking, existing models for drag only capture low-order statistics, such as the mean hydrodynamic force exerted on a particle-laden suspension. It is now recognized that suspensions will exhibit significant variance in the drag force due to interactions between particles and fluid disturbances generated by their neighbors. Variance in drag force, arising from neighbor-induced fluid velocity fluctuations (pseudo-turbulent kinetic energy; PTKE), is generally ignored in EL frameworks. However, recent works have highlighted the importance of PTKE in particle-laden flows; see the closed-form model [5] and transport equations [6].

It has become increasingly well established that standard EL methods, which employ a mean drag closure and neglect PTKE-induced drag disturbances, under-predict drag variance when compared to PR–DNS [7, 8, 9]. To this end, multiple PR–DNS studies have demonstrated that flow past a collection of monodisperse spheres yields normally distributed drag forces [10, 11, 12, 13] whose standard deviation is comparable in magnitude to its mean [10, 12]. Due to the formation of fluid wakes (PTKE), particles will interact with each other indirectly over length scales comparable to their diameter. These particle-wake interactions give rise to short-range drag perturbations that drive relative motion between neighboring particles, such as drafting-kissing-tumbling (DKT) [14]. Therefore, neglecting higher-order drag statistics, resulting from neighbor-effects, can detrimentally impact EL predictions for higher-order particle statistics (velocity variance and dispersion) [9]. To-date, few drag models have been proposed that account for neighbor-induced disturbances; which may be broadly grouped into deterministic [15, 16, 17] and stochastic approaches [8, 18, 19]. In the former, the drag force experienced by a given particle is directly mapped to its pairwise neighbor interactions, requiring that the relative position of each particle be known when computing the drag force. It is worth noting that the aforementioned information is available in an EL framework but not in an Euler–Euler (EE) framework where the solids are treated as a continuum. By contrast, stochastic approaches aim to capture higher-order particle statistics, resulting from drag fluctuations, without detailed knowledge of each particle position.

Here, a statistical approach is employed to account for neighbor-induced drag fluctuations. Specifically, we follow the mathematical theory derived by Lattanzi et al. [19] for a hierarchy of Langevin equations and employ a stochastic drag force treatment. The interested reader is referred to the cited work for a more inclusive discussion on what the force Langevin (FL) framework yields in terms of particle-phase moments. However, we note that the fluctuating drag statistics obtained from an Ornstein-Uhlenbeck (OU) process were shown to be consistent with drag statistics extracted from freely-evolving PR–DNS of homogeneous systems. Namely, the fluctuating drag is normally distributed with an exponential Lagrangian autocovariance function. As a result of capturing the drag fluctuation statistics, the steady particle velocity variance obtained with FL also agreed with PR–DNS of homogeneously distributed inertial particles at moderate Reynolds numbers. In this manner, the stochastic drag force is designed to reproduce statistics obtained from fully-resolved simulations and is not an empirically added fluctuation.

We emphasize that the present study is centered around a stochastic description of neighbor-induced drag disturbances, where the physical mechanism for drag perturbations is attributed to fluid flow disturbances generated by particles that are in close proximity. It should be noted that the concept of drag force disturbances is not intrinsic to dense suspensions but will also be present in under-resolved simulations of turbulent flows. Specifically, unresolved fluid turbulence also provides a source for particle drag disturbances in dilute flows where neighbor-induced effects are less significant. Pioneering works on turbulent single-phase flows [20, 21, 22, 23] have developed Langevin equations for reconstructing the total fluid velocity that may be applied to dilute multiphase flows [24, 25, 26]. Therefore, the fluid velocity Langevin model is akin to the force Langevin framework noted above with a crucial difference being that closures for the statistical process are appropriate for dilute turbulent flows without significant two-way coupling. Drag fluctuations arise in the former case from under-resolving the intrinsic fluid turbulence; while in the latter case, they arise from under-resolving the pseudo-turbulence generated by boundary layers around neighboring particles.

The remainder of the paper is arranged as follows. In Sec. II, a statistical description is introduced for the hydrodynamic force felt by a particle in a dynamical suspension. The drag force is decomposed into a mean and fluctuating component, where the fluctuating component is a stochastic variable that specifies higher-order statistics. In Sec. III, closure is proposed for the fluctuating drag within homogeneous suspensions of inertial particles at moderate Reynolds numbers. Details regarding the numerics are provided in Sec. IV. Homogeneous fluidization of elastic particles is considered in Sec. V as a canonical flow for comparison to the PR–DNS data of Tenneti et al. [8] and Tavanashad et al. [27]. Specifically, depending upon the initial conditions, particle velocity variance in homogeneous fluidization will grow (heat) or decay (cool) to a steady state where neighbor-induced drag disturbances are balanced by hydrodynamic dissipation. Therefore, the ability of an EL method to capture the evolution of particle velocity variance in homogeneous fluidization is a significant test of the method’s ability to capture higher-order drag statistics. We first demonstrate that the new stochastic EL framework yields convergent velocity variance in homogeneous fluidization, while the velocity variance resulting from standard EL inherently depends upon the length scale employed during two-way coupling (i.e., grid spacing or filter size). Next, the proposed stochastic EL framework is directly compared to PR–DNS data in the fluidized homogeneous heating system (FHHS) and fluidized homogeneous cooling system (FHCS). Finally, we assess the role of neighbor-induced drag fluctuations in a large-scale simulation of cluster-induced turbulence (CIT) in Sec. VI.

Particle motion follows Newton’s second law where acceleration results from the net force acting upon the body. When considering the hydrodynamic drag force exerted on particle ‘$i$,’ $\bm{F}_{\rm{inter}}^{(i)}$, an exact integration of the fluid stress tensor $\check{\bm{\tau}}$ over the surface of a particle $\Gamma$ may be evaluated as

where $\check{(\cdot)}$ denotes a microscale quantity prior to any averaging, ${\rm d}A$ is an infinitesimal area element, and $\bm{n}$ is the unit normal vector outward from the particle surface. For an isolated sphere subject to non-uniform Stokes flow, Maxey and Riley [28] consider a rigorous evaluation of Eq. (1) that yields contributions from the undisturbed fluid flow, quasi-steady drag, added mass, and Basset history: $\bm{F}_{\rm{inter}}^{(i)}=\bm{F}_{\rm{un}}^{(i)}+\bm{F}_{\rm{qs}}^{(i)}+\bm{F}_{\rm{am}}^{(i)}+\bm{F}_{\rm{uv}}^{(i)}$. However, for a general dynamic suspension, finite Reynolds number effects and neighbor disturbances do not allow a first-principles solution to be obtained. Consequently, standard EL methods generally rely upon drag correlations.

Following Anderson and Jackson [29], the fluid stress tensor in Eq. (1) may be decomposed into a filtered component $\bm{\tau}$ and residual $\bm{{\tau}}^{\prime}$, such that $\check{\bm{{\tau}}}=\bm{\tau}+\bm{{\tau}}^{\prime}$. Employing the divergence theorem and choosing a filter length scale such that $\bm{\tau}$ varies little over the volume of the particle, one obtains

where $\mathcal{V}_{p}^{(i)}$ is the volume of the particle and $\bm{\tau}\left[{\bm{X}_{p}^{(i)}}\right]$ is the filtered stress tensor evaluated at the position of the particle $\bm{X}_{p}^{(i)}$. Traditionally, the unresolved drag force $\int\bm{\tau}^{\prime}\cdot\bm{n}\hskip 2.15277pt{\rm d}A$ is closed with correlations obtained from experiments [30, 31, 32] or PR–DNS [33, 34, 35, 36]. However, these correlations only capture the average hydrodynamic force acting on a suspension. Consequently, particles that experience the same filtered hydrodynamic environment will experience the same modeled drag force, even though the neighbor-induced flow may cause significant departure from the mean contribution to drag. It should be noted that the first term on the right-hand side of Eq. (2) takes the same form as $\bm{F}_{\rm{un}}^{(i)}$ in the classical formulation of Maxey and Riley [28], but inherently contains the effects of particle disturbances due to two-way coupling. The second term involving $\bm{\tau}^{\prime}$ contains all of the other contributions (quasi-steady drag, added mass, and Basset history). Rather than attempting to tease out how each pair-wise neighbor interaction contributes to the hydrodynamic forces on a given particle, we seek a statistical representation that treats these effects stochastically. To build upon the standard EL approach and account for neighbor effects, we expand the unresolved drag into a quasi-steady contribution $\bm{F}_{d}^{(i)*}$ and fluctuating $\bm{F}_{d}^{(i)\prime\prime}$ component according to

where the double-prime notation is used here to denote a fluctuation around a modeled term while the single-prime denotes a fluctuation with respect to a filtered quantity. Generally speaking, $\bm{F}_{d}^{(i)\prime\prime}$ corresponds to a perturbation from the quasi-steady drag force. Here, we attribute force perturbations to fluid flow disturbances created by neighboring particles. Therefore, $\bm{F}_{d}^{(i)\prime\prime}$ is a stochastic variable whose statistics, such as distribution and time correlation, are designed to be consistent with PR–DNS of freely-evolving particles. In this manner, $\bm{F}_{d}^{(i)\prime\prime}$ allows higher-order drag statistics, originating from neighbor effects, to be directly enforced within a fluidized suspension (see Fig. 1). Closure of the stochastic process utilized to model $\bm{F}_{d}^{(i)\prime\prime}$ is provided in Sec. III.

When applying the decomposition in Eq. (3) to an EL method, it must be noted that there are differences in how PR–DNS studies and EL methods characterize a system. Specifically, EL methods generally employ instantaneous particle velocities and locally interpolated fluid quantities when evaluating drag force correlations; whereas PR–DNS correlations are derived from ensemble-averaged quantities. To this end, we immediately define the mean Reynolds number employed by PR–DNS

and the particle Reynolds number employed by EL

where $\phi$ is the solids volume fraction, $\rho_{f}$ is the fluid density, $d_{p}$ is the particle diameter, $\mu_{f}$ is the dynamic viscosity, $\left\langle\bm{W}\right\rangle=\left\langle\bm{u}_{f}\right\rangle-\left\langle\bm{U}_{p}\right\rangle$ is the mean slip velocity, and $\bm{u}_{f}$ and $\bm{U}_{p}$ are the fluid and particle velocity, respectively. The $\left\langle\cdot\right\rangle$ notation denotes an ensemble-average while the $\left[{\bm{X}_{p}^{(i)}}\right]$ notation is suppressed on the solids volume fraction for readability. Therefore, the $\bm{F}_{d}^{(i)*}$ contribution in EL, computed with the particle Reynolds number, will be an approximation to the ensemble-averaged mean drag force exerted on a suspension $\left\langle\bm{F}_{d}\right\rangle$. More precisely, $\bm{F}_{d}^{(i)*}$ is a stochastic model for the mean hydrodynamic force, arising from one-particle statistics (position and velocity) in a filtered realization of the fluid field, that is based on a form of average drag from PR–DNS. Therefore, $\bm{F}_{d}^{(i)*}$ is inherently coupled to the momentum filtering employed by EL and will approach a delta function centered at $\left\langle\bm{F}_{d}\right\rangle$ in the limit of infinite filter width; but for finite filter length scales, $\bm{F}_{d}^{(i)*}$ will have finite variance due to variation of interpolated quantities and particle velocity (see Sec. V.1). In the present work, focus is given to filter widths larger than a particle diameter $\delta_{f}=\mathcal{O}(10d_{p})$ that do not contribute to significant variation of $\bm{F}_{d}^{(i)*}$ in homogeneous fluidization. To avoid confusion when discussing mean drag, we draw an analogy with the Maxey-Riley equation and refer to $\bm{F}_{d}^{(i)*}$ as the quasi-steady drag force. We further motivate this analogy by noting that the present work focuses on inertial particle suspensions where added mass and Basset history effects are less significant and PR–DNS correlations are obtained from fixed particle simulations.

The force Langevin (FL) theory examined by Lattanzi et al. [19] is employed here to describe the neighbor-induced fluctuating drag force $\bm{F}_{d}^{(i)\prime\prime}$. Specifically, $\bm{F}_{d}^{(i)\prime\prime}$ follows an Ornstein-Uhlenbeck (OU) process according to

where $\tau_{F}$ is the integral time scale of the fluctuating drag force, $\sigma_{F}$ is the standard deviation of the fluctuating drag force, and $\mathrm{d}\bm{W}_{t}$ is a Wiener process increment. We refer the interested reader to Lattanzi et al. [19] where a detailed discussion regarding motivation for, and results from, the OU process are provided. Here, we briefly emphasize that the steady solution to the OU process is a normal distribution $\mathcal{N}\left[0,\hskip 2.15277pt\sigma_{F}\right]$ and a multitude of PR–DNS studies have reported normally distributed drag forces $\mathcal{N}\left[\left\langle F_{d}\right\rangle,\hskip 2.15277pt\sigma_{F}\right]$ [7, 10, 12, 18, 13]. Additionally, FL was shown to be reconcilable with PR–DNS results for granular temperature evolution since it correctly attributes the source of granular temperature to the force-velocity covariance [19].

In general, coefficients of the OU process in Eq. (6), namely the drag time scale $\tau_{F}$ and standard deviation $\sigma_{F}$, may be tensorial quantities that correlate the drag fluctuations in each direction and drive anisotropic granular temperature development. While a valuable and interesting problem, we first consider the isotropic case given in Eq. (6) so that the same force time scale and standard deviation is applied to all three fluctuating force directions. We further note that particle collisions will relax the granular temperature back towards isotropy. For statistically homogeneous suspensions, Garzó et al. [37] derived an evolution equation for the particle velocity anisotropy tensor that shows a return to isotropy occurs at steady state due to particle collisions. For dilute particle flows with small density ratio $\rho_{p}/\rho_{f}=\mathcal{O}(1)$, where viscous lubrication is significant and particle collisions are less frequent, anisotropy may play a significant role. However, a description for such flows will need to be addressed in future work as the focus here is on inertial gas-solid suspensions.

In this section, the stochastic EL framework and its discretization are summarized, using closures reported in Sec. III to model the drag force perturbations.

We consider the homogenenous fluidization of particles within a triply periodic cube of length $L$. A mean fluid flow rate $\left\langle\bm{u}_{f}\right\rangle=\left[0\hskip 4.30554pt0\hskip 4.30554ptu_{f,z}\right]^{\top}$ is imposed via $\mathcal{F}_{\rm mfr}$ in Eq. (18b) to obtain a desired mean Reynolds number ${\rm Re}_{m}$, while the gravitational body force $\bm{g}=\left[0\hskip 4.30554pt0\hskip 4.30554pt-g_{z}\right]^{\top}$ is prescribed in the opposite direction. $\mathcal{F}_{\rm mfr}$ is added as a uniform source term and is computed each timestep to enforce the constant $u_{f,z}$. The body force $g_{z}$ is set equal to the mean hydrodynamic acceleration $\left\langle\bm{F}_{d}\right\rangle/m_{p}$, computed via Eq. (15) for the mean conditions ${\rm Re}_{m}$ and $\left\langle\phi\right\rangle$. For the homogeneous cases considered here, the drag force offsets the weight of the suspension, leading to a mean particle velocity that is approximately zero throughout the simulation. In contrast, the particle velocity variance will grow or decay to a steady value that is dictated by the balance of drag variation and hydrodynamic dissipation. The simulation domain closely reflects the PR–DNS studies of Tenneti et al. [8] and Tavanashad et al. [27], which serve as benchmark data for comparison to the stochastic EL method presented here. Unless otherwise noted, the grid spacing $\Delta x/d_{p}=0.5$, kernel width $\delta_{f}/d_{p}=7$, and domain size $L/d_{p}=7$ were held fixed for the homogeneous fluidization simulations, and a summary of relevant conditions is provided in Table 1.

Within the homogeneous fluidization system, we define two canonical flows that result from different initial conditions for the particle velocity. Namely, we first examine the fluidized homogeneous heating system (FHHS) in Sec. V.2, where particles are initialized with zero velocity and particle velocity variance grows to a steady value. We then examine the fluidized homogeneous cooling system (FHCS) in Sec. V.3, where particles are initialized with an over-prescribed velocity variance that decays to a steady value. To match the FHCS simulations completed by Tenneti et al. [8], we sample the initial particle velocities from a Maxwellian distribution.

It is important to note that domain sizes $L$ considered by the PR–DNS studies of Tenneti et al. [8] and Tavanashad et al. [27] are sufficiently small that particle clustering is not observed, and the system remains homogeneous. It has been well established that large-scale fluidized systems exhibit the classic clustering instability due to two-way momentum coupling and/or dissipative collisions. We first focus on the homogeneous case and define crucial parameters to facilitate comparison between EL and PR–DNS benchmark data. Specifically, particle velocity variance $T$ is computed with a velocity fluctuation

that is defined with respect to the domain average $\left\langle\bm{U}_{p}^{(i)}\right\rangle$, leading to

Utilizing the definition for particle velocity fluctuation in Eq. (25) and particle velocity variance in Eq. (26), we define the particle Reynolds stress tensor $\bm{R}^{p}$, anisotropy tensor $\bm{b}^{p}$, and thermal Reynolds number ${\rm Re}_{T}$ as

and

We emphasize that the mean-free time $\tau_{\rm col}$ in Eq. (8) is always computed with the granular temperature $\Theta$ via Eq. (24), making it valid for inhomogeneous flows. However, when comparing to PR–DNS of homogeneous systems, we report Eqs. (27a)–(28) with the particle velocity variance $T$, computed via Eq. (26), for consistency with the definitions employed in those studies. In Sec. VI we finally consider large-scale simulations of cluster-induced turbulence (CIT) that are characterized by strong inhomogeneities in particle number density and clusters that fall faster than their terminal velocity.

Until this point, the domain size of the systems under consideration were small enough such that instabilities giving rise to heterogeneity in particle concentration were suppressed. However, large-scale flows with inertial particles will lead to the spontaneous formation of clusters due to two-way momentum coupling and/or dissipative collisions [50, 51, 47]. As discussed by Fullmer and Hrenya [51], neighbor-induced drag fluctuations act as a source of granular temperature that may contribute to cluster break up. To probe the role of neighbor-effects on cluster statistics, we consider simulations of fully-developed cluster-induced turbulence (CIT) with standard EL and stochastic EL. The CIT simulation geometry is essentially identical to the homogeneous fluidization described in Sec. V, with the exception that the domain length is increased to resolve the cluster length scale [50]

Following Capecelatro et al. [52], we set the streamwise domain length as $L_{z}=32\mathcal{L}$ to obtain converged statistics and avoid clusters falling in their own wakes. The simulation conditions are summarized in Table 3. To quantify the degree of particle segregation in homogeneous isotropic turbulence, Eaton and Fessler [53] proposed a clustering parameter $D$

where $\sigma_{p}$ is the standard deviation in solids volume fraction for a random distribution of particles (evaluated at the beginning of the simulation for the initial random particle configuration). Here, we utilize $D$ as an indicator of the degree of clustering within the entire domain, so as to probe the effect of neighbor-induced drag fluctuations on cluster break up. While $D$ provides an estimate for the degree of clustering, we note that more detailed descriptions of heterogeneous media have been proposed that depend upon the viewing window size [54, 55]. While these methods are beyond the scope of the present work, they do provide a more inclusive picture of the clustering spectrum, as opposed to a single value obtained with $D$.

Due to the presence of clusters, CIT exhibits drag reduction when compared to homogeneous systems. Specifically, entrainment of the surrounding fluid by particle clusters leads to a mean particle velocity in the streamwise direction $\left\langle U_{p,z}\right\rangle$ that is non-zero and in the direction of gravity (see FIG 8 (a)). We note that simulations are completed here in a reference frame that moves with the mean slip velocity $W_{z}$, which is specified by the mean Reynolds number ${\rm Re}_{m}$ and the assumption that $\left\langle U_{p,z}\right\rangle\approx 0$; see discussion in Sec. V. Therefore, drag reduction in CIT leads to particle acceleration in the direction of gravity and mean settling velocities in the present simulations that are $\sim 2.25W_{z}$. Examining the evolution of mean particle settling velocity over the course of the simulations shows there are negligible differences between stochastic EL and standard EL (see Fig. 8 (a)). Similarly, negligible differences are observed between stochastic EL and standard EL for the clustering parameter (see Fig. 8 (c)). Since the degree of clustering is tied to drag reduction, and by extension the mean settling velocity, $D$ and $\left\langle U_{p,z}\right\rangle$ are strongly correlated. For example, a significant reduction in $D$ would be reflected in an increase in $\left\langle U_{p,z}\right\rangle$, due to drag enhancement from homogenization of the particle phase. Since $D$ and $\left\langle U_{p,z}\right\rangle$ show consistent behavior, no appreciable change with the stochastic EL method, the model for neighbor-induced drag does not play a significant role in these domain averaged quantities, though it may alter the spectra obtained from the method of Quintanilla and Torquato [55] .

Similar to mean settling velocity and clustering parameter, the particle velocity variance also shows minor deviation between stochastic and standard EL during transients and at steady state; see in Fig. 8 (b). However, the beginning of the simulation $t/\tau_{p}\in\left[0\hskip 2.15277pt5\right]$ does show significant deviation between the velocity variance developed with stochastic EL and standard EL (see inset of Fig. 8 (b)). At this point, particles are randomly, but still homogeneously, distributed within the domain. Under these conditions, the neighbor-induced drag fluctuation is the predominate source of velocity variance since clusters with sharp solids volume fraction interfaces have yet to form. As the system enters the transient stage $t/\tau_{p}\in\left[5\hskip 2.15277pt30\right]$, meso-scale particle structuring begins and the solids accelerate in the streamwise direction. The formation of clusters leads to significant quasi-steady drag variance at the edge of clusters, where large gradients in solids volume fraction and fluid velocity occur. These quasi-steady drag variances are resolved by a standard EL method and will facilitate granular temperature transport through the cluster via shear generation and collisional conduction. Therefore, heterogeneous systems provide additional sources to particle velocity variance that can dominate over the neighbor-effect. To illustrate this point, we extract probability distributions for quasi-steady and fluctuating drag in CIT and homogeneous fluidization (see Fig. 9). For the case of homogeneous fluidization, fluctuating drag is more significant than quasi-steady drag; while in the case of CIT, the exact opposite holds. In addition, we reiterate that the stochastic force is isotropic and was shown in Sec. V.1 to yield isotropic particle velocity fluctuations in homogeneous fluidization. However, for the case of CIT, stochastic EL and standard EL yield the same anisotropy $b^{p}_{\parallel}=0.35$ and $b^{p}_{\perp}=-0.16$, which is further evidence that resolved quasi-steady drag is dominant over neighbor-induced drag.

In the present study, we examine the role of higher-order drag statistics, originating from neighbor-induced fluid flow perturbations, in Eulerian–Lagrangian (EL) methods. A model was proposed for neighbor-induced hydrodynamic forces by treating the fluctuating drag as a stochastic variable that follows an Ornstein-Uhlenbeck process. Specifically, the force Langevin (FL) method detailed in Lattanzi et al. [19] was utilized here to construct a stochastic EL framework. Closures are provided for the theoretical inputs to the FL equation, integral time scale $\tau_{F}$ and force standard deviation $\sigma_{F}$, that are appropriate for inertial particles at moderate Reynolds numbers. Specifically, the integral time scale of the fluctuating drag is approximated with the mean-free time between successive collisions, derived from the kinetic theory of non-uniform gases. The standard deviation in drag force is closed with a new correlation based upon particle-resolved direct numerical simulation (PR–DNS) of fixed assemblies. The new stochastic EL framework specifies unresolved drag statistics through the stochastic force, leading to the correct evolution and sustainment of particle velocity variance when compared to PR–DNS of freely-evolving homogeneous suspensions. Since standard EL infers drag statistics from variations in the resolved flow, it cannot replicate the higher-order particle statistics (velocity variance and dispersion) observed in PR–DNS of homogeneous suspensions. Finally, the role of neighbor-induced drag fluctuations on cluster-induced-turbulence (CIT) is considered. In contrast to homogeneous fluidization, CIT is characterized by large gradients in solids volume fraction and fluid velocity. The aforementioned gradients provide a source for drag variance in standard EL, through the quasi-steady drag closure, that can dominate over neighbor effects. Since standard EL resolves the dominate modes for generating granular temperature in heterogeneous systems — quasi-steady drag variation at cluster interface and collisional conduction — we observe negligible change when employing the stochastic EL framework.

While emphasis is placed here on inertial particle suspensions at moderate solids loading and Reynolds numbers, we stress that the proposed methodology is a general framework that may be adapted to a variety of applications. Namely, the concept of higher-order drag statistics can be readily applied to higher-order statistics in Nusselt and Sherwood correlations, employed for simulation of heat and/or mass transfer. Additionally, compressible particle-laden flows often exhibit significant drag variation and unsteady effects (added mass and Basset history). The present framework provides a stepping stone that may be leveraged by future work to account for such effects.

The authors report no conflict of interest.

This material is based upon work supported by the National Science Foundation under grant no. CBET-1904742 and grant no. CBET-1438143. The authors acknowledge the Texas Advanced Computing Center (TACC) for providing Stampede2 compute resources under Extreme Science and Engineering Discovery Environment (XSEDE) grant no. TG-CTS200008.