2. System under consideration

Three types of Langevin equations are discussed in the following sections and the evolution equations for velocity variance and dispersion are derived. The present study considers dispersion from a line source in a statistically homogeneous flow (mean particle velocity is homogeneous). We motivate the case study by considering particles in an EL method, where unresolved fluid velocity fluctuations (PTKE) give rise to drag force fluctuations that induce velocity variance and particle dispersion. In all cases, we consider an initial impulse of particles centred at $x=0$ and uniformly distributed in the y-direction. The particle velocities are initialized as either (i) all initially static i.e. $U_{p,i}=0; \ \delta (u)$ or (ii) a Maxwellian velocity distribution $\mathcal {N}[0,\sigma ^2]$; where $\mathcal {N}$ is a normal distribution with zero mean and variance $\sigma ^2$.

Analytical solutions derived here are compared to results obtained from numerical solutions of the Langevin equations. In the present study, numerical integration of the Langevin equations is considerably more straightforward than deriving the exact solution. For this reason, simulation results are considered less error prone and included here as a confirmation of the analytical results. The SDEs are integrated with an Euler–Maruyama method (Kloeden & Platen Reference Kloeden and Platen1992) with $10^{5}$ particles, and the smallest time scale of the system is resolved by 25 time steps. We do not consider particle collisions in the present work, and thus, the two relevant time scales are those associated with mean drag defined as $\tau _p = \rho _p d_p^2/(18 \mathrm {\mu })$ and the fluctuating drag force $\tau _F$; where $\rho _p$ is the particle density, $d_p$ is the particle diameter and $\mathrm {\mu }$ is the fluid dynamic viscosity. The computational domain is a doubly periodic rectangular grid with aspect ratio of 6. Probability distributions for the particle phase were obtained by decomposing the length of the domain into 501 bins and computing the particle number density within each bin. It is worth noting that particle dispersion here is one-dimensional in nature and particles do not interact with the periodic $x$-direction boundaries. The two-dimensional domain merely allows for more intuitive visualization; see figure 1 for example.

Figure 1.

(a,b) Analytical and numerical predictions for particle dispersion with the VL framework and an impulse velocity condition. (c,d) Analytical and numerical predictions for particle dispersion with the VL framework and a Maxwellian velocity condition. In all cases, $St = 5$, $Pe = 2$ and the particles (left half of images) are coloured by the magnitude of their velocity in the $x$-direction. $(a{,}c)\,\, t = 0.500, (b{,}d)\,\, t = 1.000.$

As noted above, the theory here treats particles as collisionless and emphasis is placed on hydrodynamic force fluctuations. For inelastic granular flows, a Chapman–Enskog expansion shows that the granular temperature follows a hydrodynamic equation where acceleration is balanced by generation due to shear work, diffusion and dissipation from inelastic collisions (Garzó & Dufty Reference Garzó and Dufty1999). Extending the kinetic theory to account for fluid phase effects shows that additional source and sink terms arise from hydrodynamic forces (Koch & Sangani Reference Koch and Sangani1999; Garzó et al. Reference Garzó, Tenneti, Subramaniam and Hrenya2012). For a homogeneous flow of elastic particles, the granular temperature balance only retains the hydrodynamic source and sink terms (Tenneti et al. Reference Tenneti, Mehrabadi and Subramaniam2016). It is this fluctuation dissipation relation that we probe via stochastic methods. In regard to dispersion, collisions will affect particle migration and the closures derived here should be viewed as upper bounds that are more appropriate for dilute flows. While not considered here, the method outlined by Pai & Subramaniam (Reference Pai and Subramaniam2009) would provide a path for including the effect of a uniform shear flow.

3. Position Langevin

Evolution of a particle's position follows from the definition of its velocity

(3.1)\begin{equation} \frac{\textrm{d} \boldsymbol{X}_{p,i}}{\textrm{d} t} = \boldsymbol{U}_{p,i}, \end{equation}

where $\boldsymbol {X}_{p,i}$ is the position of particle $i$ and $\boldsymbol {U}_{p,i}$ is the velocity of particle $i$. Introducing random fluctuations into (3.1) leads to the following SDE

(3.2)\begin{equation} \mathrm{d}\boldsymbol{X}_{p,i} = \boldsymbol{U}_{p,i}\, \mathrm{d}t + {\bar{\boldsymbol{\mathsf{b}}}}(\boldsymbol{X},t)\, \mathrm{d}\boldsymbol{W}_t, \end{equation}

where ${\bar {\boldsymbol{\mathsf{b}}}}(\boldsymbol {X},t)$ is the diffusion tensor and $\mathrm {d}\boldsymbol {W}_t$ is a Wiener process increment. Utilizing Itô calculus, (3.2) corresponds to the following Fokker–Planck (FP) equation (detailed derivations of the Fokker–Planck equation are not provided here but the interested reader is referred to chapter 4 of Risken & Frank (Reference Risken and Frank1996).)

(3.3)\begin{equation} \frac{\partial P(\boldsymbol{u},\boldsymbol{x};t|\boldsymbol{v},\boldsymbol{y},s)}{\partial t} + \sum_{i=1}^{N}\frac{\partial u_i P}{\partial x_i} = \sum_{i=1}^{N} \sum_{j=1}^{M}\frac{\partial D_{ij}P}{\partial x_i x_j}, \end{equation}

where $P(\boldsymbol {u},\boldsymbol {x};t|\boldsymbol {v},\boldsymbol {y},s)$ is the probability distribution conditioned on the initial condition for velocity $\boldsymbol {v}$, position $\boldsymbol {y}$ and time $s$, $u_i$ is the particle velocity and ${\bar {\boldsymbol{\mathsf{D}}}}={\bar {\boldsymbol{\mathsf{b}}}}{\bar {\boldsymbol{\mathsf{b}}}}^{\boldsymbol{\mathsf{T}}}/{2}$. For a constant diffusivity tensor ${\bar {\boldsymbol{\mathsf{D}}}}$, (3.3) follows an advection–diffusion equation.

