Nonlocality of Mean Scalar Transport in Two-Dimensional Rayleigh-Taylor Instability Using the Macroscopic Forcing Method

In this work, a two-dimensional (2D), nonreacting flow with two species—a heavy fluid over a light fluid—is considered, with gravity pointing in the negative $y$-direction. It must be noted that the behavior of 2D RTI is significantly different from three-dimensional (3D) RTI, the latter of which is more relevant to problems of engineering interest. It is well known that while 2D RTI is unsteady and chaotic, it is not strictly turbulent, since turbulence is a characteristic of 3D flows. In addition, 2D RTI has a faster late-time growth rate, develops larger structures, and is ultimately less well-mixed. These differences have been studied in RTI by Cabot (2006) and Young et al. (2001) and in Richtmeyer-Meshkov instability by Olson & Greenough (2014).

For this study, 2D RTI is chosen as the model problem instead of 3D RTI, since it is a good simplified setting for understanding nonlocality in RTI through the lens of MFM. Specifically, 2D RTI simulations are much less computationally expensive than those of 3D RTI, and MFM requires many simulations to attain statistical convergence. Thus, 2D RTI remains the focus of this work, with the hope that the understanding of nonlocality in this flow could be extended to nonlocality in 3D RTI.

In this 2D problem, $x$ is the homogenous direction. In addition, there is no surface tension, the Atwood and Mach ($Ma$) numbers are finite but small, and the Peclet ($Pe$) number is finite but large.

In this work, the effect of nonlocality on mean scalar transport is of interest, so analysis begins with the scalar transport equation under the assumption of incompressibility:

where $\mathbf{u}$ is the velocity vector and $D_{H}$ is the molecular diffusivity of the heavy fluid.

After Reynolds decomposition and averaging, this becomes

In this work, large $Pe$ (the ratio of advective transport rate to diffusive transport rate) and small $A$ are assumed. The former assumption means molecular diffusion is negligible, and the latter yields $\langle u_{i}\rangle=0$, allowing the advective term to drop. Equation 11 becomes

The term $\langle v^{\prime}Y_{H}^{\prime}\rangle$ is the turbulent scalar flux, and this is the unclosed term that needs to be modeled.

As mentioned previously, one reason the gradient diffusion approximation used to model this term is inaccurate is that it assumes locality of the eddy diffusivity. This assumption can be removed by instead considering a generalized eddy diffusivity that is nonlocal in space and time, as demonstrated by Romanof (1985) and Kraichnan (1987). For 2D RTI, such a model reduces to

Here, $y$ is the spatial coordinate in averaged space and $t$ is the time at which the turbulent scalar flux is measured, $y^{\prime}$ is all points in averaged space, and $t^{\prime}$ is all points in time. This definition is exact for passive scalar transport, including in the case studied in this work.

This nonlocal eddy diffusivity can also be viewed as a two-point correlation. This was first described by Taylor (1922) in homogeneous turbulence. Through Lagrangian statistical analysis, Taylor derived the following relation between diffusivity and velocity correlations:

Work by Shende et al. (2023) has shown that MFM recovers this Lagrangian formulation for eddy diffusivity in homogeneous flows. It should be noted that the above definition is not valid for inhomogeneous RTI (again, the exact definition of eddy diffusivity for the studied flow is the one in equation 13), but the intent here is to provide another interpretation of MFM more aligned with the well-understood two-point correlations.

The eddy diffusivity kernel can be approximated by Taylor-series-expanding the scalar gradient locally about $y$ and $t$, which results in the following Kramers-Moyal-like expansion for the turbulent scalar flux as done by Kraichnan (1987), Hamba (1995), and Hamba (2004):

Here, $D^{mn}$ are the eddy diffusivity moments; the first index, $m$, denotes order in space, while the second, $n$, denotes order in time. This is the form presented in Mani & Park (2021) and Liu et al. (2023).

When the eddy diffusivity kernel is purely local,

In this case, $D^{00}$ is the only surviving moment, while all higher-order moments in space and time are zero. Any non-zero higher-order moment therefore characterizes the nonlocality of the eddy diffusivity kernel. Thus, this expansion implies explicitly a model form for the turbulent scalar flux that incorporates nonlocality of the eddy diffusivity. Truncating the expansion provides an approximation of $\langle v^{\prime}Y_{H}^{\prime}\rangle$ but with the caveat that the expansion may not converge. This will be discussed in more detail later in §5.3.1.

