The mean momentum and heavy mass fraction, turbulent kinetic energy, and heavy mass fraction variance fields, as well as the budgets of their transport equations are examined several times during the evolution of a narrowband Richtmyer-Meshkov instability initiated by a Mach 1.84 shock traversing a perturbed interface separating gases with a density ratio of 3. The results are computed using the “quarter scale” data from four algorithms presented in the *θ*-group study of Thornber *et al.* [“Late-time growth rate, mixing, and anisotropy in the multimode narrowband Richtmyer-Meshkov instability: The θ-group collaboration,” Phys. Fluids **29**, 105107 (2017)]. The present study is inspired by a previous similar study of Rayleigh-Taylor instability and mixing using direct numerical simulation data by Schilling and Mueschke [“Analysis of turbulent transport and mixing in transitional Rayleigh-Taylor unstable flow using direct numerical simulation data,” Phys. Fluids **22**, 105102 (2010)]. In addition to comparing the predictions of the data from four implicit large-eddy simulation codes, the budgets are used to quantify the relative importance of the terms in the transport equations, and the balance of the terms is employed to infer the numerical dissipation. Terms arising from the compressibility of the flow are examined, in particular the pressure-dilatation. The results are useful for validation of large-eddy simulation and Reynolds-averaged modeling of Richtmyer-Meshkov instability.

## I. INTRODUCTION

Richtmyer-Meshkov instability (RMI) occurs when a perturbed interface between two fluids of differing properties is impulsively accelerated.^{1,2} This acceleration imparts a vorticity field near the interface, causing any perturbations on the interface to grow in time, regardless of the direction of the impulse relative to the layer. For initial perturbations with amplitude over wavelength *a*/*λ* ≪ 1, the growth rate of the instability is at first linear.^{1} As the perturbation amplitude becomes nonlinear, the perturbation growth rate decreases and the shear layers between the peaks and troughs roll up to form vortices. Shear along the interface also triggers smaller scale Kelvin-Helmholtz instabilities, which at later times contribute to the development of a range of vortical structures of differing sizes. For multimode perturbations, this process occurs simultaneously for a range of length-scales, leading to the development of an inhomogeneous turbulent mixing layer at late times. As the deposition of vorticity is impulsive, the layer with linear perturbations evolves through several growth regimes: (i) linear, (ii) nonlinear, (iii) nonlinear transitional (range of length-scales developing), and (iv) decaying inhomogeneous variable density turbulent layer.

This instability occurs in applications ranging from inertial confinement fusion^{3} and astrophysical flows^{4} to augmented mixing in scramjets and jet exhaust plumes.^{5} For a comprehensive recent review of research in this field, see Refs. 6 and 7. RMI differs significantly from Kelvin-Helmholtz instability (KHI) and closely related Rayleigh-Taylor instability (RTI) in that the driving is impulsive. With a sufficiently strong impulse, time to develop, and high Reynolds number, linear growth will be followed by a transition to turbulence, where velocity fluctuations decay in time and the layer grows at a relatively slow rate ∝ *t*^{θ} with *θ* ≈ 0.2–1 (possible if dominated by linear modes, although a fully turbulent layer is limited to^{8} *θ* ≤ 2/3) with a strong dependence on initial conditions.^{9–15} As an example *θ* ≈ 0.29 for the configuration studied in this paper,^{16} compared to KHI (∝*t*) and RTI (∝*t*^{2}). Despite the relatively small temporal exponent, the growth rate of RMI is sufficiently high that it significantly impacts the physics of the aforementioned applications, so there is substantial interest in the accurate modeling of turbulent transport and mixing in RMI.

In the context of current computational power, turbulent transport models for unsteady Reynolds-Averaged Navier-Stokes (URANS) or Large-Eddy Simulation (LES) are particularly important. URANS methods are in active development and use for design problems, principally due to the complexity of applied computations typically impacted by multiple physical phenomena, complex equations of state, and extreme conditions. Turbulence resolving computations such as Large-Eddy Simulations require resolution of the full dimensionality of the problem,^{3} whereas for rapid design evaluation and iteration, it is desirable to work with a problem reduced to one or two dimensions by appropriate averaging and application of modeling of the full three-dimensional effects. This is the key motivation for the development and use of URANS modeling.

Modeling approaches for decaying, inhomogeneous, and variable density turbulent flows where density variations do not satisfy the Boussinesq approximation are particularly complex. The majority of turbulence modeling applied to variable density flows employs a density-weighted or Favre-averaged approach.^{17} This simplifies the resulting equations describing the evolution of the density-weighted mean compared to Reynolds-averaging. Following the Favre decomposition, two principal additional complexities arise (i) additional terms associated with density fluctuations and (ii) a change in the definition of the “mean” evolved, perhaps necessitating a modification to the closures of terms from the Reynolds averaged formulation.^{18} URANS modeling approaches applicable to RMI range from single-fluid two-, three-, and four-equation^{19–24} and Reynolds stress^{25–28} to two fluid models.^{29–31} There is a strong need for data to validate and further advance these and other models.

Despite significant advances in facilities and diagnostics,^{11,32–40} experimental data are very challenging to acquire due to the relatively short time scale of RMI in terrestrial facilities, extreme conditions, or uncertainty in the initial perturbations. Such experimental research is critical and permits important validation of individual terms; however, to date, it is insufficient to fully determine the required model coefficients for URANS computations. While there is considerable literature examining turbulent transport in shear layers induced by KHI, and some prior studies of compressible RTI,^{41} to date, there has not been an equivalent study of all individual terms in the turbulent transport equations computed simultaneously for a turbulent mixing layer developing from a RMI. Thus, the objective of the present study is to address a gap in the identification of the key mechanisms responsible for variable density effects in the RMI and the terms that must be modeled to accurately represent mean transport.

Recently, a benchmark multimode, three-dimensional RMI case was computed with multiple independent computational fluid dynamics codes, establishing a well understood “*θ*-group” dataset.^{16} The principal goal of that study was to understand and quantify numerical uncertainty in the computation of such flows and to provide a foundation for basic quantities of interest. The present study computes the transport equation budgets for the mean momentum, mean heavy fluid mass fraction, heavy fluid mass fraction variance, and specific turbulent kinetic energy from the results of the four codes presented in that paper: Triclade, Turmoil, Flamenco, and Flash. The analysis of these data identifies the terms of greatest importance in the mean transport equations and provides a reliable benchmark for LES and URANS model development and validation.

This paper is organized as follows: Sec. II summarizes the numerical methods employed to generate the *θ*-group dataset, the test case examined, the transport equations examined, and numerical methods employed to compute the terms. Section III presents and analyses the turbulent transport budgets for four representative times ranging from just postshock to a weakly turbulent state. Finally, the key conclusions are summarized in Sec. IV.

## II. GOVERNING EQUATIONS AND COMPUTATIONAL DETAILS

### A. Summary of the numerical data employed

This paper utilizes the “quarter scale” data from four algorithms presented in the *θ*-group study of Thornber *et al.*^{16} In that paper, the full details of the initial conditions, results, and grid convergence studies for each individual algorithm are presented. It established an initialization to explore the physics of narrowband Richtmyer-Meshkov instability formulated in a way, which could be simulated by a wide range of algorithms. The configuration consisted of a light and a heavy gas, with a density ratio of 3, separated by a perturbed interface. A shock wave of strength Mach 1.84 impinges on the perturbed layer, triggering Richtmyer-Meshkov instability. The perturbation on the interface is defined using a narrowband power spectrum with constant amplitude between a defined minimum and maximum wavenumber, where modal amplitudes and phases are defined using deterministic random numbers, thus generating an exactly reproducible initial condition.

The initial conditions in the shocked heavy fluid domain 0.0 m < *x* < 3.0 m are (*ρ*, *u*, *p*) = (6.375 kg/m^{3}, −61.49 m/s, 400 kPa). Given a surface perturbation *A*(*y*, *z*), the unshocked heavy fluid lies between 3.0 m < *x* < 3.5 m + *A*(*y*, *z*) with initial properties (*ρ*, *u*, *p*) = (3.0 kg/m^{3}, −291.6 m/s, 100 kPa). The third region contains the unshocked light fluid that lies between 3.5 m + *A*(*y*, *z*) < *x* < 2.8*π* m and is initialized with properties (*ρ*, *u*, *p*) = (1.0 kg/m^{3}, −291.6 m/s, 100 kPa). The configuration is shown schematically in Fig. 1.

The perturbation power spectrum is constant between length scales of *L*/32 (*λ*_{min}) and *L*/16, where *L* is the cross section, and an initial diffuse interface of error function form and thickness *δ* = *L*/128 is applied as follows:

where *f*_{1} is the volume fraction of the heavy gas and *S*(*y*, *z*) = 3.5 + *A*(*y*, *z*). The standard deviation of the perturbation amplitude is *σ*_{0} = 0.1*λ*_{min} which places the shock-interface interaction initially at the large-amplitude end of the linear regime at the highest wavenumber. The diffuse initial condition ensures that the initial condition was accessible for all numerical schemes, as some algorithms are unstable in the presence of sharp interfaces. Together, the addition of a diffuse interface and the larger amplitude will impact the initial growth rates, particularly for the highest wavenumbers, which could be expected to be ≈20% lower than that predicted from linear theory.^{42} Although diffusion is not explicitly represented in the computation, an implicit assumption is that the gases are miscible. Turbulent layers evolving in immiscible fluids may have different growth rates and mixing parameters (see, for example, Refs. 43–45).

The heavy and light unshocked gases are at the same temperature in the undisturbed flow, and the specific heats for each component (*C*_{v1} and *C*_{v2}) are chosen to ensure this. Finally, the problem description is completed by assuming an ideal gas equation of state with a ratio of specific heats *γ* = 5/3.

The Cartesian computational domain is *L*_{x} × *L*_{y} × *L*_{z} = 2.8*π* × 2*π* × 2*π* m^{3}. The simulations were run from *t* = 0 to 5 s, which is equivalent to a dimensionless time *τ* = 246 (the nondimensionalization is detailed in Sec. II C). However, at times later than 2 s, individual spikes were observed to exit the domain; thus, only data up to 2 s are typically employed. Prior analysis has indicated that by the final time an asymptotic approximately self-similar state is approached—this assumption is examined in more detail here. Four time instants are examined, where the first time is immediately following shock passage (0.01 s), the second is shortly after this time where layer mixedness as measured by the mix parameter Θ has approximately doubled (0.025 s), and the third and fourth are approaching a self-similar state (1.0 and 2.0 s). See the discussion in Sec. III A and Fig. 4 for more details.

The configuration was chosen such that the low order statistics for the problem are just converged on the finest grid level using modern compressible algorithms. This introduces a degree of uncertainty on the influence of numerical dissipation on each transport term. This has been addressed in two ways: the first is through a grid convergence study and the second is by computing all transport terms using data produced from four independent codes/algorithms: Flamenco (Godunov method with the fifth order low-Mach-number-corrected MUSCL scheme^{46,47}), Flash (compressible Euler with the Piecewise Parabolic Method and Monotonized Central limiter^{48–50}), Triclade (conservative finite difference using the wave propagation algorithm of LeVeque^{51} with high order accurate corrections^{52}), and Turmoil (Lagrange remap method with high-order artificial viscosity^{53–55}). The boundary conditions are periodic in the homogeneous directions, but the best practice is adopted for each code in the shock-direction, which permits the transmission of the reflected rarefaction and transmitted shock through the boundaries with minimal reflections. The implicit numerical dissipation mechanisms in each of the algorithms differ substantially. Thus, agreement between the results either implies that for a given quantity (i) numerical dissipation has been minimized or (ii) numerical dissipation has acted to dissipate the quantity in a statistically similar manner. All simulations have no model for physical viscosity, diffusivity, or conduction.