Here, we consider the particle diameter $d_p$, Stokes response time $\tau _p$ and terminal velocity $V_T=\tau _p g(1- \rho _f/\rho _p)$ as reference scales, where $\rho _f$ is the fluid density. From these scales, we define the Stokes number $St = V_T \tau _p/d_p$ and Péclet number $Pe = V_T d_p/D$. Taking $\bar {b}_{ij} = \sqrt {2D}\delta _{ij}$, we obtain a dimensionless form for (3.2)

(3.4)\begin{equation} \mathrm{d}\boldsymbol{X}_{p,i} = St\, \boldsymbol{U}_{p,i} \,\mathrm{d}t +\sqrt{2 St\,Pe^{-1}} \,\mathrm{d}\boldsymbol{W}_t \end{equation}

that leads to the following Fokker–Planck equation

(3.5)\begin{equation} \frac{\partial P}{\partial t} + St\, \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} P) = St\,Pe^{-1} \Delta(P). \end{equation}

In general, the Fokker–Planck equation poses a significant challenge to analytical methods. However, closed-form solutions are highly instructive when possible. For one dimension and constant velocity, a change of reference frame ($\eta = x - St\, u t$) yields the diffusion equation

(3.6)\begin{equation} \frac{\partial P}{\partial t} = St\,Pe^{-1}\frac{\partial^{2} P}{\partial \eta^{2}}. \end{equation}

An initial impulse at the origin, $\delta (x)$, leads to the classic heat kernel solution

(3.7)\begin{equation} P(\eta;t) = \frac{1}{\sqrt{4{\rm \pi} St\,Pe^{-1} t}}\exp\left(-\frac{\eta^2}{4St\,Pe^{-1}t}\right), \end{equation}

or $\mathcal {N}[0, 2 St\,Pe^{-1}t]$ for shorthand. Equation (3.6) will serve as the starting point in the following subsections and connections will be made to (3.7).

3.2. Dispersion (PL)

The variance in position given by $\mathrm {Var}(x) = \langle (x - \langle x \rangle )^2 \rangle$ is utilized here to characterize particle dispersion, for which the first two moments $\langle x \rangle; \ \langle x^2 \rangle$ are of great significance. The $\langle \cdot \rangle$ notation employed here indicates an ensemble average over all the particles and at each instance in time, although different definitions for the mean may be adopted depending upon the solution method. Higher dimension distributions will arise in subsequent sections, leading to conditional averages – e.g. $\langle x \langle u|x\rangle \rangle$ reads as the ensemble average of the product of position and velocity conditioned on position. When computing moments of the Fokker–Planck equation, the following relations are useful

(3.8a)\begin{gather} \int x^n \frac{\partial P}{\partial t} \,\textrm{d}x = \frac{\textrm{d} \langle x^n \rangle}{\textrm{d} t}, \end{gather}

(3.8b)\begin{gather} \int x^n \frac{\partial P}{\partial x} \,\textrm{d}x \equiv \int \left[\frac{\partial x^nP}{\partial x} - nx^{n-1}P\right] \,\textrm{d}x = -n\langle x^{n-1} \rangle, \end{gather}

(3.8c)\begin{gather} \int x^n \frac{\partial^{2} P}{\partial x^{2}} \,\textrm{d}x \equiv \int \left[\frac{\partial }{\partial x}\left(x^n \frac{\partial P}{\partial x}\right) - n\frac{\partial x^{n-1}P}{\partial x} + n(n-1)x^{n-2}P\right] \,\textrm{d}x = n(n-1)\langle x^{n-2} \rangle, \end{gather}

where integration is performed over the real line and the last equalities in (3.8b) and (3.8c) are obtained with integration by parts and decay of the distribution function ($\lim _{x\rightarrow \pm \infty } x^n P = 0$). (Similar approaches are employed when taking moments of the Boltzmann equation; see chapter 9 of Vincenti & Kruger (Reference Vincenti and Kruger1975). Additionally, the Gaussian functions obtained here can be shown to decay faster than a polynomial of order $n$.) Utilizing (3.8a)–(3.8c), the first and second moments of (3.6) follow

(3.9a)\begin{gather} \frac{\textrm{d} \langle x \rangle}{\textrm{d} t} = 0, \end{gather}

(3.9b)\begin{gather} \frac{\textrm{d} \langle x^2 \rangle}{\textrm{d} t} = 2St\,Pe^{-1}. \end{gather}

Employing the impulse initial condition when integrating (3.9a)–(3.9b) allows for evaluation of the position variance

(3.10)\begin{equation} \langle \left(x - \langle x \rangle \right)^2 \rangle =2St\,Pe^{-1}t. \end{equation}

Equation (3.10) shows that the dispersion obtained with a PL treatment is consistent with the known solution in (3.7) – i.e. random increments in particle position yield a diffusion process. Equation (3.10) connects the mean-squared displacement to diffusion in one dimension but will be multiplied by the dimensionality $d$ for isotropic diffusion in higher dimensions.

Returning (3.10) to dimensional form may be completed in a straightforward manner. We first note that $St\,Pe^{-1} = D \tau _p/d_p^2$. Multiplying the right-hand side of (3.10) by $d_p^2$ and replacing $t$ with $t/\tau _p$, we obtain the classic result: $\mathrm {Var}(x) = 2Dt$.

4. Velocity Langevin

Since collisions are not considered here, the relevant forces acting on the particle are due to interphase drag, $\boldsymbol {F}_{{d},i}$. In general, the drag force is proportional to the slip velocity $\boldsymbol {F}_{{d},i} = f(\phi ,Re_m) (\boldsymbol {u}_f - \boldsymbol {U}_{p,i})$, where $f(\phi ,Re_m)$ is a function of the solids volume fraction $\phi$ and particle Reynolds number $Re_m=((1-\phi ) d_p \rho _f |u_f - U_{p,i}|)/\mathrm {\mu }$. Established correlations for the drag force at varying Reynolds number and solids loading may be found in Gidaspow (Reference Gidaspow1994), Hill, Koch & Ladd (Reference Hill, Koch and Ladd2001), Beetstra, van der Hoef & Kuipers (Reference Beetstra, van der Hoef and Kuipers2007), Tenneti, Garg & Subramaniam (Reference Tenneti, Garg and Subramaniam2011) and Rubinstein, Derksen & Sundaresan (Reference Rubinstein, Derksen and Sundaresan2016). For Stokes flow, the system of SDEs for position and velocity is

