A model for drift velocity mediated scalar eddy diffusivity in homogeneous turbulent flows

Abstract

Low Stokes number particles at dilute concentrations in turbulent flows can reasonably be approximated as passive scalars. The added presence of a drift velocity due to buoyancy or gravity when considering the transport of such passive scalars can reduce the turbulent dispersion of the scalar via a diminution of the eddy diffusivity. In this work, we propose a model to describe this decay and use a recently developed technique to accurately and efficiently measure the eddy diffusivity using Eulerian fields and quantities. We then show a correspondence between this method and standard Lagrangian definitions of diffusivity and collect data across a range of drift velocities and Reynolds numbers. The proposed model agrees with data from these direct numerical simulations, offers some improvement to previous models in describing other computational and experimental data and satisfies theoretical constraints that are independent of Reynolds number.

1. Introduction

The motion of particles in a turbulent carrier fluid subject to exogenous body forces manifests in myriad applications such as ash settling (Mingotti & Woods 2020), energy production (Ishii 1977; Guet & Ooms 2006) and bubbly wakes (Carrica et al. 1999). These particles may experience multiphysics including nucleation, coalescence, dissolution and growth, but their dispersion by energetic eddies of the background flow naturally lends itself to the field of turbulence modelling. In the framework of population balance equations, as reviewed in Shiea et al. (2020), particles can be segregated into size classes associated with a drift velocity, $u_d$, relative to the background flow with root-mean-square velocity $u_{rms}$. In nature, $u_d/u_{rms}$ can be large; a Hinze-scale 1 mm bubble, the largest size that should not undergo breakup, may rise at $12\,{\rm cm}\,{\rm s}^{-1}$, while a characteristic upper ocean turbulence velocity is closer to $2\,{\rm cm}\,{\rm s}^{-1}$ (Detsch & Harris 1989; Deane & Stokes 2002; D'Asaro 2014).

For dilute, low Stokes number ($St$) particle-laden flows with negligible particle–particle interactions, the particle concentration – or void concentration – fields can be described with the passive scalar advection–diffusion equation (Moraga et al. 2003). A low $St$ particle flow is one for which the particle momentum relaxation time is far smaller than the flow time scale, and ‘dilute’ flows should have particle volume fractions ($\phi$) less than ${\sim }O({10^{-3}}$) (Brandt & Coletti 2022). Each class can be evolved separately and forces like buoyancy and gravity are accounted for with additional velocity components from relations like the drag law of Schiller & Naumann (1935). For sedimenting particles, this drift is aligned with the gravity vector; in the case of bubbly flows, it is generally anti-parallel.

Solving for full-resolution scalar evolution when only a mean state is required to determine quantities of engineering interest is prohibitively expensive. Averaging to find the mean occurs in homogeneous spatio-temporal dimensions, defined through Reynolds averaging or filtering in the large-eddy simulation context. When averaging is applied to the Navier–Stokes and advection–diffusion equations that govern the evolution of scalar-laden fluid flows, the turbulent scalar flux, given by $\overline {u_i c}$, appears. Here, $u_i$ and $c$ are the fluctuating velocity and scalar fields, respectively, and $\bar {\bullet }$ is an average. This term represents unresolved scalar transport by turbulence and its closure is essential to making transport simulations tractable.

A common model form for the turbulent scalar flux is gradient diffusion, written analogously to Fickian diffusion as $\overline {u_i c} = - D_{ij} ({\partial \overline {C}}/{\partial x_j})$, with $D_{ij}$ representing eddy diffusivity and $C$ the full scalar field. Foundational work has shown eddy diffusivity decays with increased drift, but few extant models for capturing the flux are algebraic closed-form expressions (Yudine 1959; Moraga et al. 2003; Reeks 2021). An exception is Csanady (1963), which derives an expression for particle diffusivity scaling as a function of $u_d$ from theoretical arguments about the form and relevant parameters of the velocity autocorrelation, but Squires & Eaton (1991) showed disparities between it and measured experimental and computational turbulent data.

As Csanady (1963) serves as a starting point for more complex models (e.g. Wang & Stock 1993) and particle-laden turbulence is still a ripe topic (e.g. Berk & Coletti 2021), we wish to revisit this problem using the recently developed macroscopic forcing method (MFM) to calculate eddy diffusivity (Mani & Park 2021; Shirian & Mani 2022). In general, the eddy diffusivity is a non-local spatial and temporal operator, but when there is separation of scales between the large-eddy size and the scalar cloud size, as is relevant for this problem, measurement of a single local coefficient of eddy diffusivity, denoted $D^0$ and sometimes called a dispersion coefficient, suffices. Mani & Park (2021) showed that $D^0$ is the leading-order moment of the full eddy diffusivity operator and Shende & Mani (2022) used it to predict scalar transport. Furthermore, it is the first term in the Kramers–Moyal approximation of the integro-differential kernel that defines the true eddy diffusivity operator.

In this work, we propose a model using measurements of eddy diffusivity by MFM in numerical simulations of scalar fields driven by homogeneous, isotropic turbulence (HIT) over a range of drift velocities. The MFM has more favourable computational costs when compared with a similar Lagrangian method for this a priori modelling approach. The proposed model better fits empirical data and theoretical constraints for tracers subject to drift.

2. Model problem

Consider a triply periodic cubic domain of HIT. When $u_d = 0$, the diffusivities along the three principal axes – $D^0_{11}$, $D^0_{22}$ and $D^0_{33}$ – are all equal to $D^0$. If we impose a non-negative drift velocity in the $x_1$ direction, the symmetry of the set-up is broken such that only the diffusivities in the $x_2$ and $x_3$ directions are equal and all are affected by drift. While the derivation herein reaches conclusions similar to others (e.g. Csanady 1963; Squires & Eaton 1991; Mazzitelli & Lohse 2004), we distinctly adopt an Eulerian perspective.

For $u_d = 0$, we begin with the Lagrangian formulation of eddy diffusivity of Taylor (1922), with the scalar represented by tracer particles with position $X_i(t)$ and velocity $V_i(t)$. Taylor (1922) writes the eddy diffusivity in the $x_1$ direction of a ensemble of such particles, in the long time limit. If the velocity field is statistically stationary, this can be expressed as

(2.1)\begin{equation} D^0_{11} = \int^\infty_0 \langle (V_1(\tau+t) V_1(t)\rangle\,{\rm d}\tau, \end{equation}

where $\langle \bullet \rangle$ represents an average over the ensemble and $\tau$ is a temporal offset. In the high Péclet number ($Pe$) limit, where molecular diffusivity is far smaller than eddy diffusivity, $V_i(t)$ is the flow velocity at a given particle's position. Therefore, we can express the underlying Eulerian velocity field component, $u_1$, as a function of the full three-dimensional tracer position, $\boldsymbol {X}$. This yields

(2.2)\begin{equation} D^0_{11} = \int^\infty_0 \langle u_1(\boldsymbol{X}(t),t) u_1(\boldsymbol{X}(t+\tau),t+\tau) \rangle\,{\rm d}\tau. \end{equation}

We define a standard complementary characteristic turbulence length scale as

(2.3)\begin{equation} L_{11} = \frac{1}{u^2_{rms}} \int^\infty_0 \langle u_1(x_1, x_2, x_3,t) u_1(x_1+r,x_2,x_3,t) \rangle\,{\rm d}r, \end{equation}

where $r$ is a spatial offset, $\langle \bullet \rangle$ now represents an average over all independent variables. Here, $u_{rms}$ normalizes the velocity autocorrelation at zero displacement, namely $\sqrt {\langle u_1u_1 \rangle }$.

Classical Brownian motion analysis posits that scalar diffusivity at infinite time scales as $u_{rms}^2 \tau _k$, where $\tau _k \equiv L_k/u_{rms}$ denotes some time scale of the underlying flow development and $L_k = 2L_{11}$ is the large-eddy length scale. Consider the scalar field now with some constant drift velocity, $u_d$, with respect to the flow: for very small $u_d$, the fundamental turbulence statistics felt by a scalar parcel are not largely affected and the classical model holds. In the opposite limit that $u_{rms} \ll u_d$ such that $L_{11}/u_d \ll \tau _k$, however, this picture is not appropriate.

In this limit, the scalar field transits the turbulence field at a very fast time scale $\tau _d = L_{11}/u_d$, such that scalar particles drift before the local flow has evolved. The turbulence can therefore be considered frozen compared to the evolution of the scalar field for computing (2.2), and a leading-order approximation for the differential translation of a scalar parcel is $\Delta X_1 = {\rm d}r \approx u_d\,{\rm d}t$. Thus, the addition of a very large drift is equivalent to examining a translating inertial coordinate frame with respect to the frozen flow, akin to Taylor's hypothesis. These premises allow us, for very large $u_d$, to rewrite (2.2) as

(2.4)\begin{equation} \lim_{u_d \to \infty} D^0_{11} = \int^\infty_0 \langle u_1(x_1, x_2, x_3,t) u_1(x_1+u_d\tau,x_2,x_3,t) \rangle\,{\rm d}\tau. \end{equation}

We have here posited that the change in particle position is purely due to drift, and over the time of $O(L_{11}/u_d) \ll \tau _k$ where the kernel is non-zero, the field $u_1$ does not change. If we now use a change of variable between the displacement and the drift velocity, we can write

(2.5)\begin{equation} \lim_{u_d \to \infty} D^0_{11} = \int^\infty_0 \langle u_1(x_1, x_2, x_3,t) u_1(x_1+r,x_2,x_3,t) \rangle d(r/u_d) = u^2_{rms} \frac{L_{11}}{u_d}. \end{equation}

Here, $L_{11}$ is the correlation length scale in (2.3) and it requires the spatial correlation drop to zero in the domain. For the transverse diffusivities, $D^0_{22} = D^0_{33}$. We can therefore similarly write

(2.6)\begin{equation} \lim_{u_d \to \infty} D^0_{22} = \int^\infty_0 \langle u_2(x_1, x_2, x_3,t) u_2(x_1+r,x_2,x_3,t) \rangle d(r/u_d) = u^2_{rms} \frac{L_{22}}{u_d}. \end{equation}

Appealing to the isotropy of the underlying velocity fields, we also can show that

(2.7)\begin{equation} \lim_{u_d \to \infty} D^0_{22} = \int^\infty_0 \langle u_1(x_1, x_2, x_3,t) u_1(x_1,x_2+r,x_3,t) \rangle d(r/u_d) = u^2_{rms} \frac{L_{22}}{u_d}. \end{equation}

This asymptotic limit jibes with intuition, as a particle with infinite $u_d$ samples zero mean velocity over every time horizon in homogeneous turbulence. This decay of $D^0_{ii}$ with increasing $u_d/u_{rms}$, the ‘crossing trajectories’ effect of Yudine (1959) and Csanady (1963), persists in relatively high Reynolds number ($Re$) bubble experiments (Mathai et al. 2018).

We seek a model form that analytically matches the exact behaviour at the infinite and zero drift velocity limits. In the intermediate regime, the model should be a smooth and monotonic function of the drift velocity, so a simple model form that captures the asymptotic limits of infinite and zero drift is

(2.8)\begin{equation} \frac{D^0_{ii}}{D^0} \approx \Bigg( 1 + \left(\frac{u_dD^0}{L_{ii}u^2_{rms}}\right)^{\alpha_{ii}} \Bigg)^{{-}1/\alpha_{ii}} = \left( 1 + \left(\frac{u_d}{u^*_{ii}}\right)^{\alpha_{ii}} \right)^{{-}1/\alpha_{ii}}. \end{equation}

Here, $\alpha _{ii}$ is a free parameter and $u_{ii}^*$ is an Eulerian ‘diffusion velocity’ that competes with the drift and is defined using Eulerian measures as

(2.9)\begin{equation} u^*_{ii} = \frac{u^2_{rms} L_{ii}}{D^0}. \end{equation}

In incompressible flows, $L_{22} = L_{33} = L_{11}/2$ (Csanady 1963). The eddy diffusivity at large drift velocity can be written as $D^0_{ii}/D^0 = u^*_{ii}/u_d$. In contrast, Csanady (1963) proposed

(2.10)\begin{equation} \frac{D^0_{ii}}{D^0} \approx \Bigg( 1 + \left(\frac{\gamma_{ii}~\beta~ u_d}{u_{rms}}\right)^2 \Bigg)^{{-}1/2}, \end{equation}

where $\beta$ is defined with Lagrangian and Eulerian statistics and $\gamma _{22} = 2\gamma _{11} = 2$.

The proposed model form of (2.8) is very similar to (2.10), suggesting the space of candidate functions appropriate for this problem is narrow. As the adopted problem set-up (i.e. one-way coupled, dilute particle limit, negligible inertial effects) is a simplification, even such an approximate model should offer substantial predictive value.

3. Numerical set-up

Direct numerical simulation (DNS) of forced incompressible HIT in a triply periodic box of edge length $L_{box}$ provides the data for this test. Solved with continuity are the incompressible Navier–Stokes momentum equations with a linear forcing term. They are written as

(3.1)\begin{equation} \frac{\partial{u_i}}{\partial{t}} + \frac{\partial{u_i u_j}}{\partial x_j} ={-} \frac{1}{\rho} \frac{\partial p}{\partial x_i} + \nu \frac{\partial^2 u_i}{\partial x_j \partial x_j} + A\widetilde{u_i}, \end{equation}

where $A$ is a controller that maintains the turbulent kinetic energy, $k \equiv u_i u_i/2$, at a prescribed level, $\rho$ is the fluid density and $\nu$ is the kinematic viscosity. This approach is similar to Bassenne et al. (2016), but with a high-pass filtered velocity described in spectral space as $\widehat {\widetilde {u_i}} = G({|\boldsymbol {{k}}|})\widehat {u_i}$, where $|\boldsymbol {k}|$ is the wavenumber magnitude and

(3.2)\begin{equation} G({|\boldsymbol{k}|}) = \begin{cases} 0, & {|\boldsymbol{k}|} \leq 2, \\ \dfrac{1}{2} - \dfrac{1}{2}\cos({\rm \pi} ({|\boldsymbol{k}|}-2)), & 2 < {|\boldsymbol{k}|} \leq 3,\\ 1, & {|\boldsymbol{k}|} > 3. \end{cases} \end{equation}

In figure 1, plotted instantaneous energy spectra show decay at the largest scales for all considered $Re$. The values of $L_k$, $L_{11}$, $L_{22}$ and $L_{33}$ are all less than $L_{box}/2$ by construction due to the filtered energy injection. As such, the simulations are not influenced by the cubic box shape or orientation, unlike the standard linear forcing method, as in Rosales & Meneveau (2005).

Figure 1. Instantaneous unnormalized velocity autocorrelations (a) and normalized energy spectra $E(kL_{11})/u_{rms}^2L_{11}$ (b) for the cases in a $2{\rm \pi} ^3$ box, with filtered forcing preventing energy growth in the largest modes. Note that $u_{rms} = 1$ for all cases.

To quantify eddy diffusivity, the MFM formulation of Mani & Park (2021) is used with code adapted from Pouransari, Mortazavi & Mani (2016) for solving the governing equations on a staggered $N^3$ grid with finite volume operators and a Runge–Kutta time-advancement scheme. The MFM measures the response of the scalar field to an imposed macroscopic forcing by solving an additional equation for a scalar field with molecular diffusivity $D_m$. Following the procedure of Shirian & Mani (2022) for finding $D^0$, we decompose the scalar field as $C = \overline {C} + c$ and add a macroscopic source term, $s$, to the scalar equation so that $c$ is governed by

(3.3)\begin{equation} \frac{\partial{c}}{\partial{t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( ( \boldsymbol{u} + \boldsymbol{u}_d ) c \right) = D_m \nabla^2 c - (\boldsymbol{u} + \boldsymbol{u}_d ) \boldsymbol{\cdot} \boldsymbol{\nabla} \overline{C} + s(\boldsymbol{x},t). \end{equation}

Following the method of moments, the forcing maintains ${\partial \overline{C}}/{\partial x_i} = 1$ (Mani & Park 2021). Setting $i = 1$ allows measurement of the axial diffusivity in the direction of drift and $i = 2, 3$ allows quantification of transverse diffusivity in directions perpendicular to drift. Table 1 summarizes parameters swept to measure eddy diffusivity in the $x_1$ and $x_2$ directions.

Table 1. Summary of measured values relevant for the computation of $D^0_{11}$ and $D^0_{22}$.

In all cases, $\nu =D_m$ and $u_{rms} = \sqrt {2k/3} = 1$. The underlying field being HIT means the root-mean-square value for each of the three velocity components is this $u_{rms}$. Once the velocity and scalar fields are fully developed, $D^0_{ii} = -\overline {u_ic}$ is post-processed from the turbulent scalar flux for a time of $O(200\unicode{x2013}500)\tau _k$. Confidence intervals are calculated using the standard error of each mean statistic by constructing decorrelated samples of appropriate length compared with $\tau _k$ (Shirian, Horwitz & Mani 2023). For each simulation, $\varDelta /\eta = \varDelta /\lambda _B \lessapprox 2.1$, where $\varDelta$ is the grid spacing, $\eta \equiv \nu ^{3/4}\epsilon ^{-1/4}$ is the Kolmogorov length scale, $\lambda _B$ is the Batchelor scale and $\epsilon$ is the energy dissipation rate. We also report $Re_\lambda \equiv \sqrt {{15 u^4_{rms}}/{\epsilon \nu }}$. This resolution, which can also be reported as $k_{max}\eta \gtrapprox 1.5$, where $k_{max}$ is the maximum resolved flow wavenumber, is sufficient to resolve the low-order statistics studied herein and in line with canonical guidelines (Pope 2000), studies of flow energetics (Kaneda et al. 2003; Bassenne et al. 2016) and other works studying particle diffusivity (Squires & Eaton 1991; Mazzitelli & Lohse 2004; Jung, Yeo & Lee 2008).

3.1. Comparison of MFM with Lagrangian formulae

We show equivalence in the MFM and Lagrangian definitions of eddy diffusivity. The MFM measures $D^0$ by calculating the correlation of simulated scalar and velocity fields. The Eulerian–Lagrangian method (ELM), in contrast, simulates scalar particles in a background flow to calculate (2.1). Following Falkovich, Gawędzki & Vergassola (2001), the inertialess scalar particles are governed by ${\rm d}\boldsymbol {X} = (\boldsymbol {u}(t) + \boldsymbol {u_d})\,{\rm d}t + \sqrt {2D_m}\,{\rm d}\boldsymbol {W}$. Here, $\boldsymbol {W}(t)$ is a zero-mean, unity-variance Weiner process that allows this Langevin equation to correspond to the advection–diffusion equation solved by MFM.

To assess this equivalence, both methods measure eddy diffusivity from the same HIT case of $L_{box} = 2{\rm \pi}$ and $Re_\lambda = 14.4$ from table 1 with $u_d = 0$. In addition to the Eulerian DNS, MFM requires solution of (3.3) on $O(10^5)$ mesh points. The ELM requires the same flow-field DNS, along with Lagrangian trajectory simulations for $O(10^6)$ particles. Interpolation of velocity to particle locations uses modified Akima piecewise cubic Hermite functions.

For the candidate turbulent flow, figure 2 compares the measured estimate of the mean value of the eddy diffusivity as a function of the number of ensembles considered. Each independent ensemble is of length $\approx 10 \tau _k$, collected from a fully developed DNS. The mean estimate from both methods is approximately the same for this realization, but the confidence intervals are quite dissimilar in their size. The MFM provides more confident mean estimates for this application and is more scalable for simulations, as Eulerian fields are easier to distribute and solve with parallel computing than Lagrangian particles. An additional advantage of MFM not utilized here is that MFM can find higher non-local moments of the diffusivity kernel beyond the local-limit leading-order moment (Mani & Park 2021).

Figure 2. Mean estimates of $D^0$ for the $Re_\lambda =14.4$ case using MFM and ELM with 95 % confidence intervals showing convergence to the final estimate from each method (– –).

4. Results and discussion

Equation (2.8) is fitted to the mean diffusivity data from the cases of table 1 as a function of $u_d$. An iterative method is used to determine $\alpha _{11}$ and $\alpha _{22}$ and the forthcoming section will show that the value of $\alpha _{ii}$ is largely invariant to the tested Reynolds numbers.

In figure 3, $D^0_{11}$ and $D^0_{22}$ values are plotted for the $Re_\lambda = 14.4$ cases as a function of drift. The eddy diffusivity decay trends for the two directions differ for three values of $L_{box}$. At low drift velocities, the model accurately describes the data; however, at high drift, the $2{\rm \pi} ^3$ box data diverge from the analytic scaling. For a fixed box size, the computational value of $D^0_{ii}$ asymptotes to a non-zero value in the limit of large drift velocity, in violation of the analytical derivation presented earlier. This effect can be explained: in a finite periodic domain, a large drift velocity, $u_d > L_{box}/\tau _k$, subjects a scalar particle to see the same turbulence field in a characteristic advection time. Therefore, the autocorrelation becomes non-zero and the diffusivity value saturates. Figure 3 shows that increasing the computational domain size ameliorates this effect and improves data convergence to the model prediction for the highest drift velocities studied. This effect, however, has minimal effect on parameter choice: for further details, see Appendix B.

Figure 3. Convergence of $D^0_{11}$ and $D^0_{22}$ with respect to box size for $Re_\lambda = 14.4$. Dashed lines represent (2.8) for $\alpha _{11} = 3.9$ (a,b) and $\alpha _{22} = 1.9$ (c,d). Certain MFM data are shown with representative 95 % confidence intervals, which consider only the statistical error and exclude finite box-size error; arrows indicate convergence with increasing box size. Panels (a,c) show standard axes, while panels (b,d) show log–log axes. The data and Appendix B suggest criteria for which the box size of $2{\rm \pi}$ is sufficient.

We now examine the effect of $Re$. In figure 4, $D^0_{11}$ and $D^0_{22}$ at three $Re_\lambda$ values for $u_d/u_{rms} < 5$ are shown and the data collapse when normalized by $u^*_{ii}/u_{rms}$. There is good agreement between the computational data and model predictions over the drift velocities and Reynolds numbers explored for $\alpha _{11} \approx 4$ and $\alpha _{22} \approx 2$, based on the values from figure 3. Since $\alpha _{22} = 2$ is used in the literature and Appendix C suggests that $\alpha _{11} \leq 4$, we adopt the closest integer approximations of $\alpha _{11} = 4$ and $\alpha _{22} = 2$ for figure 4 and the remainder of this work. This observed agreement is better and has narrower scatter about the model prediction than presented in Squires & Eaton (1991), where $\alpha _{11}=\alpha _{22}=2$. This may be due to the long simulation time and Eulerian MFM method for quantification of $D^0_{ii}$.

Figure 4. The MFM-measured $D^0_{11}$ and $D^0_{22}$ as a function of $u_d/u^*_{ii}$ for different flow Reynolds numbers for $L_{box} = 2{\rm \pi}$. The dashed lines represent the analytical curve given by (2.8) for $\alpha _{11} = 4$ (top) and $\alpha _{22} = 2$ (bottom). Certain MFM data are shown with representative 95 % confidence intervals, which consider only the statistical error and exclude finite box-size error. Panel (a) shows standard axes, while panel (b) shows log–log axes.

Absent drift, Shirian & Mani (2022) used MFM to similarly show that scale-dependent eddy diffusivity, normalized by $u_{rms}$ and an eddy length scale, was largely invariant to $Re$. When considering $\overline {u_ic}$, this should not be surprising, as taking an average of a multi-scale field that decays rapidly in magnitude at large wavenumbers ensures that the large-scale content dominates. Bos, Touil & Bertoglio (2005), for example, show that even at $Re_\lambda = 28$, a relatively sharp scaling exponent means the large-wavenumber content of the passive scalar flux is far smaller in magnitude than that at the near-integral scales. As $Re$ increases, figure 1 shows us that $L_{ii}$ decreases: so $D^0_{ii}$ decreases with $Re$, ceteris paribus. This comment concerns the $Re_\lambda$ cases simulated here, but § 4.1 will discuss more general, higher Reynolds number conclusions. However, an increase in $Re$, properly normalized, should not greatly affect small-wavenumber quantities, which dominate the normalized measure of eddy diffusivity.

The value of $\alpha _{22} = 2$ corresponds to the transverse diffusivity results of Csanady (1963) and matches the overall scaling of Wang & Stock (1993), but the value of $\alpha _{11} = 4$ has not previously been reported. Squires & Eaton (1991) noted that their computational data and measurements of glass beads by Wells & Stock (1983) differed from the predictions of (2.10). We plot the corresponding data from Squires & Eaton (1991) for $D^0_{11}$ in figure 5 (cf. figure 11b of that work), and their version of (2.10), inferring $\beta = 1.1$ from their plot. If $\alpha _{11} = 4$ from this work is used, we better predict their presented data.

Figure 5. The present model compared with: (a) data from this work with parameter choices from table 1; from Wells & Stock (1983) and Squires & Eaton (1991) with $\beta = 1.1$; from He et al. (2005) with $\beta = 1.0$; and from Mazzitelli & Lohse (2004) with $\beta = 0.71$ with annotation for a particle motion case with and without lift. (b) Bounds at infinite $Re$ specified by Yudine (1959).

Squires & Eaton (1991) hypothesized that discrepancies between their data and Csanady's model might be due to assumptions about the form of velocity autocorrelation and a differing values of Lagrangian and spatial measures. The eddy length scale control and consistent use of Eulerian correlations mean the data in figures 3 and 4 are not affected by these considerations.

We can plot $D^0_{11}$ markers from figure 4 and observe that they fall along the same curve as the previous data, showing the lower scatter of MFM data about the model prediction. To build additional confidence, we can also plot computational bubble dispersion data from figure 2 of Mazzitelli & Lohse (2004) on the same axes of figure 5. For this one-way coupled bubble data, we approximate the $D^0$ value at zero drift to compute $u^*_{11}/u_{rms} = 1.4$. Most of the data come from bubble simulations where the particles experience a lift force. At the highest drift velocity, however, two values of diffusivity are provided, for which the lower value comes from a simulation where lift was neglected in the bubble dynamics. The data of figure 11(b) of He et al. (2005) are also included as another non-inertial point of comparison at $Re_\lambda = 42.4$. The markers of figure 5 come from a wide range of $Re$ values, $Re_ \lambda =[14-62]$, and yet all five datasets are better described by the model with $\alpha _{11} = 4$ for the examined range of swept drift velocities.

The choices of Reynolds numbers are broadly indicative of scenarios of $Re_\lambda = [38\unicode{x2013}60]$ that have been examined in the literature for target oceanic contexts (Shim et al. 2020). Specifically, for particle dispersion, Elghobashi & Truesdell (1992) found that their simulations of low-inertia particle mean-square displacement at even $Re_\lambda \leq 18$ agreed with corresponding experiments at $Re_\lambda \approx 48.5$.

So the value of $\alpha _{11} = 4$ is supported by data, but it also satisfies Csanady's arguments, namely that it corresponds to an autocorrelation that supports Taylor's frozen flow hypothesis. More fundamentally, Yudine (1959) calculated bounds for the value of the axial eddy viscosity as a function of drift based on the Kolmogorov–Obukhov structure function description of homogeneous turbulence. For the second-order structure function required, the refined Kolmogorov hypothesis does not significantly alter the conclusions (Pope 2000). These bounds are pictured in figure 5 for a $u^*_{11}$ imputed from that work and it is clear that $\alpha _{11} = 4$ is close to the envelope of realizable curves, which represent the limit state of infinite $Re$ turbulence. The value $\alpha _{11} \approx 1$ closely describes the corresponding lower bound. As the transverse structure function differs from the axial one by only a multiplicative constant, the bounding exponents are the same for the transverse diffusivity, so $\alpha _{22} = 2$ also satisfies theoretical constraints of Yudine (1959).

4.1. Generalizing to higher $Re$

Since the theoretical bounds are in the infinite $Re$ limit, the model in this work should be only a weak function of $Re$, as is suggested by figures 4 and 5.

To use (2.8) directly, a straightforward correlation of $u_{ii}^*$ with $u_{rms}$ would be useful. In the definition of $u_{ii}$, given by (2.9), $D^0$ could be written in terms of the $k$ and $\epsilon$ of the underlying turbulence; Shirian & Mani (2022) attempted to establish such a correlation and showed $D^0 = D_c u^4_{rms}/\epsilon$, where $D_c$ is an order-unity constant. Furthermore, Sreenivasan (1998) and Kaneda et al. (2003) offer evidence towards an ultimate regime of $\epsilon = \epsilon _c u_{rms}^3/L_{11}$, where $\epsilon _c$ is also a constant, allowing the establishment of

(4.1)\begin{equation} \frac{u^*_{11}}{u_{rms}} = \frac{u_{rms} L_{11}}{D^0} \approx \frac{\epsilon_c}{D_c} . \end{equation}

For $Re_\lambda = 67$, Shirian & Mani (2022) report $D_c = 0.77$ and we can estimate $\epsilon _c \approx 0.8\unicode{x2013}0.86$ (Sreenivasan 1998), which means (4.1) can be evaluated to $\approx 1.04\unicode{x2013}1.12$, in comparison with the value of $1.13$ obtained from the highest $Re_\lambda$ realized in this work. Beyond this empirical study of convergence, we can alternatively use equation (8) of Sawford, Yeung & Hackl (2008) as a model for the Lagrangian integral time scale to confirm a limiting case exists. Combining this expression with the correlation of Kaneda et al. (2003), we can write

(4.2) \begin{align} \frac{u^*_{11}}{u_{rms}} = \frac{u_{rms} L_{11}}{D^0} \approx \frac{\epsilon_c \dfrac{u_{rms}^4}{\epsilon}}{u^2_{rms} \sqrt{\dfrac{\nu}{\epsilon}} \Bigg( 4.77 + \left( \dfrac{Re_\lambda}{12.6} \right)^{{4}/{3}} \Bigg)^{{3}/{4}}} = \frac{0.258 \epsilon_c Re_\lambda}{\Bigg( 4.77 + \left( \dfrac{Re_\lambda}{12.6} \right)^{{4}/{3}} \Bigg)^{{3}/{4}}} . \end{align}

As $Re_\lambda$ increases, this expression approaches $\approx 3.25\epsilon _c$, bolstering support for our parameter reaching a constant value; for $Re_\lambda = 1201$, $\epsilon _c \approx 0.41$ such that (4.2) evaluates to $\approx 1.33$ (Kaneda et al. 2003). So while our data in table 1 and figure 5 suggest this ratio is an $O(1)$ constant, further investigations are needed to establish its convergence vis-à-vis that of $\epsilon _c$ with $Re$.

5. Conclusions

We propose an algebraic model for capturing the effect of drift velocity on the turbulent dispersion of passive scalars. This model, presented in (2.8), captures the asymptotic limits of zero and infinite drift velocities exactly and the single free parameter of $\alpha _{ii}$ captures the effects of intermediate drift velocities. By measuring $\alpha _{11} \approx 4$ and $\alpha _{22} \approx 2$ from the data and predicting eddy diffusivity across the span of drift velocities, this work performs a priori modelling without assuming a form for the underlying velocity structure. In so doing, we have shown a correspondence between MFM, an efficient Eulerian method for determining eddy diffusivity and other similar transport operators, and the classical Lagrangian definition. In particular, we show that MFM offers superior statistical convergence and computational properties for the efficient calculation of eddy diffusivity.

In the model form, $D^0$ and $u^*$ capture the $Re$ dependence, and can be measured independent of drift. As the $\alpha _{ii}$ that best describes the two directions are different, the transition between the limiting asymptotic behaviours appears to occur more rapidly in the transverse directions than in the axial direction, even accounting for the differences in eddy size.

There is cause to hope these results, which show improvements from previous models, can be applied to a wider range of situations than expected. As $St$ and $\phi$ increase, the velocity field seen by a particle no longer resembles the undisturbed Eulerian field and the fluid and particles become two-way coupled. Regarding the first parameter, Reeks (1977) and others have shown that an increase in Stokes number may increase diffusivity, but theory and evidence from simulations suggest this effect is relatively small for $St < 0.1$ (Jung et al. 2008). For example, Wang & Stock (1993) directly consider $St\neq 0$ and find diffusivity relations independent of $St$ outside of measurements at $u_d = 0$, viz. equations ((2.26)–(2.27)). With finite particle loading, Mathai et al. (2018) found axial diffusivity exceeded transverse diffusivity at $\phi$ as high as $5\times 10^{-4}$. Some of the particles of Wells & Stock (1983) had dynamics significant enough to affect the reported flow and Mazzitelli & Lohse (2004) simulated microbubbles with added mass, lift and drag, but they are still acceptably described by the model in figure 5. It may, therefore, be fair to conclude that the basic effects of anisotropic decreases in diffusivity with increasing drift obscure the effects of other complex nonlinear interactions. For transverse diffusivity, Loth (2023) shows that data with a range of negligible to finite Stokes numbers, such as from Groszmann, Fallon & Rogers (1999), can still be described by the Csanady (1963) and $\alpha _{22}$ models, albeit with deviations of the order of those seen in figure 5. For such cases, revisiting calculations to determine $D_{ii}^0$ and $L_{ii}$ in the presence of inertial particles might improve data collapse.

The predictive nature of this model, enabled by MFM, can be applied to other domains. The inverse problem can also be examined, in which dispersion measurements of natural particles or bubbles with known constant drift velocity can be used to diagnose turbulence parameters of the underlying flow. While this work's computational box set-up is conventionally used for Reynolds-averaged equation closure, this model can provide subgrid-scale closures in the large-eddy simulation context, as computational cells far from a wall should represent HIT.