### B. Mean and fluctuating components

This study explores the individual transport terms in the mean momentum and mass fraction, mass fraction variance, and fluctuating kinetic energy equations. A detailed prior study of Rayleigh-Taylor mixing of Schilling and Mueschke^{41} presented the required transport equations which are considered here. Referring to the equation numbering in Schilling and Mueschke,^{41} this paper examines the transport of mean momentum, Eq. (6), mean mass fraction, Eqs. (7) and (8), turbulent kinetic energy, Eqs. (14)–(15e), and mass fraction variance, Eqs. (19)–(20c). Averaging processes are defined in Eqs. (2) and (3) and accompanying text. In this contribution, the key differences are that the homogeneous directions are *y* and *z*, and the shock-direction is *x*, with a velocity component *u*, different from the Rayleigh-Taylor case where the principal direction and velocity are labeled *z* and *w*, respectively.

The Reynolds average of a field $\varphi x , t $ over the statistically homogeneous plane is

which enables a decomposition into mean and fluctuating components, where $\varphi ( x , t ) = \varphi \xaf ( x , t ) + \varphi \u2032 ( x , t ) $ and the fluctuating component is denoted by a prime. The Favre average is

where *ρ*(**x**, *t*) is the density and $ \varphi \u0303 ( x , t ) $ is defined such that $\varphi ( x , t ) = \varphi \u0303 ( x , t ) + \varphi \u2033 ( x , t ) $. This averaging reduces the fields to being one-dimensional functions of *x* only.

The individual terms which may be extracted from a single ILES are listed below, along with the notation used in the figures and analysis. For mean momentum transport,

where $ \tau 11 = \rho \xaf u \u2032 \u2032 2 \u0303 $, *u* is the *x*-direction velocity component, *p* is the pressure, and *ND*_{u} represents the effective numerical dissipation or diffusion for this equation. For mean mass fraction transport,

where *m*_{1} is the mass fraction of the heavy fluid. For turbulent kinetic energy transport,

where kinetic energy is $ E \u2033 = u i \u2033 u i \u2033 /2$. For mass fraction variance transport,

The mean translational velocity *u*_{trans} of the layer

accounts for the small nonzero net velocity of the mixing layer in some simulations, for example, in the quarter scale case of Flamenco, there is a residual *ũ* ≈ 0.5 m/s attributed to imperfect boundary conditions (shock wave evacuation outside the computational domain). Each term is computed principally via sixth order central differences, but with biased stencils near the boundaries.

### C. Nondimensionalization

The nondimensionalization follows the analysis outlined by Thornber *et al.*^{16} These include estimates of the initial growth rate $ W \u0307 0 $ using a power-weighted mean wavelength $ \lambda \xaf =2\pi / k \xaf $ and mean postshock density $ \rho c + $, which are defined as