Each $D^{mn}$ provides more information about the eddy diffusivity kernel with increasing order. For example, $D^{00}$ represents the volume of the kernel in space-time. The coefficient corresponding to one higher-order in space, $D^{10}$, provides information about the centroid of the kernel in space. $D^{20}$ contains information about the moment of inertia of the kernel in space, $D^{01}$ contains information about the centroid of the kernel in time, and so on.

MFM is a method for numerically determining closure operators in turbulent flows (Mani & Park, 2021). Much like a rheometer measures the molecular viscosity of a fluid by imposing a shear force on the flow, MFM forces the transport equation in a turbulent flow and extracts the closure operator from its response. Unlike the molecular viscosity, which is a material property, the turbulent closure operator is a property of the flow, so MFM measurements of one flow cannot be generalized for all flows; the MFM-measured closure for one flow cannot be applied exactly as it is to a different flow. However, MFM measurements of one flow can reveal characteristics of the turbulent closure that are expected be true for a family of similar flows.

Specifically, MFM can be used to determine the RANS closure operator, as shown in the pipeline diagram in figure 1. In MFM, two simulations are run at once: the donor and the receiver simulation. In this work, the donor simulation numerically solves the multicomponent Navier-Stokes equations in equations 38 - 41. The receiver simulation “receives” $u_{i}$ from the donor simulation and uses it to solve the scalar transport equation with a forcing $s$:

Ultimately, forcings on the receiver simulation effect a response from the flow, and measuring this response allows for determination of the eddy diffusivity. In particular, these forcings are macroscopic. Here, macroscopic quantities are defined as fields that are unchanged by Reynolds-averaging. Mathematically, the macroscopic forcing is such that $s=\overline{s}$. This macroscopic nature is crucial to the method, since it does not disturb the underlying mixing process, which allows for measurement of the closure operator without changing it. For details, see Mani & Park (2021).

In actuality, the Inverse Macroscopic Forcing Method (IMFM) is used to determine eddy diffusivity moments. That is, instead of the forcings being chosen, certain mean mass fraction fields are chosen. Numerically, mean mass fractions are enforced in each realization, so the averages (denoted by $\overline{*}$) described here are in $x$, the homogeneous direction in space. The forcing needed to maintain the chosen $\overline{Y_{H}}$ is determined implicitly along the process and is not directly used in the analysis.

As an illustration, the measurement of $D^{00}$ can be considered. According to equation 15, choosing $\overline{Y_{H}}=y$ (for $y$ between $-1/2$ and $1/2$) results in $\frac{\partial\overline{Y_{H}}}{\partial y}=1$, and all other higher-order derivatives are zero. Thus, choosing this $\overline{Y_{H}}$ in each realization results in the realization- and spatially-averaged measurement $-\langle v^{\prime}Y_{H}^{\prime}\rangle=D^{00}$.

Measurement of higher-order moments involves similar choices of $\overline{Y_{H}}$ but requires information from lower order moments. For example, measuring $D^{10}$ involves choosing $\overline{Y_{H}}=y^{2}$, which results in $-\langle v^{\prime}Y_{H}^{\prime}\rangle=yD^{00}+D^{10}$. Here, $D^{00}$ comes from the simulation using $\overline{Y_{H}}=y$. Thus, $D^{10}$ is computed by subtracting $yD^{00}$ from the $\langle v^{\prime}Y_{H}^{\prime}\rangle$ measurement from the simulation using $\overline{Y_{H}}=y^{2}$.

Specifically, the following desired mean mass fractions are used for each moment for $y$ between $-1/2$ and $1/2$:

From these $\overline{Y_{H}}$, the needed forcing in each timestep is numerically determined:

where the superscript $k$ denotes the timestep number, $\overline{Y_{H}}_{\text{desired}}$ is the mean mass fraction desired as outlined in equations 22 - 25, and $\Delta t$ is the timestep size.