(4.1a)\begin{gather} \mathrm{d}\boldsymbol{X}_{p,i} = \boldsymbol{U}_{p,i}\,\mathrm{d}t, \end{gather}

(4.1b)\begin{gather} \mathrm{d}\boldsymbol{U}_{p,i} = \frac{\boldsymbol{u}_f - \boldsymbol{U}_{p,i}}{\tau_p} \mathrm{d}t + {\bar{\boldsymbol{\mathsf{b}}}}(\boldsymbol{V},\boldsymbol{X};t) \mathrm{d}\boldsymbol{W}_t.\end{gather}

In the absence of collisions, the velocity SDE in (4.1b) contains a restoring force that relaxes the particle velocity towards the fluid velocity $\boldsymbol {F}_{{d},i}/m_{p,i} = (\boldsymbol {u}_f-\boldsymbol {U}_{p,i})/\tau _p$ and a random fluctuation ${\bar {\boldsymbol{\mathsf{b}}}}\,\mathrm {d}\boldsymbol {W}_t$.

When defining the driving force for drag, we utilize the instantaneous particle velocity $\boldsymbol {U}_{p,i}$, which is resolved in an EL method. When considering a VL approach, Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016) utilized the mean particle velocity $\left \langle \boldsymbol {U}_{p}\right \rangle$ to define the driving force for mean drag. By doing so, Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016) was able to introduce a different time scale for the relaxation of particle velocity fluctuations $\gamma$. This new time scale implies that the viscous dissipation of particle velocity fluctuations is fundamentally different than the drag force for mean particle velocity. Similarly, the dissipation of granular temperature due to viscous effects has been probed via multi-pole (Sangani & Mo Reference Sangani and Mo1996) and lattice Boltzmann simulations (Wylie, Koch & Ladd Reference Wylie, Koch and Ladd2003). These works also suggest that the viscous dissipation of particle velocity fluctuations can depart from mean drag and will scale with solids volume fraction and Reynolds number $Re_T=(\rho _f d_p \sqrt {T})/\mathrm {\mu }$ based on granular temperature $T$ (see $R_{{diss}}(\phi )$ in Sangani et al. Reference Sangani, Mo, Tsao and Koch1996; Koch & Sangani Reference Koch and Sangani1999; Wylie et al. Reference Wylie, Koch and Ladd2003). Here, we decompose the instantaneous particle velocity into a mean and fluctuating component $\boldsymbol {U}_{p,i} = \left \langle \boldsymbol {U}_{p}\right \rangle + \boldsymbol {U}_{p,i}^{\prime }$ and substitute into the drag force correlation. The process employed here is consistent with treatments employed by EL methods but requires that the velocity fluctuations inherit the mean drag time scale. The dissipation time scale obtained via direct substitution $1/\tau _p$ will be compared to results with the time scale of Sangani et al. (Reference Sangani, Mo, Tsao and Koch1996) and Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016) in § 6; $R_{{diss}}(\phi )/\tau _p$ and $\gamma$ respectively.

For a reference frame moving with the average particle velocity (similar to § 3), the drag force term in (4.1b) becomes $-\boldsymbol {U}_{p,i}/\tau _p\, \mathrm {d}t$, where we drop the $(\cdot )^{\prime }$ notation for the sake of readability and it is understood that the variables are fluctuating quantities. Taking the diffusion tensor to be $\bar {b}_{ij} = \sqrt {2D/\tau _p^2} \delta _{ij}$, the dimensionless system of SDEs is given by

(4.2a)\begin{gather} \mathrm{d}\boldsymbol{X}_{p,i} = St \boldsymbol{U}_{p,i}\,\mathrm{d}t, \end{gather}

(4.2b)\begin{gather} \mathrm{d}\boldsymbol{U}_{p,i} = -\boldsymbol{U}_{p,i}\, \mathrm{d}t + \sqrt{2\left(Pe\,St \right)^{-1}}\, \mathrm{d}\boldsymbol{W}_t. \end{gather}

The corresponding Fokker–Planck equation is

(4.3)\begin{equation} \frac{\partial P(\boldsymbol{u},\boldsymbol{x};t|\boldsymbol{v},\boldsymbol{y},s)}{\partial t} + St \boldsymbol{\nabla}_{x} \boldsymbol{\cdot} \left(\boldsymbol{u}P\right) - \boldsymbol{\nabla}_{u} \boldsymbol{\cdot} (\boldsymbol{u}P) =(Pe\,St)^{-1}\Delta_{u}(P), \end{equation}

where the subscripts on the differential operators denote the coordinates of phase space that the derivatives are taken with respect to (position or velocity). Integration of (4.3) over the spatial coordinate yields the Fokker–Planck equation for the marginal velocity distribution. In one dimension, the equation is given by

(4.4)\begin{equation} \frac{\partial P(u;t|v,s)}{\partial t} - \frac{\partial \left(uP\right)}{\partial u} = \left(Pe\,St \right)^{-1}\frac{\partial^{2} P}{\partial u^{2}}. \end{equation}

Given an impulse initial condition for the particle velocity $\delta (u-V_0)$, the solution to (4.4) may be obtained via Fourier transform and method of characteristics (Risken & Frank Reference Risken and Frank1996; Pope Reference Pope2000; Gardiner Reference Gardiner2009)

(4.5)\begin{equation} P(u;t) = \frac{1}{\sqrt{2{\rm \pi} (Pe\,St)^{-1} [1 - \exp(-2t)]}}\exp\left(-\frac{( u - V_0\exp(-t) )^2 }{2 (Pe\,St )^{-1}[1 - \exp(-2t)]}\right), \end{equation}

or $\mathcal {N}[V_0\exp (-t), (Pe\,St )^{-1} [1 - \exp (-2t)]]$ for shorthand. We note that the restoring force in the Ornstein–Uhlenbeck process yields a bounded Gaussian distribution for velocity in the long time limit ($t \gg \tau _p$) with variance of $(Pe\,St )^{-1}$; as opposed to the unbounded variance for position obtained for the diffusion process in (3.10).

4.1. Fluctuating velocity (VL)