where the postshock standard deviation of the perturbation amplitude is $ \sigma 0 + =C \sigma 0 $, an approximate compression factor^{1} is *C* = (1 − Δ*u*/*U*_{i}), where *U*_{i} = 434 m/s is the velocity of the incident shock, Δ*u* = 291.575 m/s is the velocity impulse imparted by the shock; postshock densities were 5.22 and 1.80 kg/m^{3} for the heavy and light fluids, respectively, giving $ \rho c + =3.51$ kg/m^{3}, *u*_{c} = $ W \u0307 0 $ = 12.649 m/s, and $ \lambda \xaf =0.2571$ m. The layer center is defined as the location at which $ m \u0303 1 ( x , t ) =0.5$.

Three measures of self-similarity may be employed. The first is to examine whether the layer properties at multiple time instants collapse when normalized by characteristic length, time, and mass scale. To examine this, some figures are plotted using axes scaled with the integral width

and $ W \u0307 $(*t*), the values of which are given in Table I. This table also gives a second measure, $ log 2 ( 3 W / \lambda \xaf ) $, which is an indicator of the number of modal generations. Elbaz and Shvarts proposed that a mixing layer must achieve ≥3 generations to reach self-similarity.^{10} Here, approximately 2.29 generations are reached, an indication that the layer may not be in the fully self-similar regime. A third measure of self-similarity is an integral mixing measure such as the “molecular mixing”

which varies by less than 1.5% in the last 3 s of the Flamenco simulation. The first measure will be explored further within this paper. For the rest of the paper, quantities are plotted against *x*/*W*, where *x* = 0 is the center of the mixing layer.

t (s)
. | τ
. | W (m)
. | $ log 2 ( 3 W / \lambda \xaf ) $ . | $ W \u0307 $ (m/s) . |
---|---|---|---|---|

0.01 | 0.49 | 0.076 05 | −0.17 | 4.12 |

0.025 | 1.23 | 0.110 7 | 0.37 | 1.436 |

0.5 | 24.6 | 0.235 | 1.46 | 0.122 |

1 | 49.7 | 0.283 3 | 1.73 | 0.077 9 |

2 | 98.5 | 0.344 3 | 2.0 | 0.049 2 |

4 | 197 | 0.419 5 | 2.29 | 0.029 56 |

t (s)
. | τ
. | W (m)
. | $ log 2 ( 3 W / \lambda \xaf ) $ . | $ W \u0307 $ (m/s) . |
---|---|---|---|---|

0.01 | 0.49 | 0.076 05 | −0.17 | 4.12 |

0.025 | 1.23 | 0.110 7 | 0.37 | 1.436 |

0.5 | 24.6 | 0.235 | 1.46 | 0.122 |

1 | 49.7 | 0.283 3 | 1.73 | 0.077 9 |

2 | 98.5 | 0.344 3 | 2.0 | 0.049 2 |

4 | 197 | 0.419 5 | 2.29 | 0.029 56 |

### D. Grid convergence

Grid convergence of the mean flow quantities has been demonstrated by Thornber *et al.*^{16} for a highly resolved simulation where all length scales are four times larger (relative to the domain size). In the current study, higher order statistics are examined; thus, it is important to revisit grid convergence to ensure that these metrics are also reasonably converged. The Appendix plots all of the individual terms of each transport equation for the code Flamenco at grid resolutions ranging from 180 × 128^{2} to 720 × 512^{2}, in Figs. 27 and 28. The overall agreement at the two highest grid resolutions for all transport terms is excellent, indicating that the results are reasonably grid converged.

The spikes of heavy fluid are intense vortical structures of relatively small size (compared to the grid spacing) and as such are the most challenging “large” scale features to converge. Given even small differences in the trajectory at early times, strong interactions of the spike’s ring-shaped vortices with other spikes and the mixing layer can cause large differences at late times. This can be seen as a region of relatively lower convergence at *x*/*W* ≈ 8 for the kinetic energy and kinetic energy balances in Figs. 28(c) and 29(c). Such an intense impact of the spikes on the kinetic energy budget has also been observed in an analysis of the RMI spectral kinetic energy budget.^{56}

The levels of noise in each of the quantities are a good indicator of statistical convergence, as no smoothing has been applied. However, the statistical noise does not mask the overall trends for the majority of the results presented. Based on the previous study,^{16} it is expected that results from Turmoil and Triclade would converge at a similar grid resolution to Flamenco. Flash, however, has historically demonstrated slower convergence compared to the other algorithms, but given the good agreement with the other algorithms presented here, and the fact that Flamenco is reasonably converged at the intermediate grid level (a factor of eight times fewer points than the highest resolution), the Flash results are also expected to be reasonably converged.

### E. Note on the treatment of numerical results

Compiling the data from each code highlighted numerical and physical issues at different times. Flamenco exhibits postshock oscillations at the earliest times which are clear in the plots of *ũ*, Triclade and Turmoil show highly fluctuating acoustic pressure correlations at late times (terms such as Π^{E″}), and Flash has difficulties at the spike side boundary at late times (visible in the spike side $ E \u2033 \u0303 / u c 2 $ at *t* = 2 s, for example). An effort has been made here to present all results except where (i) the result is so noisy that it obscures interpretation of the rest of the data or (ii) it is clearly dominated by a nonphysical effect. There are some cases which are marginal (perhaps unphysical) but can fit on the plot without obscuring the data. Where data are deliberately excluded due to such effects, the omission and reason for omission are stated in the figure caption.

## III. RESULTS

### A. Flow features at four time instants

Visualization of the distribution of the volume fractions within the flow field is shown for the four time instants chosen using slices in the *y*-plane in Fig. 2 and three-dimensional isosurfaces of volume fraction in Fig. 3 from the Flamenco code.

At the initial time *t* = 0.01 s, the individual modes at the highest wavenumber have already formed characteristic mushroom shaped bubbles (light fluid penetrating the heavy) and spikes (heavy fluid penetrating the light). As expected at this Atwood number, several spikes have propagated substantially beyond the bulk of the mixing layer. Two of these are visible on the left-most image in Fig. 2; however, due to the perspective, they are not as clearly seen in the three-dimensional rendering.

At *t* = 0.025 s, the layer has begun to transition, where the spikes contain relatively well mixed fluid, but the bubble heads still contain some unmixed material. This provides a preliminary indication that turbulent fluctuations are initially stronger on the spike side of the mixing layer. There is visible evidence of the vortex projectiles propagating far from the layer on the spike side; however, their intensity and size have reduced. At the two latest times, there is a core of well mixed fluid, but the flow is still relatively nonuniform in the homogeneous direction. The layer width grows by ≈20% between *t* = 1 and 2 s, and low wavenumber modes dominate at *t* = 2 s. Strong three-dimensionality and a range of vortex length-scales are apparent at *t* = 1 and 2 s in Fig. 3; however, there is still evidence of coherent spike structures such as the ring vortex at the upper left of the three-dimensional visualization at *t* = 2 s.

Figure 4 plots the integral width and molecular mixing measure Θ at early and late times. Superimposed on this figure are the chosen sampling times. The first two sampling times are chosen to be just after the shock passes through the interface and at a short time later. At this stage, the layer is expected to be formed of coherent structures and a wide range of turbulent length scales is not expected. The two latest times are chosen to be at a time where the layer may be approaching self-similarity, but with a compromise between having the widest possible range of turbulent length scales and a moderate impact of a lack of statistical accuracy in the homogeneous direction and domain-size constraint.^{57} From *t* = 1 to 4 s, Θ changes by only 1.5% in Flamenco simulations, an indication that it is approaching a self-similar state (but not yet reached, as discussed earlier).

### B. Momentum transport

Figure 5 plots the profiles of *ũ*, *ṽ*, and $ w \u0303 $ normalized by *u*_{c} immediately following shock interaction (0.01 s), a short time later where mixing is increasing rapidly (0.025 s), at the onset of approximate self-similarity (1.0 s), and at the latest time (2.0 s). The overall agreement for *ũ* is very good across the codes, particularly, at late times, with the main discrepancies at early times where the peak *ũ* varies by 14% from the average peak value. The Favre-averaged velocity is positive, representing the net flux of heavy material into the light material as expected. At early times, the profile of *ũ* is asymmetric, with a larger gradient on the bubble side (*x* < 0) and more gradual decrease on the spike side (*x* > 0); however, at late times, it is reasonably symmetric and the peak velocity is at the centerline (defined by mean mass fraction of one half). Asymmetries are expected at this moderate Atwood number.

The Favre-averaged in-plane velocities *ṽ* and $ w \u0303 $ are relatively large compared with *ũ* at late times, where it could be expected that these quantities should average to zero in the limit of an infinite domain. As the domain size is finite, these components may take on values roughly equivalent to the noise observed in the evaluation of *ũ*. This is clearly not the case as the profile of *ũ* is relatively smooth. The large magnitude of *ṽ* and $ w \u0303 $ also extends well into the spikes (*x*/*W* > 4), where the plane-averaged heavy mass fraction is <1%, and these fluctuations are associated with the propagation of “vortex projectiles” exiting the mixing layer. As this is a turbulent mixing layer, it is not surprising to see differences in the locations of the peaks and troughs for the in-plane components between the algorithms, especially given that the mean *ṽ* and $ w \u0303 $ are 0.6% of the peak values at a given *x* location.