This MFM forcing bears some resemblance to other forcings used in the literature, such as interaction by exchange with the mean (IEM) (Sawford, 2004; Pope, 2001). One main difference between forcings in such methods and MFM is that the purpose of the latter is to drive the flow to a specified mean gradient, which allows for measurement–not enforcement –of the eddy diffusivity moments. In other words, in MFM for scalar transport, the input is a mean scalar gradient, and the output is the eddy diffusivity moment; in IEM and similar methods, the input is a desired moment (e.g., in IEM, the input moment is $\langle c^{2}\rangle$) and the output is a mixing model. In addition, methods such as IEM use microscopic forcings, while MFM uses macroscopic forcings, which is a distinguishing characteristic of the latter method.

To determine $D^{00}$, $D^{10}$, $D^{01}$, and $D^{20}$, four separate simulations are needed. For each of these simulations, the moments can be calculated using measurements of the turbulent scalar flux as follows:

where $F^{mn}$ denotes the $-\langle v^{\prime}Y_{H}^{\prime}\rangle$ measured from the receiver simulation using the forcing corresponding the moment $D^{mn}$.

We perform our analysis in the self-similar regime. First, we define a self-similar coordinate:

so that $\langle Y_{H}\rangle$ is only a function of $\eta$. Note that $\eta$ requires a definition of $h(t)$. From the previous discussion on the self-similarity of RTI, an appropriate definition is $h(t)=\alpha Agt^{2}$.

Through self-similar analysis of equation 15, the eddy diffusivity moments and turbulent scalar flux can be normalized. Details of this process can be found in the Appendix.

Recall that $h(t)=\alpha Agt^{2}$ is used in the self-similarity analysis. This is valid only for late time, so the subsequent analyses in this work are all done in this self-similar timeframe. Usually, $\alpha$ can be determined from $\frac{h(t)}{Agt^{2}}$, where $h(t)$ is computed from the simulation via equation 5. However, due to the convergence and statistical errors as well as the existence of a virtual time origin, $\alpha Agt^{2}$ is not a good representation of $h(t)$ measured in the DNS. Instead, a fitting coefficient $\alpha^{*}$ and virtual time origin $t^{*}$ are determined to make a shifted quadratic fit to $h(t)$ from the simulation:

With this fit, the normalizations of the turbulent scalar flux and moments become

For exact self-similarity, plots of the measured $\widehat{D^{mn}}$ against $\eta$ must be independent of time. This expectation sets a criterion to assess the extent to which ideal self-similarity is achieved. Plots and assessment of the self-similar collapse of the measurements presented in this work are in the Appendix.

The governing equations solved in this work are the compressible multicomponent Navier-Stokes equations, which involve equations for continuity, diffusion of mass fraction $Y_{\alpha}$ of species $\alpha$ (characterized by its binary molecular diffusivity $D_{\alpha}$), momentum transport, and transport of specific internal energy $e$:

Here, $\frac{D}{Dt}$ is the material derivative $\frac{\partial}{\partial t}+u_{i}\frac{\partial}{\partial x_{i}}$, $\rho$ is density, $u$ is velocity, $p$ is pressure, and $g$ is gravitational acceleration, active in the $-y$ direction. The viscous stress tensor $\sigma_{ij}$ and heat flux vector $q_{j}$ are respectively defined as

Here, $\mu$ is the dynamic viscosity, $\kappa$ is the thermal conductivity, $T$ is temperature, and $h_{\alpha}$ is the specific enthalpy of species $\alpha$.

Component pressures and temperatures of each species are determined using ideal gas equations of state. Under the assumption of pressure and temperature equilibrium, an iterative process is performed to determine volume fractions $v_{\alpha}$ that allow for computation of partial densities and energies. More details on the hydrodynamics equations and computation of component quantities can be found in Morgan et al. (2018).

Finally, total pressure is determined as the weighted sum of component pressures:

In general, in these compressible equations, $Y_{\alpha}$ are not passive scalars. However, the component equations of state are scaled so that a consistent hydrostatic pressure gradient is maintained across the mixing layer. Thus, in this work, $Y_{\alpha}$ are effectively passive.