Applying the moment relations in § 3.2 to the Fokker–Planck equation given in (4.4) leads to

(4.6a)\begin{gather} \frac{\textrm{d} \langle u \rangle}{\textrm{d} t} = - \langle u \rangle, \end{gather}

(4.6b)\begin{gather} \frac{\textrm{d} \langle u^2 \rangle}{\textrm{d} t} = -2 \langle u^2 \rangle + 2\left(Pe\,St \right)^{-1}. \end{gather}

Equation (4.6b) is the one-dimensional analogue to the granular temperature balance obtained by Koch & Sangani (Reference Koch and Sangani1999) and Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016). The first term on the right-hand side is the velocity dissipation, defined as $\varGamma _{{vis}}=2R_{{diss}}(\phi )T/\tau _p$ by Koch & Sangani (Reference Koch and Sangani1999) and directly modelled from PR-DNS data by Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016). Due to the non-dimensionalization here, the $1/\tau _p$ cancels on our dissipation term, and thus, the difference between our dissipation and that of Koch & Sangani (Reference Koch and Sangani1999) can be attributed to $R_{{diss}}(\phi )$. The second term is the source due to drag fluctuations, which VL models as a Wiener process increment. With a $V_0$ initial condition, integration of (4.6a) and (4.6b) yields

(4.7a)\begin{gather} \langle u \rangle = V_0 \exp(-t),\end{gather}

(4.7b)\begin{gather} \langle u^2 \rangle = (Pe\,St)^{-1} + (V_0^2 - (Pe\,St)^{-1} )\exp(-2t). \end{gather}

The time evolution of the velocity variance is then given by

(4.8)\begin{equation} \langle (u - \langle u \rangle )^2 \rangle = (Pe\,St)^{-1}[1 - \exp(-2t)], \end{equation}

which is again consistent with the known solution given in (4.5).

Following the arguments given at the end of § 3.2, (4.5) may be cast in dimensional form by replacing $(Pe\,St )^{-1}$ with $D/(V_T^2\tau _p)$, multiplying by $V_T^2$ and setting $t = t/\tau _p$. We note that the diffusion tensor for an Ornstein–Uhlenbeck process is often set as $\sqrt {2\sigma _v^2/\tau _p}$, rather than $\sqrt {2D/\tau _p^2}$ as is done here; where $\sigma _v^2$ is a specified velocity variance. We set the diffusion tensor such that a canonical form is obtained for the dispersion solutions, thereby allowing straightforward comparison across different frameworks (PL, VL and FL). If one wishes to employ a diffusion tensor of $\sqrt {2\sigma _v^2/\tau _p}$ for the VL framework, then the variance solutions (4.8) and (4.12a) may be converted to dimensional form by substituting $\sigma _v^2\tau _p$ in place of $D$ – i.e. in dimensional variables, the steady granular temperature in (4.8) would be directly set by $\sigma _v^2$.

4.2. Dispersion (VL)

The spatial marginal distribution (integration of (4.3) over velocity space) leads to

(4.9)\begin{equation} \frac{\partial P(\boldsymbol{x};t|\boldsymbol{y},s)}{\partial t} + St \boldsymbol{\nabla}_{x} \boldsymbol{\cdot} (\langle u | x \rangle P) = 0, \end{equation}

where $\langle u|x \rangle$ is the average velocity conditioned on position. Taking the first and second spatial moments of (4.9) leads to

(4.10a)\begin{gather} \frac{\textrm{d} \langle x \rangle}{\textrm{d} t} = St \langle u \rangle, \end{gather}

(4.10b)\begin{gather} \frac{\textrm{d} \langle x^2 \rangle}{\textrm{d} t} = 2 St \langle x \langle u | x \rangle \rangle, \end{gather}

(4.10c)\begin{gather} \frac{\textrm{d} \langle (x - \langle x \rangle)^2 \rangle}{\textrm{d} t} \equiv 2 St (\langle x \langle u | x \rangle \rangle - \langle x \rangle \langle u \rangle) = 2 St \int_0^{s_2} K_{uu}(s_1,s_2) \, \mathrm{d}s_1, \end{gather}

where $K_{uu}=Pe^{-1} \{\exp (-(s_2 - s_1)) - \exp (-(s_2 + s_1)) \}$ is the velocity autocovariance function. Obtaining the second equality in (4.10c) is non-trivial and the detailed derivation is given in appendix A.1. Integrating (4.10c), we arrive at the evolution equation for position variance.

An important simplification to the present derivation may be observed for the case of a constant velocity variance – i.e. if the particle velocity is initialized with a Maxwellian distribution $\mathcal {N}[0, (Pe\,St )^{-1}]$ ((4.5) at steady state). For the Maxwellian initial condition, the velocity variance is independent of time and the autocovariance may be written in terms of a single time lag $K_{uu}(s_1) = (Pe\,St )^{-1} \exp (-s_1)$; leading to different dispersion behaviour at early times (see derivation in appendix A.2 and discussion in § 4.3). Evaluating the integral for mean position (4.10a) leads to

(4.11)\begin{equation} \langle x \rangle = V_0 St [1 - \exp(-t)], \end{equation}

while the integral for position variance (A 5) gives

(4.12a)\begin{gather} \langle (x - \langle x \rangle)^2 \rangle = 2 St\,Pe^{-1} \left\{t - 2\mathbb{E}_1 + \frac{1}{2}\mathbb{E}_2\right\}, \end{gather}

(4.12b)\begin{gather} \mathbb{E}_1 = [1 - \exp(-t)], \end{gather}

(4.12c)\begin{gather} \mathbb{E}_2 = [1 - \exp(-2t)]. \end{gather}

4.3. Preliminary discussion

Comparing the position variance arising from PL and VL (Maxwellian and impulse initial condition), viz.

(4.13a)\begin{gather} \langle (x - \langle x \rangle)^2 \rangle = 2 St\,Pe^{-1} t, \end{gather}

(4.13b)\begin{gather} \langle (x - \langle x \rangle)^2 \rangle = 2 St\,Pe^{-1} \{ t - \mathbb{E}_1\}, \end{gather}