The relatively large magnitude is most likely due to the mean gradient of *ũ* in the *x* direction being constrained to be close to zero due to approximate incompressibility and the corresponding divergence-free constraint. Due to numerical miscibility, the current simulations will not satisfy ∇ · **u** = 0, as ∇ · **u** = −*D* ln *ρ*/*Dt*, where *D*(·)/*Dt* is the material derivative. Although the flow at late times is at low Mach number and nearly incompressible, the advected density changes due to numerical mixing. The periodic boundary conditions in the *y* and *z* directions do not constrain the net flow in that direction, thus leading to the observed profiles.

Late time self-similarity is examined in Fig. 6, where the Flamenco data have been nondimensionalized by the integral width *W*(*t*) and the time rate of change of the integral width $ W \u0307 $(*t*). An additional point at *t* = 4 s has been added, where at that time, spikes have left the computational domain, but the bulk of the mixing layer may be assumed to be reasonably well resolved. The momentum profiles collapse well for the bulk of the layer under the chosen scalings for *t* > 1 s; however, the spike side is taking longer to achieve approximate self-similarity.

Figure 7 shows the individual components of the mean momentum transport equation for all codes, Fig. 8 shows the same time instants but only the Flamenco results so that the relative magnitudes of the individual terms may be more clearly discerned, and Fig. 9 plots the *F*^{u} term alone, where the data have been smoothed by a simple five point moving average. The codes agree well for all terms except for the pressure gradient term *F*^{u}, which shows good agreement on the bubble side of the mixing layer but relatively poor agreement on the spike side. Although the bulk of the mixing layer is contained in −4 ≤ *x*/*W* ≤ 4, there is a substantial contribution to the mean momentum equation at *x*/*W* > 4 which may be attributed to highly energetic ejected spikes and features which combine both very high velocities and low pressure cores in a localized feature. Such features have been observed in previous analysis of spectral energy transfer in the RMI^{58} and RTI.^{59}

As has been shown for a small Atwood number Rayleigh-Taylor instability,^{15} the time rate of change of mean momentum *t*^{u} is mostly generated as the remainder of two larger terms, the mean pressure gradient *F*^{u} and the Reynolds stresses *R*^{u}, where *F*^{u} ≈ −*R*^{u} at all times. At early times, both *R*^{u} and *F*^{u} are approximately antisymmetric across the centerline, but at later times, the bubble side amplitudes are larger and closer to the core of the mixing layer than the spike side amplitudes, as the gradient of the Reynolds stress *τ*_{11} is much sharper on the bubble than the spike side (see Sec. III E). Mean advection is nearly negligible at all times, which is expected as the mean layer position is approximately stationary with respect to the chosen reference frame.

The mean momentum transport equation has contributions from numerical dissipation *ND*_{u} representing the under-resolved components of the Reynolds stresses. As the grid resolution increases, the magnitude of these implicitly modeled effects should decrease. Figure 10 plots the nondimensional numerical dissipation term for Flash, Flamenco, and Triclade, computed as the residual of the mean momentum equation. Although the overall magnitude is similar to the other terms in the momentum equation, the magnitude is largely due to the presence of noise in the computation of the individual physical terms as shown clearly in Fig. 7. Examining only the Flamenco results (which have relatively low noise) shows that the peak magnitude is less than 3.5 × 10^{−5}, which is ten times lower than the dominant physical terms, and the rms magnitude defined over −5 < *x*/*W* < 5 is 7.6 × 10^{−6}. A similar conclusion can be reached if a simple moving average is applied to the results. It should be noted however that as the net change of mean momentum is a result of the difference between the two largest physical terms *F*^{u} and *R*^{u}, the numerical dissipation represents approximately one quarter of the magnitude of the overall rate of change of *ũ* and is not negligible.

### C. Mass fraction transport

Figure 11 presents the profile of the mean heavy fluid mass fraction $ m \u0303 1 $. The profile is approximately linear and at late times exhibits no discernable kinks indicating a region of well mixed flow, with very small mean diffusive transport. All codes are in good agreement barring some late time differences on the spike side. Figure 12 shows the scaled $ m \u0303 1 $ profile from Flamenco at four times, demonstrating an excellent collapse for the three latest times on the bubble side in particular. Figure 12 also shows a close-up of the spike side, where the size of the spikes relative to the mixing layer is still decreasing at late times, a further indication that a fully self-similar state is not yet attained. This is consistent with breakup and dissipation of these strong vortical structures.

Figure 13 plots all terms in the $ m \u0303 1 $ transport equation, where all code results collapse well. Turbulent transport *T*^{m1} is approximately antisymmetric about the center of the mixing layer and dominates for *x*/*W* < 0; however, due to mean advection *A*^{m1} being negative throughout the layer, the time rate of change *t*^{m1} dominates for *x*/*W* > 0. Mean advection *A*^{m1} is non-negligible, symmetric at the two latest times, and negative throughout the layer due to the net transport from the heavy to light side, but is less than half the magnitude of turbulent transport. As a result, the time variation of $ m \u0303 1 $ largely follows *T*^{m1} but is further skewed toward the heavy side due to the contribution from *A*^{m1}.

Similar to the mean momentum equation, the results presented here also include numerical dissipation *ND*_{m1}, which is estimated as the remainder in the governing equations and shown in Fig. 14 for *t* = 1 s. Again, the result is quite noisy, and the peak of 9 × 10^{−5} is at least one order of magnitude lower than the peak of the dominant physical terms. There is no discernable structure to the numerical dissipation profile, exhibiting as symmetric noise about zero.

### D. Mass fraction variance transport