Simulations for 2D RTI are run using the Ares code, a hydrodynamics solver developed at Lawrence Livermore National Laboratory (LLNL) (Morgan & Greenough, 2015; Bender et al., 2021). Ares employs an arbitrary Lagrangian-Eulerian (ALE) method based on the one by Sharp & Barton (1981), in which the governing equations (equations 38 to 41) are solved in a Lagrangian frame and then remapped to an Eulerian mesh through a second-order scheme. The spatial discretization is a second-order non-dissipative finite element method, and time advancement is a second-order explicit predictor-corrector scheme.

The Reynolds number (more specifically, the kinematic viscosity $\nu$) is set through a numerical Grashof number, such that

Here, $\Delta$ is the grid spacing; in the simulations, a uniform mesh is used, and $\Delta=\Delta x=\Delta y$. To ensure that the unsteady structures are properly resolved and for the simulation to appropriately be considered a DNS, $Gr$ should be kept small. A $Gr$ that is too large results in a simulation with dissipation dominated by numerics rather than the physics. Morgan & Black (2020) found that past $Gr\approx 12$ in the Ares code, numerical diffusivity dominates molecular diffusivity. For our simulations, we use $A=0.05$ and $Gr=1$, the latter of which is in line with the DNS by Cabot & Cook (2006). These choices give a $\nu$ of $10^{-9}m^{2}/s$. The Schmidt number $Sc$, defined as $\nu/D_{M}$, is set to unity, so $D_{M}=10^{-9}m^{2}/s$.

The Mach number, $Ma=\frac{u}{c}$, where $c$ is the speed of sound, characterizes compressibility effects of the flow. $Ma$ is set by the specific heat ratio $\gamma$, which is $5/3$ in the simulations im this work. The maximum $Ma$ is measured at the last timestep to be approximately $0.03$, which is ascertained to be small enough to assume incompressibility.

The Peclet number $Pe$ characterizes the advection versus diffusion rate and is defined as $ReSc$. Here, a $Pe_{L}$ and a $Pe_{T}$ are reported, which use a large-scale $Re_{L}$ and the Taylor Reynolds number $Re_{T}$, respectively. In the presented simulations, $Sc=1$. The two $Pe$ are computed in post-processing: $Pe_{L}$ is approximately $8,000$, and $Pe_{T}$ is approximately $54$. Both are below the criterion established by Dimotakis (2000), suggesting that the simulated flow is transitional or pre-transitional.

The number of cells in each simulation is $2049\times 2049$. The width $L_{x}$ of the domain is $1$, and the height $L_{y}$ is $1$. The boundary conditions are periodic in $x$ and no slip and no penetration in $y$.

Initially, the velocity field is zero, temperature is 293.15 K, and pressure is 1 atm. A tophat perturbation based on the ones used by Morgan & Greenough (2015) and Morgan (2022) is imposed on the density field at the interface of the heavy and light fluids:

where $\phi_{1,k}$ and $\phi_{2,k}$ are phase shift vectors randomly taken from a uniform distribution, and $L_{y}$ is the length in $y$ of the domain. Here, the minimum and maximum wavenumbers are set to $\kappa_{\text{min}}=8$ and $\kappa_{\text{max}}=256$, respectively.

The stop condition of the simulations is when $h$ is approximately half the domain size in $y$. This corresponds to the nondimensional simulation time $\tau$ of $30.84$. $\tau$ is defined as $\frac{t}{t_{0}}$, where $t_{0}=\sqrt{\frac{h_{0}}{Ag}}$ and $h_{0}$ is the dominant lengthscale determined by the peak of the initial perturbation spectrum.

Before the MFM analysis was conducted, the results of the donor simulations were examined. In figure 2(a), mixedness is observed to reach a value of around $0.6$ but appears to not have converged yet. Figure 2(b) shows the two definitions of $\alpha$ over time. The first definition, $\alpha=h/Agt^{2}$, reaches a value of about $0.05$ by the end of the simulation, but it does not appear to be converged. The second definition, $\alpha=\dot{h}^{2}/4Agh$, is oscillatory, due to the sensitivity of the time derivative to noise, and it appears to fluctuate about a value of approximately $0.04$. It is acknowledged that this behavior indicates that the RTI simulated here is only weakly turbulent. However, it is observed that the flow is still self-similar at late times. The contour plot of $\langle Y_{H}\rangle$ in figure 3 exhibits parallel contour lines after $\tau\approx 17$, indicating self-similarity at those times. It is also shown in the Appendix that the mean concentration and normalized turbulent scalar flux profiles exhibit self-similar collapse after $\tau\approx 17$, so the presented self-similar analysis is valid.