(4.13c)\begin{gather} \langle (x - \langle x \rangle)^2 \rangle = 2 St\,Pe^{-1} \left\{ t - 2\mathbb{E}_1 + \frac{1}{2}\mathbb{E}_2\right\}, \end{gather}

demonstrates the canonical form prefaced in § 4.1. Namely, PL gives pure diffusion while VL gives diffusion in addition to exponential time scales $\mathbb {E}_i$ that arise from the drag force and the initial condition. At sufficiently long times, VL dispersion converges to pure diffusion. If the particle velocities are initialized with the steady Maxwellian distribution in VL, then the velocity variance has no time dependence and the velocity autocovariance function depends upon a single time lag $s_1$. As a direct result, the integral in (A 5) simplifies and the classic ballistic–diffusive solution of Taylor (Reference Taylor1922) is recovered; see (4.13b). By contrast, if the particles are initialized with zero velocity, then the granular temperature evolves in time and additional time scales are introduced into the ballistic regime that act to suppress dispersion; see (4.13c). Thus, the solution in (4.13c) demonstrates the impact of a temporally evolving velocity variance on particle dispersion with the VL framework. Since ballistic dispersion terms scale with $\tau _p^2$, they are most significant for highly inertial particles.

The solutions given in (4.13b) and (4.13c) are first verified against numerical integration of (4.2a) and (4.2b) with $St=5$ and $Pe=2$. Specifically, the analytical distribution $P(x;t) = \mathcal {N}[0,\sigma ^2]$, with $\sigma ^2$ given in (4.13a)–(4.13c) is compared to distributions obtained on the binned computational domain; see figure 1. Strong agreement is observed between the analytical solutions and numerical results for both the impulse (figure 1a,b) and Maxwellian (figure 1c,d) velocity conditions. Additionally, the ensemble averaged velocity variance and position variance, at $Pe=5$ and $St = 1, 10, 100$, are compared to (4.8) and (4.13c) in figure 2. As the Stokes number is reduced, particles behave more like fluid tracers and their velocity fluctuations undergo rapid relaxation, thereby yielding reduced dispersion. The classic ballistic $(t \leq 1)$ and diffusive $(t \geq 5)$ regimes may be seen in the dispersion curves in figure 2.

Figure 2.

$(a)$ Analytical and numerical predictions for velocity variance with the VL framework and an impulse velocity condition. $(b)$ Analytical and numerical predictions for position variance with the VL framework and an impulse velocity condition. $Pe=5$ for each Stokes number.

5. Force Langevin (FL)

As noted in § 1, force Langevin treatments are far less pervasive within the particle-laden flow literature; the works by Esteghamatian et al. (Reference Esteghamatian, Euzenat, Hammouti, Lance and Wachs2018) and Andrews et al. (Reference Andrews, Loezos and Sundaresan2005) being noted exceptions. Since the application of FL schemes has been empirical in nature, the velocity variance and dispersion resulting from such a framework has not been emphasized in the literature. While not a force Langevin, the fluctuating fluid velocity model of Pozorski & Apte (Reference Pozorski and Apte2009) has some interesting analogies to the FL scheme. Specifically, an additional SDE is introduced for the fluctuating fluid velocity which is commonly taken to follow an Ornstein–Uhlenbeck process. We consider a similar approach for the force Langevin, motivated by the Gaussian drag forces reported by PR–DNS studies. While we do not consider the task here, solutions obtained within this section for FL may be readily adapted to fluid VL schemes.

For the FL treatment we begin by splitting the total drag force into a mean and fluctuating component $\boldsymbol {F}_{d,i} = \left \langle \boldsymbol {F}_{d}\right \rangle + \boldsymbol {F}^{\prime }_{d,i}$ (Esteghamatian et al. Reference Esteghamatian, Euzenat, Hammouti, Lance and Wachs2018). Random increments introduced on the particle force lead to

(5.1a)\begin{gather} \mathrm{d}\boldsymbol{X}_{p,i} = \boldsymbol{U}_{p,i}\,\mathrm{d}t \end{gather}

(5.1b)\begin{gather} \mathrm{d}\boldsymbol{U}_{p,i} = -\frac{1}{\tau_p} \boldsymbol{U}_{p,i}\,\mathrm{d}t + \frac{\boldsymbol{F}_{d,i}}{m_{p,i}}\,\mathrm{d}t, \end{gather}

(5.1c)\begin{gather} \mathrm{d}\boldsymbol{F}_{d,i} = -\frac{1}{\tau_F} \boldsymbol{F}_{d,i}\,\mathrm{d}t + \hat{\boldsymbol{b}}\, \mathrm{d}\boldsymbol{W}_t, \end{gather}

where the $(\cdot )^{\prime }$ notation is again omitted for readability. In order to keep the fluctuating force distribution bounded, a restorative component must be present in the force SDE. Similar to Pozorski & Apte (Reference Pozorski and Apte2009) and Esteghamatian et al. (Reference Esteghamatian, Euzenat, Hammouti, Lance and Wachs2018), we consider an Ornstein–Uhlenbeck process with relaxation time scale $\tau _F$. Setting the diffusion tensor as $\bar {b}_{ij} = \sqrt {2Dm_{p,i}^2/\tau _p^2 \tau _F^2} \delta _{ij}$, we arrive at a dimensionless form of

(5.2a)\begin{gather} \mathrm{d}\boldsymbol{X}_{p,i} = St \boldsymbol{U}_{p,i}\,\mathrm{d}t, \end{gather}

(5.2b)\begin{gather} \mathrm{d}\boldsymbol{U}_{p,i} = -\boldsymbol{U}_{p,i}\,\mathrm{d}t + \boldsymbol{F}_{d,i}\,\mathrm{d}t, \end{gather}

(5.2c)\begin{gather} \mathrm{d}\boldsymbol{F}_{d,i} = -St_F \boldsymbol{F}_{d,i}\,\mathrm{d}t + \sqrt{2 St_F^2\left(Pe\,St \right)^{-1}}\,\mathrm{d}\boldsymbol{W}_t, \end{gather}