Profiles of the heavy species mass fraction variance $ m 1 \u2032 \u2032 2 \u0303 $ are shown in Fig. 15 for four representative times. At early times, the profiles peak moderately on the spike side of the mixing layer and are relatively symmetric. There is a substantial difference in the peak mass fraction variance from each code, bounded by Flamenco as the lowest and Turmoil as the highest. This is expected as at early times the layer is in a state of transition, where the initially smooth density variation across the interface is being sharpened to a grid-scale discontinuity by shear flows driven by the growth of RMI. The minimum thickness of the layer is a function of the dissipation within each code, hence the increased uncertainty. At the latest two times, the mean gradients within the layer are smoother and the codes are in better agreement, although the peaks are still only in agreement to within ±10%. As time progresses, the variance decreases toward a self-similar profile as stirring reduces the gradients through mixing.

Figure 16 shows the collapse of the $ m 1 \u2032 \u2032 2 \u0303 $ profiles at several times. Except for the earliest “transitional” time of 0.5 s, the profiles at the other three times collapse well, having a similar peak value of ≈0.046. Mass fraction variance gradients are steeper on the spike side of the layer compared to the bubble side.

Individual terms in the transport equation for $ m 1 \u2033 \u0303 $ are plotted in Fig. 17 for all codes. At the two earliest times, all terms but mean advection are important. At the latest time, mean production dominates all other physically resolved terms by a factor of three, although it must be noted that mean numerical dissipation is of a similar magnitude. Turbulent transport increases variance at the bubble and spike fronts but reduces it within the core. Although this also occurs at the later times, it is dominated by mean production, which is asymmetrically weighted toward the bubble side of the layer.

The numerical dissipation *ND*_{m1″} has been computed for *t* = 1 s and is plotted in Fig. 18 for Flamenco, Flash, and Triclade, alongside the production over dissipation for Flamenco. The agreement between the codes is excellent, particularly, given the different algorithmic approaches. As described above, numerical dissipation is of similar magnitude to the largest physical terms and is purely dissipative. The ratio of production to dissipation peaks at ≈1.5, somewhat larger than the Rayleigh-Taylor instability value of 1.25 reported by Schilling and Mueschke^{41} and takes a volume weighted average of 0.92 over the region −5 < *x*/*W* < 5.

### E. Turbulent kinetic energy transport

The dimensionless turbulent kinetic energy $ E \u2033 \u0303 $ is shown in Fig. 19 for four times. At the earliest time, all codes show a “double peak” structure, where the peaks are focused on the bubble and spike heads, respectively. The highly energetic spikes persist at late times, causing a broadening of the profile into the light fluid (*x*/*W* > 0). The bulk of the turbulent kinetic energy content in the bubbles is restricted to the region *x*/*W* > −6 and increases steeply on the bubble side. This is consistent with the visualizations in Fig. 2, which shows an earlier onset of mixing on the spike side (an indicator of stirring by a range of length scales) than on the bubble side. The advection of spikes away from the layer at early times smooths the spatial contribution of the spike to the turbulent kinetic energy profiles, and turbulent transport acts to merge the early “double bump” profile into a single peak. The peak is slightly biased to the spike side of the mixing layer at late times, where most of the kinetic energy is located. For example, at *t* = 1 s and *x*/*W* = 5.5, there is substantial $ E \u2033 \u0303 $ even though $ m \u0303 1 <0.1$%. At late times, the profile is asymmetric. In comparing the codes, there are consistent peak positions in Triclade, Turmoil, and Flamenco. Flash shows large values on the spike side at the latest time, which is attributed to unphysical interactions with the boundaries.

The scaling of the profiles of $ E \u2033 \u0303 $ with the mixing layer width and its growth rate is plotted in Fig. 20. Overall, the collapse is excellent for the core of the layer, but notable differences can be seen in the prediction of the peaks in $ E \u2033 \u0303 $ in the spikes. As observed previously, the approximate vortex ring structures do not scale with the same power law as the rest of the mixing layer and generate isolated peaks and troughs in the $ E \u2033 \u0303 $ profile which are not consistent from one time to the next when scaled using the integral layer properties. However, the overall shape of the profile is in good agreement given the factor of 8 reduction in dimensional kinetic energy between the first time and the last shown in Fig. 20.

There are six terms in the kinetic energy budget, of which the time variation *t*^{E″}, pressure-dilatation correlation Π^{E″}, and turbulent transport *T*^{E″} terms are relatively large, and mean advection *A*^{E″}, buoyancy production $ P b E \u2033 $, and shear production $ P s E \u2033 $ are relatively small. The large terms are plotted together in Fig. 21 and the small terms in Fig. 22. To aid a clear interpretation of the relative magnitude of each of the terms, the full budget is plotted for the Flamenco results only in Fig. 23. Finally, $ P b E \u2033 $ and *T*^{E″} have strong fluctuations; thus, Fig. 24 plots those two terms independently at *t* = 1 and 2 s where a five point moving average has been applied to smooth the data.

First, a note on some numerical issues in computing specific terms. The overall level of fluctuations in the results is substantially larger than in the transport equations for $ m \u0303 1 $ and *ũ*, which is partly expected as $ E \u2033 \u0303 $ is a second order quantity; however, there were also specific difficulties in the computation of terms involving correlations of pressure fluctuations. This led to a large variation in results from code to code for the pressure-dilatation Π^{E″} and turbulent transport *T*^{E″} terms.

Several different stencils were employed to compute the pressure-dilatation term, ranging from second to sixth order to explore the convergence of the numerical results. For Flamenco, it was noted that both the magnitude and sign changed for the pressure-dilatation correlation from second order estimation of derivatives to sixth order. It is encouraging that the third and fourth order stencils are converging toward the sixth order result.

Now, consider the individual terms in Fig. 21. At the first and second times, turbulent transport *T*^{E″} forms two symmetric peaks at the bubble and spike fronts, indicating a transport/redistribution of turbulent kinetic energy from the core of the layer to the edges. Turbulent transport is the largest term contributing to the time rate of change of $ E \u2033 \u0303 $. The peak in *T*^{E″} at the bubble side remains at later time; however, the peak on the spike side advects away from the layer, becoming flatter and indicating a net transport of $ E \u2033 \u0303 $ from the spike to bubble side.