Figure 4 shows a plot of the algebraic fit for $h$, described in equation 32. For the simulations presented here, $\alpha^{*}$ is $0.0046$ and $t^{*}$ is $-1.6\times 10^{3}$. The plot shows a strong quadratic dependence of $h$ on $t$ at late time, as $h_{\text{fit}}$ matches DNS at $\tau\gtrapprox 17$, validating the self-similar ansatz of $h\sim t^{2}$.

To further ensure the simulations are working as desired, the flow fields of the donor and receiver simulations can be examined qualitatively. The $Y_{H}$ contours at the last timesteps of each simulation are shown in figure 5. The receiver simulation shown is the one used to compute $D^{00}$ (where $\langle Y_{H}\rangle=y$). Self-similar RTI turbulent mixing is observed at this timestep, where the characteristic heavy spikes are sinking into the lighter fluid and the light bubbles rise into the heavier fluid. Both simulations have the same velocity fields, since the receiver simulation “receives” the velocity field from the donor simulation. In contrast with the donor simulation, which has a stark black-and-white difference between the heavy and light fluids, there is a grey gradient of density in the receiver simulation due to the imposed mean scalar gradient. The fluctuations of $Y_{H}$ in each simualtion are also compared in figure 6. The $Y_{H}^{\prime}$ contours are not identical but are qualitatively very similar. In both simulations, $Y_{H}^{\prime}$ is constrained to the mixing layer. Based on these observations, it is concluded that the simulations are visually working as intended.

Figure 7 shows normalized MFM measurements of the eddy diffusivity moments $D^{00}$, $D^{10}$, $D^{01}$, and $D^{20}$ averaged over $1,000$ realizations and the homogeneous direction $x$. Some expected characteristics of the measured moments are observed:

1. 1.

   The leading order moment is over two magnitudes larger than the molecular diffusivity ($10^{-9}m^{2}/s$). The scaled higher order moments shown are all at least one magnitude larger than the molecular diffusivity.

2. 2.

   $D^{00}$ is symmetric and always positive.

3. 3.

   $D^{10}$ is antisymmetric. This antisymmetry can be understood by interpreting $D^{10}/D^{00}$ as the centroid of the eddy diffusivity kernel. Physically, for $\eta>0$, it is expected that the mean scalar gradient at the center of the mixing layer (at a negative distance away) has more influence on the turbulent scalar flux than the mean scalar gradient at the outer edges, since the mixing layer is growing outwards. This makes the eddy diffusivity kernel skewed more towards the center of the domain, so $D^{10}<0$ for $\eta>0$. A similar effect occurs for $\eta<0$, which results in $D^{10}>0$.

4. 4.

   $D^{01}$ is symmetric and always negative. The latter must be true for the flow to depend on its history (it does not violate causality).

5. 5.

   $D^{20}$ is symmetric and always positive, as is characteristic of moment of inertia of a positive kernel.

Based on the magnitudes of the normalized moments, some initial observations on importance of each moment can be made. $D^{00}$ has the highest magnitude of all the moments, which is expected since it is the leading-order moment. The magnitudes of $D^{10}/h$ and $D^{01}/t$ are on the order of $10\%$ of the magnitude of $D^{00}$, which suggests that they are non-negligible. On the other hand, the magnitude $D^{20}/h^{2}$ is on the order of $1\%$ that of $D^{00}$, so $D^{20}$ is likely not an important moment to include in modeling RTI. More robust studies will be presented in the following sections to determine importance of each of the eddy diffusivity moments.

It is also observed that there is statistical error in the measurements. Due the chaotic nature of RTI, the moment measurements contain statistical error, but this error can be reduced by averaging many realizations. To demonstrate statistical convergence of the measurements, plots of $D^{00}$ averaged over different numbers of realizations are included in figure 8. As the number of realizations increases, the plots become smoother, and it is found that after $1,000$ realizations, the rate of reduction in statistical error slows down significantly. Averaging over this number of realizations results in a smooth $D^{00}$ measurement and higher-order moment measurements with acceptably less statistical error.