where $St_F = \tau _p/\tau _F$ relates the Stokes response time to the integral time scale of the fluctuating drag force. For large fluctuating Stokes numbers $St_F$, the fluctuating drag forces have a short lifespan in comparison to the particle response time. Under these conditions, one would expect the FL treatment to be similar to the VL treatment, where the fluctuating drag force contains no memory; see (4.1b). In the opposite limit of small $St_F$, the fluctuating forces persist for time scales larger than the Stokes response time and significant departure from the VL treatment should be expected. The corresponding one-dimensional Fokker–Planck equation is given by

(5.3)\begin{equation} \frac{\partial P}{\partial t} + St \frac{\partial }{\partial x} (uP) + \frac{\partial }{\partial u} [(\,f - u)P] - St_F \frac{\partial }{\partial f} (\,f P)= St_F^2(Pe\,St)^{-1} \frac{\partial^{2} P}{\partial {f}^{2}}. \end{equation}

Details derived for the Ornstein–Uhlenbeck process in § 4 ((4.6a)–(4.8)) may be directly applied to the fluctuating force – i.e. a mean of $f_0\exp (-St_F t)$ and variance of $St_F(Pe\,St )^{-1}[1 - \exp (-2St_F t)]$.

In the limit of stationary particle assemblies, previous works have shown that a drag force distribution with characteristic variance will be present; see cited works in § 1. For this reason, we take the fluctuating force to be fully developed $\mathcal {N}[0,St_F(Pe\,St )^{-1}]$ and do not consider the temporal growth of its variance. Thus, in the present derivation, the fluctuating force distribution is stationary with a constant variance and relaxation time. In general, FL is not restricted to a constant variance or relaxation time, and previous results with particle-resolved simulations suggest that these parameters will vary within a general multiphase flow (Koch & Sangani Reference Koch and Sangani1999; Huang et al. Reference Huang, Wang, Zhou and Li2017; Esteghamatian et al. Reference Esteghamatian, Euzenat, Hammouti, Lance and Wachs2018). Connections between the solutions derived here and PR–DNS are made in § 7, with emphasis being given to how the model coefficients may be informed.

For the sake of completeness, we note that the diffusion tensor has again been set to yield the canonical form for dispersion. If the diffusion tensor is set in the manner discussed at the end of § 4.1$\sqrt {2\sigma _F^2/ \tau _F}$, then $\sigma _F^2$ characterizes the steady variance in fluctuating force. For this case, the solutions derived in this section may be cast in dimensional form by replacing $D$ with $\sigma _F^2\tau _p^2\tau _F/m_{p,i}^2$.

5.1. Fluctuating velocity (FL)

Integrating (5.3) over the spatial $(x)$ and fluctuating force $(\,f)$ coordinates allows the marginal velocity distribution to be obtained

(5.4)\begin{equation} \frac{\partial P(u;t|v,s)}{\partial t} + \frac{\partial }{\partial u} [(\langle\, f|u\rangle - u)P] = 0. \end{equation}

The first two moments of (5.4) are given by

(5.5a)\begin{gather} \frac{\textrm{d} \langle u \rangle}{\textrm{d} t} = - \langle u \rangle + \langle\langle\, f|u\rangle\rangle, \end{gather}

(5.5b)\begin{gather} \frac{\textrm{d} \langle u^2 \rangle}{\textrm{d} t} = -2 \langle u^2 \rangle + 2\langle u\langle\, f|u\rangle\rangle. \end{gather}

Comparing the granular temperature balance with FL (5.5b) to VL (4.6b) shows that the functional form is unaltered but the source term is described by the force–velocity covariance with FL. This is a direct result of modelling the force statistics via FL as well as the substitution of the instantaneous particle velocity into Stokes drag. Since the fluctuating force is defined to be stationary with variance $St_F(Pe\,St )^{-1}$, the unconditional mean of force $\langle \langle\, f|u \rangle \rangle$ is zero, leading to zero mean velocity $\langle u \rangle$. We quickly recognize that $\langle u\langle\, f|u \rangle \rangle$ must be handled via Lagrangian variables in a similar manner to $\langle x \langle u | x \rangle \rangle$ in § 4.2. Starting with Newton's second law in dimensionless form

(5.6)\begin{equation} \frac{\textrm{d} U_{p,i}}{\textrm{d} t} = -U_{p,i} + F_{p,i}, \end{equation}

the particle velocity with time is given by

(5.7)\begin{equation} U_{p,i}(t) = \exp(-t) \int_0^t \exp(t^{\prime}) F_{p,i}(t^{\prime}) \, \mathrm{d}t^{\prime}, \end{equation}

leading to

(5.8)\begin{equation} \langle u\langle\, f|u\rangle\rangle = \exp(-t) \int_0^t \exp(t^{\prime}) \langle F_{p,i}(t^{\prime}) F_{p,i}(t) \rangle \, \mathrm{d}t^{\prime}. \end{equation}

Since the fluctuating force is fully developed, the two-time Lagrangian force correlation $\left \langle F_{p,i}(t^{\prime })F_{p,i}(t) \right \rangle$ only depends upon a single time lag $K_{uu}(s_1) = St_F (Pe\,St )^{-1}$$\exp (-St_F s_1)$; $s_1 = t - t^{\prime }$. Therefore, the source term in (5.5b) follows

(5.9)\begin{equation} \langle u\langle\, f|u\rangle\rangle = St_F(Pe\,St )^{-1}\exp(-t) \int_0^t \exp(t^{\prime})\exp(-St_F(t-t^{\prime})) \, \mathrm{d}t^{\prime}. \end{equation}

Integrating (5.9), we arrive at

(5.10)\begin{equation} \langle u\langle\, f|u\rangle\rangle = (Pe\,St )^{-1} \frac{St_F}{St_F + 1} [1 - \exp(-(St_F + 1) t) ]. \end{equation}

Substituting (5.10) into (5.5b) allows an integrating factor solution for $\langle u^2 \rangle$ to be found

(5.11a)\begin{gather} \langle u^2 \rangle = 2 (Pe\,St)^{-1} \frac{St_F}{St_F + 1} \left\{ \frac{1}{2} \mathbb{E}_2 + \frac{1}{St_F - 1} (\mathbb{E}_2 - \mathbb{E}_3 ) \right\}, \end{gather}

(5.11b)\begin{gather} \mathbb{E}_3 = [ 1 - \exp(-(St_F + 1) t)]. \end{gather}