The pressure-dilatation term is symmetric, peaking at the center of the layer, and increases $ E \u2033 \u0303 $. As pointed out by Schilling and Mueschke,^{41} although the Mach number is low, mean dilatation is non-negligible due to mixing of the fluids; thus, ∇ · **u** = −∇ · (*D*∇ ln *ρ*), with *D* being an effective numerical diffusion, the form of which may be more complex than Fickian.

Both Turmoil and Triclade results exhibited large fluctuations in pressure fluctuation-based terms which precluded the extraction of clear trends from their data. This is due to the treatment of acoustic waves in both of those algorithms, which both dissipate acoustic fluctuations less strongly than Flamenco and Flash. Acoustic waves add a “fast” component to the kinetic energy budget, but are a mostly isentropic process, and thus should not impact the overall time rate of change of $ E \u2033 \u0303 $. Figure 19 shows that this is the case: the peak kinetic energies have reduced by a factor of 1000 over the simulation time, yet all codes remain in good agreement; thus, net dissipation must be of a similar order of magnitude for all algorithms and very likely to be scaling relative to a physical (rather than numerical) dissipation rate. This division of the pressure-dilatation into compressible and incompressible components is in agreement with previous observations for homogeneous decaying turbulence and shear layers,^{18,60,61} where it was noted that although the fast component has a large magnitude, its global impact is small compared with the “incompressible” correlation.

Figure 25 plots the spatial distribution of the rms pressure variance $ \u27e8 p \u2032 2 \u27e9 $ at several times using Flamenco and the temporal variation at the center of the mixing layer [defined by *x* such that $ m \u0303 1 ( x , t ) =0.5$] and the integral average in the heavy bulk fluid [*x* such that $ m \u0303 1 ( x , t ) >0.99$] and light fluid [*x* such that $ m \u0303 1 ( x , t ) <0.01$]. At late times, the pressure fluctuations within the mixing layer in Flamenco clearly follow an incompressible scaling ≈21.5*ρ*_{av}$ W \u0307 2 $ with $ \rho a v = ( \rho H + + \rho L + ) /2$. This is a reasonable result noting that the integral width is substantially smaller than the actual width of the layer which leads to the large constant of proportionality. However, for Triclade and Turmoil, the pressure fluctuations are substantially larger.

Focusing on the fluctuations in the pure heavy fluid in Fig. 25(b) and pure light fluid in Fig. 25(d), the fluctuations are modeled as a simple spherical source deposited at shock-interaction, where given conservation of acoustic energy, the pressure fluctuations should scale with $ \u27e8 p \u2032 2 \u27e9 ( t ) \u2248 \u27e8 p \u2032 2 \u27e9 ( t 0 ) r 0 ( t 0 ) / [ r 0 ( t 0 ) + a H , L t ] $. Here, *a*_{H} and *a*_{L} are the speed of sound in the shocked heavy/light fluid, respectively, *t*_{0} is a sampling time measured from shock interaction, and the initial acoustic radius is modeled as $ r 0 ( t 0 ) = \lambda \xaf /2+ a H t 0 $.

The overall agreement is very good from *t* = 0.01 to 0.1 s between Flamenco, Triclade, and Turmoil, despite the complexities of the acoustic waves and the influence of the exiting rarefaction at early times (*t* < 0.1 s). At later times, numerical dissipation in Flamenco reduces the magnitude of the acoustic waves in the heavy and light fluids, to below the pressure fluctuations produced by turbulent motion in the layer. This is not seen in Turmoil and Triclade, which preserve acoustic fluctuations better. One-dimensional tests indicate that Turmoil damps acoustic waves less strongly than Flamenco, which is consistent with this result, especially noting that by *t* = 1 s the wave has traveled through the cross section more than 400 times and so numerical dissipation has had ample opportunity to cause the observed differences. This explains the discrepancy between the three codes in the layer center, where pressure fluctuations in the heavy fluid penetrate the layer and are at a larger magnitude than pressure fluctuations generated through vortical motion alone.

At even later times, acoustic pressure fluctuations will dominate the vortical fluctuations by orders of magnitude. This also explains the difficulty in obtaining agreement among the codes for transport terms including pressure gradients, and correlations including pressure fluctuations and gradients. Despite this, the impact on the overall transport balance is small.

The time derivative *t*^{E″} is plotted in Fig. 21 and shows an increase in turbulent kinetic energy at both the bubble and spikes at the earliest times. However, at *t* = 1 and *t* = 2 s, the time rate of change is negative throughout the core of the mixing layer, with positive variations only at the isolated vortex “projectiles.” Note that the overall balance is impacted by dissipation provided by the implicit numerical model, which is discussed at the end of this section.