Additionally, the higher the order of the moment, the slower its rate of statistical convergence. Recall that determination of higher-order moments requires information from lower-order moments. For example, in determining $D^{01}$, $tD^{00}$ is subtracted from $F^{01}$, the turbulent scalar flux measurement in the simulation associated with $D^{01}$. Naturally, there is statistical error associated with both $D^{01}$ and $D^{00}$. However, the error in $D^{00}$ is amplified by $t$, so the overall statistical error of $D^{01}$ increases with time. This statistical error “leakage” occurs for all higher-order moments. The higher the order of the moment, the worse the statistical error, since information from more lower-order moments is needed, and so more statistical error is accumulated and amplified. The relatively high statistical error of the higher-order moments makes it challenging to study their importance. Particularly, taking derivatives of quantities with high statistical error amplifies the error, so smoother measurements are desired. In this work, the moment measurements are smoothed using a Savitzky-Golay filter function in Matlab with a polynomial order of unity and window size of $191$. These smoothed moments are shown in figure 9. While it is possible to design an alternative formulation of MFM that removes leakage of statistical error from low-order moments to higher-order moments (see Lavacot et al., 2022), for this 2D study and for the order of moments considered here, the statistical convergence is sufficient.

Using these measurements, nonlocal timescales and lengthscales ($t_{NL}$ and $L_{NL}$, respectively) can be defined:

Note that this analysis can only be done for $-1\leq\eta\leq 1$, since the moments are analytically zero outside the mixing layer.

Nondimensionally, the nonlocal timescale is $\tau_{NL}=t_{NL}/t_{0}$, and the nonlocal lengthscale is $\eta_{NL}=L_{NL}/h$. Contour plots of the nondimensionalized nonlocal timescale and lengthscale are in figure 10. Note that $\tau_{NL}$ scales as $\tau$, so profiles of $\tau_{NL}/\tau$ against $\eta$ are also plotted in figure 11 in the self-similar time regime ($\tau>17$). The scaled profiles collapse and have a centerline value of approcimately $0.1$. This means that the mean fluxes at some time $\tau$ are affected by mean scalar gradients $0.1\tau$ earlier. Figure 12 shows the minimum nonlocal lengthscale is at the centerline, where $\eta_{NL}\approx 0.09$. The maximum lengthscales occur near the outer edges of the mixing layer: at around $\eta=\pm 0.87$, $\eta_{NL}\approx 0.27$. This indicates that mean fluxes at the mixing layer edges depend mostly on mean scalar gradients about a quarter of a mixing width away, while mean fluxes at the centerline depend on mean scalar gradients about one tenth of a mixing width away; nonlocality appears to be stronger at the mixing layer edges than at the centerline. These nonlocal properties of the eddy diffusivity for RTI could not be predicted without direct measurement of the eddy diffusivity moments, which has been made possible through MFM.

In this section, two RANS operator forms using information about the nonlocality of the eddy diffusivity are presented. These are the explicit and implicit operator forms; the former is a truncation of the turbulent scalar flux expansion 15, and the latter will be presented shortly. It must be stressed that the intention of the following studies is not to propose a new RANS model. Ultimately, a RANS model should not depend on direct MFM measurements that can only be retrieved from impractically many DNS. Instead, these studies are performed to further assess the importance of each of the eddy diffusivity moments, determine which combinations of moments best enhance the performance of a RANS model, and examine the differences between the explicit and implicit operator forms. The aim of these studies is to inform development of more predictive RANS models for RTI, not to suggest that these are the exact models that should be used.

In addition, $D^{20}$ will not be included in the following studies. This is mainly due to the high statistical error in the measurement that makes it difficult to ascertain whether errors in the results are due to this statistical error or solely the addition of the moment to the model. From the comparison of terms in §5.2.1, it is expected that $D^{20}$ is not as important as $D^{10}$ and $D^{01}$ to include in a RANS model. This should be tested in future work when a more statistically-converged measurement is achieved for $D^{20}$, ideally in a 3D analysis.