For the special case that $St_F = 1$ ($\tau _F = \tau _p$), the integrating factor for (5.5b) exactly offsets the exponential term in (5.10) and a simplified solution is obtained

(5.12)\begin{equation} \langle u^2 \rangle = (Pe\,St)^{-1} \left\{ \frac{1}{2} \mathbb{E}_2 - t\exp(-2t) \right\}. \end{equation}

Granular temperatures computed from simulations with the FL scheme ((5.2a)–(5.2c)) are in excellent agreement with (5.11a)–(5.12); see figure 3(a). In figure 3(a) the Stokes number is held fixed while the fluctuating Stokes number and Péclet number are varied. For demonstration purposes, the Péclet numbers are chosen such that the force variance is conserved in all three simulations, thereby allowing the effect of the force memory, $St_F$, to be probed. Reducing $St_F$, at a fixed force variance, causes the force fluctuations to decay at a slower rate and ultimately giving rise to enhanced granular temperatures. Therefore, fluctuating drag forces with significant inertia are characterized by small $St_F$. It is important to note that the VL scheme inherently assumes that the fluctuating drag force is non-inertial, just as the PL scheme assumes the particles themselves are non-inertial. This hierarchy further motivates the Stokes number analogy employed for the fluctuating force and allows a connection to be drawn to the VL scheme. Specifically, FL predicts larger granular temperatures as the fluctuating forces become more inertial, smaller $St_F$, just as the VL scheme predicts greater dispersion as the particles become more inertial, larger $St$.

Figure 3.

$(a)$ Analytical and numerical prediction for transient growth in the velocity variance with the FL framework and a $\delta (u)$ initial condition. $(b)$ Analytical and numerical predictions for the same system but with particle velocities initialized as a Maxwellian. In all cases, $St=10$.

It is important to note that the solution obtained in (5.11a) converges to a steady velocity variance of $\langle u^2 \rangle _{\infty } = (Pe\,St )^{-1} St_F /(St_F + 1)$. However, the FL framework yields a time dependent granular temperature even if the system is initialized with the steady velocity variance. By contrast, the VL framework yields a constant granular temperature if the system is initialized with the steady velocity variance; see (4.8) with an initial condition of $(Pe\,St )^{-1}$. To emphasize this point, we consider imposing an IC of $\langle u^2 \rangle _{\infty }$ when solving (5.5b) to obtain

(5.13)\begin{equation} \langle u^2 \rangle = 2 \left(Pe\,St \right)^{-1} \frac{St_F}{St_F + 1} \left\{ \frac{1}{2} - \frac{1}{St_F - 1} (\mathbb{E}_2 - \mathbb{E}_3) \right\}, \end{equation}

and for the case $St_F=1$

(5.14)\begin{equation} \langle u^2 \rangle = (Pe\,St)^{-1} \left\{ \frac{1}{2} - t\exp(-2t) \right\}. \end{equation}

Comparing (5.13) to numerical results shows the time dependence of the granular temperature; see figure 3(b). At long times, the system converges to a velocity variance of $(Pe\,St )^{-1} St_F /(St_F + 1)$. However, an initial reduction in granular temperature is observed, due to the $\mathbb {E}_3$ time scale (force–velocity covariance), and the recovery time depends upon the fluctuating Stokes number (competition between $\mathbb {E}_3$ and $\mathbb {E}_2$ terms).

The aforementioned results may be physically understood by recognizing that the Maxwellian distributions employed for the fluctuating force and particle velocity are sampled independently when the system is initialized. Therefore, the fluctuating forces, which drive the granular temperature, are not correlated with the particle velocities when the system begins evolving. Since the force–velocity correlations develop over time scales associated with the fluctuating Stokes number, the effect of the initial force condition persists longer as $St_F$ is decreased and leads to the accentuated dips in figure 3(b). We stress that the initial reduction in granular temperature, due to force–velocity covariances, is indeed physical and was observed in the PR–DNS simulations of Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016) with freely evolving particles (see green curve in figure 5a of Tenneti et al. Reference Tenneti, Mehrabadi and Subramaniam2016). As we increase the fluctuating Stokes number, dips in granular temperature are less prominent; and at sufficiently large $St_F$, the granular temperature remains constant. This behaviour is consistent with what would be obtained by a VL framework, which treats the fluctuating forces as non-inertial. In summary, the fluctuating Stokes number characterizes the inertia of the fluctuating drag force and gauges the departure of FL from VL. The potential dependencies of $St_F$ and connections with various types of flows are discussed in greater detail in § 6.

5.2. Dispersion (FL)

Integrating (5.3) over the fluctuating force $f$ and velocity $u$ coordinates allows the marginal spatial distribution to be obtained

(5.15)\begin{equation} \frac{\partial P(x;t|y,s)}{\partial t} + St \frac{\partial }{\partial x} [ \langle u|x \rangle P ] = 0. \end{equation}

The first two moments of (5.15) are given by

(5.16a)\begin{gather} \frac{\textrm{d} \langle x \rangle}{\textrm{d} t} = St \langle\langle u|x\rangle\rangle, \end{gather}

(5.16b)\begin{gather} \frac{\textrm{d} \langle x^2 \rangle}{\textrm{d} t} = 2St \langle x\langle u|x\rangle\rangle. \end{gather}

The unconditional velocity $\langle \langle u|x\rangle \rangle$ was shown to be zero in § 5.1 (stationary distribution for $f$), and thus, the mean position will also be zero. The $\left \langle x\left \langle u|x\right \rangle \right \rangle$ term is addressed in the same manner as appendix A.1 and the details are given in appendix A.3. Additionally, for the sake of completeness, the velocity autocovariance function for a fully developed velocity distribution is derived in appendix A.4. Integration of the velocity autocovariance function gives a closed-form solution for position variance

(5.17)\begin{equation} \langle x^2\rangle = 2 St\,Pe^{-1}\{ t - C_1 \mathbb{E}_1 + C_2 \mathbb{E}_2 - C_3 \mathbb{E}_3 + C_4 \mathbb{E}_4\}, \end{equation}

where $\mathbb {E}_4 = [1 - \exp (-St_F t ) ]$ and the constants $(C_i)$ are given by