Examining the smaller terms plotted in Fig. 22, mean advection *A*^{E″} shows a transport of kinetic energy from the bubble to the spike side of the layer. Both buoyancy and shear production are negative on the bubble side and positive on the spike side, and all three terms are only nonzero within the core of the mixing layer. Some oscillations in $ P s E \u2033 $ may be seen at *x*/*W* > 8 at *t* = 2 s for Flash, due to the unphysical interaction of the spikes with the boundary condition.

The relative magnitude of all six terms is shown clearly in Fig. 23 for the Flamenco results only. It is clear that the time evolution of $ E \u2033 \u0303 $ is dominated by two physical effects, the turbulent transport and the pressure-dilatation correlation, and by numerical dissipation with a mean magnitude of the same order as the resolved physical terms.

The numerical dissipation *ND*_{E″} has been computed for *t* = 1 s and is plotted in Fig. 26. This figure compares the numerical dissipation computed for Flamenco and Flash, where the Triclade results are not shown as they are dominated by the fluctuations in the pressure-dilatation term due to the acoustic field. As mentioned in the previous paragraph, the numerical dissipation of turbulent kinetic energy is of a similar order to the dominant terms in the transport balance. This large order of magnitude is to be expected. Turbulent kinetic energy is conserved until it reaches the viscous scales, where it is dissipated. In an implicit LES, turbulent kinetic energy is dissipated by numerical dissipation as it reaches the grid scale. The mechanism of turbulent dissipation is highly code-dependent, with some algorithms having dissipation in the remap phase for low Mach number flows (e.g., Turmoil), and others using upwind schemes to dissipate on the individual wave strengths (e.g., Flamenco, Flash, and Triclade).

However, implicit LES relies on the concept that the actual dissipation rate is driven by the largest scales of turbulence, and independent of the detailed mechanism of dissipation as long as the dissipation does not impact the large scales. Thus, when a sufficiently wide range of large scales is resolved by a numerical algorithm, the flow would evolve in a statistically identical manner regardless of whether the smallest scales of the flow are resolved or not. The word “statistical” is used since turbulent flows are very sensitive to small perturbations; hence, an exact match is not expected.

That the four codes used in the current study have such good agreement in the turbulent kinetic energy profiles at a very late time, where peak values have reduced by a factor of 1000, indicates that the decay rates are very similar from algorithm to algorithm. Thus, although only results from two codes are plotted in Fig. 26, there is reasonable confidence that a similar magnitude would be observed for the other two codes. Finally, the grid convergence of *ND*_{E″} for Flamenco as shown in Fig. 27 also indicates that the simulations have achieved the necessary scale separation such that the total dissipation rate is determined by processes at resolved scales, rather than by numerical dissipation.

## IV. CONCLUSIONS

An analysis of turbulent transport in an inhomogeneous compressible turbulent mixing layer induced by Richtmyer-Meshkov instability has been presented. The analysis used data from four independent ILES codes which were applied to the quarter-scale problem presented in the *θ*-group collaboration.^{16} All terms in the high-Reynolds number (inviscid) limit of the Reynolds-averaged transport equations for momentum, species mass fraction, turbulent kinetic energy, and mass fraction variance were computed in the inhomogeneous direction. Numerical dissipation provided an implicit subgrid model in all algorithms, which was computed from the balance of the transport equations. The objective of this study was to provide insight into the dominant physical mechanisms influencing turbulent transport and to facilitate the future development and validation of reduced order modeling of those quantities. The inclusion of four numerical methods also provided an understanding of the confidence in predicting each of the terms, and to compare and contrast the impact of numerical dissipation.

For the momentum transport equation, there is high confidence in the mean momentum in the turbulent stage, with codes agreeing to within ≈±4% in the shock-direction. The time rate of change of momentum is dominated by the balance of the pressure gradient and Reynolds stresses, and the numerical dissipation is very low. The pressure gradient term was particularly noisy due to acoustic modes. Advecting spikes which escape the mixing layer cause large variations in the individual transport terms, the position and magnitude of which were not easily captured numerically. The results here may also be impacted by the small sample size of such spikes which escape the developing mixing layer at an early time. The individual transport terms have a spread of ≈±10% due to statistical errors and acoustic fluctuations, but the agreement between codes is good and could be employed to validate URANS and LES modeling approaches.

For mass fraction transport, the mean profiles are nearly identical at all times; thus, there is a high level of confidence. The time variation of mass fraction distribution largely follows the turbulent transport term, moderated by mean advection, and numerical dissipation is small. Here, there is very good agreement between all algorithms, and thus a high level of confidence in the interpretation of the dominant terms. The mass fraction variance shows a higher variation at early times, but again agreement is very good when the layer transitions to turbulence (≈±5%). All terms in the transport equation are important, but mean advection is lower by a factor of approximately two. Numerical dissipation is purely negative as expected and is of a similar magnitude to the largest resolved terms. All terms considered are in sufficient agreement to be useful for validation of lower order modeling results.

The analysis of the turbulent kinetic energy transport equation was more challenging. The largest error in turbulent kinetic energy occurred at the earliest, nonturbulent times (≈±10%), but with excellent agreement at later times (≈±3%). Temporal evolution is dominated by turbulent transport and numerical dissipation (implicit subgrid dissipation), with pressure-dilatation the third largest. The computed pressure-dilatation term was strongly impacted by acoustic modes, rendering the confidence in the results for that term relatively low compared to the rest of the data, given the numerically challenging problem of evolving sound waves for relatively long periods. Despite this uncertainty, the impact on the overall development of kinetic energy in the layer is concluded to be low, as evidenced by the excellent agreement in late time kinetic energy for algorithms which have substantially different mean pressure fluctuations at late times.

## ACKNOWLEDGMENTS

The authors greatly acknowledge Dr. Karnig Mikaelian who suggested the original “*θ*-group” collaboration at the Fourteenth International Workshop on the Physics of Compressible Turbulent Mixing in 2014. This research was supported under the Australian Research Council’s Discovery Projects funding scheme (Project No. DP150101059). The lead author would like to acknowledge the computational resources at the National Computational Infrastructure through the National Computational Merit Allocation Scheme which were employed for all Flamenco cases presented here. Flash was developed by the DOE-sponsored ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344.

### APPENDIX: GRID CONVERGENCE

Due to the large number of individual terms computed, convergence is demonstrated here through the examination of the Flamenco results alone, as detailed in Sec. II D, using simulations ranging from 180 × 128^{2} to 720 × 512^{2}. Figure 28 shows the convergence of mean properties at *t* = 1 s. The two highest grid levels are in very good agreement. Figure 29 shows the grid convergence of the individual terms in each transport equation for Flamenco. Again, the two highest resolutions agree well, particularly considering the higher-order nature of the metrics under examination. Finally, Fig. 27 shows the convergence of computed numerical dissipation for Flamenco. For the mean momentum and mass fraction, numerical dissipation is small and dominated by noise. For the mass fraction variance and kinetic energy, numerical dissipation is significant and the difference between the two highest grid resolutions is again small.