(5.18a)\begin{gather} C_1 = \frac{2 St_F - 1}{St_F - 1}, \end{gather}

(5.18b)\begin{gather} C_2 = \frac{St_F}{2(St_F - 1)}, \end{gather}

(5.18c)\begin{gather} C_3 = \frac{1}{(St_F + 1)(St_F - 1)}, \end{gather}

(5.18d)\begin{gather} C_4 = \frac{1}{St_F(St_F-1)}. \end{gather}

Repeating the process for $St_F=1$ leads to

(5.19)\begin{equation} \langle x^2\rangle = 2 St\,Pe^{-1} \left\{ t - 3 \mathbb{E}_1 + \frac{3}{4} \mathbb{E}_2 + t \exp(-t) - \frac{1}{2} t \exp(-2t) \right\}. \end{equation}

The analytical solution for dispersion, Gaussian distribution $P(x;t) = \mathcal {N}[0,\sigma ^2]$ with the variance specified by (5.17), is in strong agreement with numerical results; see figure 4. The conditions in figure 4 are given by: $St=100$, $Pe=0.1$, and $St_F=0.3$. Ensemble averaging of the position variances, at the same conditions considered in figure 3 (constant force variance), yields the curves in figure 5. Examination of figures 5 and 3 shows that, for a constant force variance, smaller $St_F$ lead to larger granular temperature and greater dispersion. Furthermore, it becomes clear that the force variance and the relaxation time scale allow the FL framework to control the granular temperature and dispersion. By contrast, the VL framework with a mean drag closure only allows for control of the granular temperature (dispersion is a result of the specified velocity variance), and the PL framework only allows for control of the long time diffusion behaviour.

Figure 4.

(a,b) Analytical and numerical predictions for particle dispersion with the FL framework. An impulse velocity condition was utilized in all cases and the conditions are given by $St = 100$, $Pe = 0.1$, $St_F=0.3$. Particles (left half of images) are coloured by the magnitude of their velocity in the $x$-direction. $(a)\,\, t = 0.500, (b)\,\, t = 1.000.$

Figure 5.

Analytical and numerical predictions for dispersion with the FL framework. An impulse velocity condition was utilized in all cases and $St=10$.

7. Model development for Euler–Lagrange methods

In general, FL requires closure of two model parameters: (i) the drift vector (integral time scale of fluctuating force; $\tau _F$) and (ii) the diffusion tensor (force variances; $\sigma _F^2$). Examining the fluid–particle force statistics within PR–DNS of freely evolving particles allows both parameters to be probed. The relaxation time scale(s) may be examined via the two-time force autocovariance; see appendix B here and equation (28) in Esteghamatian et al. (Reference Esteghamatian, Euzenat, Hammouti, Lance and Wachs2018). Preliminary examinations of the force autocovariance have been made by Esteghamatian et al. (Reference Esteghamatian, Bernard, Lance, Hammouti and Wachs2017), but correlations for this quantity are not readily available in the literature. Similarly, the drag force variance has been reported in a variety of PR–DNS works (Koch & Sangani Reference Koch and Sangani1999; Wylie et al. Reference Wylie, Koch and Ladd2003; Akiki et al. Reference Akiki, Jackson and Balachandar2016; Huang et al. Reference Huang, Wang, Zhou and Li2017), but again, correlations for this quantity are not well established. Since force time scales extracted from PR–DNS may lead to small fluctuating Stokes numbers, where differences between FL and VL occur, it appears that the FL framework may have a significant impact on certain multiphase flows. Therefore, well characterized inputs to the FL framework, obtained from statistical analysis of the drag forces present in freely evolving suspensions, would be useful.

There are a variety of means by which the present analysis may be extended. For large Stokes numbers, Koch & Sangani (Reference Koch and Sangani1999) approximate the time scale for the drag force with the mean-free time between successive collisions (Chapman, Cowling & Burnett Reference Chapman, Cowling and Burnett1970)

(7.1)\begin{equation} \tau_{col} = \frac{d_p}{24 \phi \chi} \sqrt{\frac{{\rm \pi}}{T}}, \end{equation}

where $\chi$ is the radial distribution function (RDF) at contact. The mean-free time $\tau _{col}$, computed with the correlation of Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016) and the RDF of Ma & Ahmadi (Reference Ma and Ahmadi1988), is observed to be in good agreement with the $\tau _F$ time scales extracted from PR–DNS; see table 1. It is worth noting that $\tau _{col}$ is ${\sim }13\,\%$ larger than $\tau _F$ for the $Re_m=20$ case and actually improves agreement between FL theory and the correlation of Tenneti et al. (Reference Tenneti, Mehrabadi and Subramaniam2016). Therefore, taking $\tau _{col} \approx \tau _F$ and correlating the force variances $\sigma _F^2$ in static particle assemblies may constitute an appropriate model for large Stokes numbers. Physically, larger granular temperatures will lead to rapid alteration of particle structuring and shorter force memories. While the data in table 1 only constitute two points, values of $St_F$ are observed to increase with increasing $Re_T$. Informing the force time scale with the granular temperature would provide valuable connections between the velocity and force equations, which are presently one-way coupled. Since the steady, fluid-mediated, granular temperature is defined by the $(\rho _f/\rho _p, \ \phi , \ Re_m )$ triple point (Tenneti et al. Reference Tenneti, Mehrabadi and Subramaniam2016), an implicit relation between the flow conditions and FL inputs is expected. Under this modelling paradigm, space–time locality should also be explored since variation in the flow conditions will lead to variation in the FL inputs. Additionally, sub-grid, pseudo-turbulent fluid velocity fluctuations (Mehrabadi et al. Reference Mehrabadi, Tenneti, Garg and Subramaniam2015) may prove useful for informing the drag variance. Specifically, the use of rotation matrices, as in Peng et al. (Reference Peng, Kong, Zhou, Sun, Passalacqua, Subramaniam and Fox2019), may allow for anisotropic behaviour. For this case, fluctuating forces in a body-fixed coordinate system aligned with the slip velocity would be constructed and then rotated back to the laboratory reference frame. Such a framework could account for stochastic forces in the direction of mean slip, as well as the transverse directions.