Late-time growth rate, mixing, and anisotropy in the multimode narrowband Richtmyer–Meshkov instability: The θ-group collaboration Free

Richtmyer–Meshkov instability (RMI) occurs when a perturbed interface between two materials is accelerated impulsively, usually by an incident shock wave.^1,2 During the passage of the shock-wave, a layer of vorticity is deposited by baroclinic source terms at the interface between the gases and in the nearby region as the shock returns to planarity. This deposition causes a net growth of the layer, which during the linear stage for a single mode may be estimated by $ȧ≈kAt+a0+Δu$⁠,¹ where $ȧ$ is the rate of increase of the amplitude of the perturbation, k is the wavenumber, At⁺ is the post-shock Atwood number, $a0+$ is the post-shock amplitude, and Δu is the velocity jump imparted by the incident shock.

The linear growth rate reduces as the amplitude of the perturbation becomes large relative to the wavelength, and characteristic bubble/spike (mushrooms) features develop. At later times, shear layers along the contact surface break down due to Kelvin–Helmholtz instabilities that drive a transfer of kinetic energy from the shock-parallel direction to the perpendicular components. It is reasonably well established that a multimode perturbation will transition to a turbulent mixing layer with a width $∝tθ$⁠; however, the exact behavior of the layer is complicated by a strong dependence on initial conditions and chaotic turbulence physics. For a comprehensive and up-to-date survey of research on Richtmyer–Meshkov and the closely related Rayleigh–Taylor instability, see Zhou.³ 

RMI occurs in man-made events such as inertial confinement implosions (e.g., Ref. 4), explosions (e.g., Ref. 5), fuel injection in supersonic flows,⁶ and in large astrophysical events such as supernovae (e.g., Ref. 7). Such events are by their nature inaccessible to detailed experimental diagnostics due to short time scales, extremes of temperature and pressure, or very small or very large observational distances. Thus, the understanding of RMI relies on a judicious use of theory informed by careful experimental and numerical studies, the latter of which is the focus of this paper.

There have been multiple numerical studies of three-dimensional multimode RMI that have improved our understanding of this complex flow.^8–17 Each study has explored the integral properties of the flow field, such as width, mixedness, and fluctuating kinetic energy, along with higher order quantities. Even for the lower order quantities, a consensus on the measured physical behavior has not been achieved, given the sensitivity of the late time behavior to the initial conditions and the challenging flow physics for even state-of-the-art algorithms.

As summarized by Zhou,³ the growth exponent θ for the mixing layer width lies between 0.213 and 2/3 in theoretical, experimental, and numerical studies. Of the ten prior studies, four of them have similar idealised narrowband initial conditions and thus may be compared;^9,11,14,17 however, the other studies have individually tailored the initial conditions to match specific experiments or for numerical reasons. Thus, a direct comparison of integral properties between all prior studies is not possible.

Furthermore, of all prior studies, only four include computations using more than one algorithm for the same initial condition,^11,13,14,18 each of them using two independent algorithms. It is not possible to obtain a robust estimate of numerical uncertainty with such a limited sample. As all other studies used differing initial conditions, information about numerical uncertainty cannot be obtained from the existing literature. This is particularly important as all algorithms are challenged by the need to (i) capture shock waves, (ii) capture moving contact surfaces, and (iii) represent fine scale turbulent motion.

This paper considers RMI generated from an initial high wavenumber perturbation on a contact surface. Although most realistic problems have a larger range of initial perturbation lengthscales, this perturbation form is of interest as it represents a lower limit on the growth rate due to RMI, since the late time growth at large lengthscales is driven solely through mode coupling. Based on the previous literature, it is expected that the growth exponent of the layer will be $≈0.213≤θ≤0.295$⁠. Linked to the growth rate is the mixedness of the layer, where generally a higher θ implies a less well-mixed layer. Mixedness is particularly challenging to resolve since it depends on resolving a large enough portion of the sub-inertial range such that sub-cell variations in concentration are negligible (without a sub-cell representation of scalar variance). Finally, previous studies have highlighted significant asymmetry in the self-similar turbulent kinetic energy components; however, late time self-similar anisotropy has not yet been conclusively demonstrated.

Here, a new initial condition is proposed, which combines the most favourable features of prior studies such that simulations may be run using a wide range of algorithms ranging from Arbitrary-Lagrangian-Eulerian (ALE) and Godunov methods to algorithms based on compact differences. Simulations are then run with eight independent algorithms for an initial condition for which a grid converged solution may be achieved.

Comparing results from many different codes simulating the same mixing flow has already proved fruitful in the past. The α-group collaboration¹⁹ analysed the results of seven different codes in the context of Rayleigh–Taylor turbulent mixing. The conclusions were instrumental in understanding the effect of the initial perturbation spectrum on the value of the turbulent Rayleigh–Taylor growth parameter α [in $W(t)=αAgt2$ where W(t) is the mixing zone width, $A$ is the Atwood number, and g is the acceleration of the interface].

A further set of simulations was then run for a more challenging condition that aimed to achieve later dimensionless times than have been presented to date. The results are then employed to give the current best-estimate of the key integral properties of the flow including the growth rate, integral mix measures, and kinetic energy. The best estimates are accompanied by uncertainty estimates for each, which for the first time incorporate algorithmic uncertainty.

The layout of this paper is as follows. Section II details the initial conditions, domain sizes, boundary conditions, algorithms, diagnostics, and the grid convergence study. Results for the standard case are presented in Sec. III and the one-quarter length-scale case in Sec. IV. Finally, the key conclusions are summarised in Sec. V.

The objectives of the formulation of the initial conditions are to provide a platform from which an effective collaboration could be constructed to address gaps in understanding or regions of uncertainty and to provide a reliable estimate of the growth exponent θ for a specific initial condition. Key to the formulation was to ensure that the initial condition could be run consistently in a range of numerical algorithms spanning Eulerian, Lagrange-remap through to compact differences. It should also be computationally feasible, i.e., a grid converged solution may be achieved at a resolution accessible by a wide range of codes. The choice of initial conditions builds on previous studies^{11,13,18,20,21} and has the following features:

• An initial range of length scales in the perturbation from L/8 to L/4 where L is the cross section. This should ensure that the layer growth is reasonably converged at $≈1283$ and permits a future Direct Numerical Simulation (DNS) at an achievable resolution $≈10243$⁠.

• A diffuse initial interface of error function form and thickness L/32. This leaves open a future option to run a DNS of this configuration and is accessible for Large-eddy simulations (LES) that need an initial resolved interface thickness, whilst being narrow enough to not impact greatly the evolution of the instability.

• A power spectrum with constant power at all initialised wavelengths and an overall amplitude of $0.1λmin$⁠.

• Mode amplitudes and phases defined using deterministic random numbers. These are the same at all grid resolutions such that a grid refinement study can be achieved with the same interface shape.

• A moderate density ratio of 3:1, γ = 5/3 for both gases, and a Mach 1.84 shock.

The test case uses the initial conditions derived by Ref. 9, which are shown schematically in Fig. 1. The flow field consists of a heavy and light gas separated by a perturbed interface, where the perturbation satisfies a given power spectrum and mean amplitude. The heavy and light unshocked gases are at the same temperature in the undisturbed field, and the specific heats for each component (⁠$Cv1$ and $Cv2$⁠) are chosen to ensure this.

The incident shock wave has a Mach number of 1.8439, equivalent to a four-fold pressure increase. The initial conditions in the shocked domain and left and right gases are

An initial velocity −Δu = −291.575 is given to the gas interface such that the centre of the interface is stationary after passage of the shock wave, where Δu is determined from an unperturbed computation. The pre-shock densities of the heavy and light fluid are 3 and 1 kg/m³, respectively, giving an Atwood number A = 0.5. Both fluids have an initial pressure of 100 kPa. Post-shock densities were 5.22 and 1.8 kg/m³ for the heavy and light fluids, respectively, giving a post-shock Atwood number of At⁺ = 0.487. Assuming that the layer amplitude has been reduced by compression $a0+=(1−Δu/Ui)a0$⁠, where U_i is the velocity of the incident shock, gives the Richtmyer velocity $ȧ=kmaxA+a0+Δu=29.36$ m/s.

Here, the initial conditions use $λmin=L/8$ and $λmax=L/4$ and a constant power spectrum. An outline of the derivation of the perturbation in the initial condition is detailed in the  Appendix and is more extensively discussed by Youngs and Thornber et al.^9,11 The initial perturbation is defined by a perturbation power spectrum that is constant between the minimum and maximum wavelengths within the problem. The phase and amplitude of each individual mode are computed from a Gaussian distribution of deterministic random numbers where the mean amplitude satisfies the specified power spectrum. Given the perturbed interface, each algorithm requires volume fractions f₁ and f₂ = 1 − f₁ within the mixed cells. Here, a diffuse interface of gradient thickness of δ = L/32 is employed, thus f₁ is computed as

where S(y, z) = 3.5 + A(y, z), A(y, z) is the perturbation satisfying the specified power spectrum (defined in the  Appendix) and f₁ = 1 − f₂ is the volume fraction of the heavy gas.

To improve the quality of the initial condition on coarse meshes, the initial condition is evaluated in $4×4×4$ sub-cells for each computational cell and then summed to give the volume fraction for each computational cell. This was deemed to be simpler than integrating the initial condition and gives satisfactorily converged results. A visualisation of the isosurface of f₁ = 0.5 is shown in Fig. 2 for the $180×1282$ resolution. The diffuse initial condition was compared to a sharp interface in the Flamenco code, and it was verified that it does not have a large impact on the overall behavior of the mixing layer, marginally increasing the growth rate and decreasing the mix parameter Θ. Simulations using double the domain length highlighted that some spikes, which have formed strong coherent vortex rings, are expected to advect out of the domain at $t≈0.5$ s; thus, simulations that are run to t = 0.5 s should give a reasonable estimate of the asymptotic mix parameters and be reasonably independent of the boundaries (for lower order parameters). At the latest time, the integral width is less than 10% of the domain size, and it should be reasonably statistically converged and not be box constrained.¹⁷ As an additional check, the normalised two-point, in-plane longitudinal correlation,

is computed for $0≤s≤L/2$ at the plane closest to the centre of the mixing layer, defined by $⟨f1⟩=0.5$⁠, where $⟨.⟩$ indicates a plane-averaged quantity. Figure 3 plots this for five representative times during the computation, where the first zero crossing point of the correlation is below s/L = 0.16 for all cases. This confirms that the choice of the domain size is sufficient for the chosen physical time span of the simulation. However, it is important to emphasize that it may not be self-similar as it may not have had sufficient time to become fully developed.

The computational domain is Cartesian as shown in Fig. 1 and measures $x×y×z=2.8π×2π×2π$ m³. All codes but one (Hesione) in this study undertook a grid convergence study including mesh sizes $180×1282$⁠, $360×2562$⁠, and $720×5122$⁠, and TURMOIL and Triclade were also run at a grid resolution of $1440×10242$⁠. The boundary conditions are periodic in the y and z directions and outflow in the x direction. It should be noted that the boundary conditions in the direction of shock propagation vary considerably between the implementations, as it was decided that this should employ the “best practice” for each code. For example, the ALE algorithms would use a moving mesh option, and some Eulerian codes included an extended “buffer” region to damp out unphysical waves reflected from the outflow and inflow boundaries. As detailed previously, simulations were run from t = 0 s to t = 0.5 s that is equivalent to a dimensionless time τ = 6.15, where the basis for non-dimensionalisation is detailed in Sec. II C.

Note that the computation is intended as a test of the high Reynolds number limit of turbulent mixing.^22,23 At early stages where the layer is effectively a highly corrugated, strained, laminar layer, physical diffusion would be important in determining mixing, for example; however, at later times when the flow is fully developed, it is assumed that scalar dissipation is well represented on the chosen grid resolutions. This parallels the research by Dimotakis,²⁴ who examined the “mixing transition,” showing that for Re$>104$⁠, the dissipation rates are insensitive to the actual values of viscosity and diffusivity. Thus, the effects of physical viscosity and diffusivity on the resolved scales are assumed to be zero and all simulations presented here are nominally inviscid.

The grid convergence study demonstrated that several algorithms had integral width W with variations $<1%$ between the lowest and highest grid resolutions on the standard problem. Quantitatively, the change in predicted integral widths at t = 0.5 between the lowest and the highest grid resolutions for each algorithm were Ares (5.8%), Flamenco (0.5%), FLASH (5.1%), Miranda (3.8%), NUT3D (16%), Triclade (0.2%), and TURMOIL (1.6%). This indicated that several of these algorithms could reliably compute a converged integral width for a problem where all physical length scales were divided by four, with the computational domain size maintained, at the highest grid resolution employed in the standard problem. This new case improves statistical fidelity, may be run to much later dimensionless times (⁠$τ≈125$⁠) and thus have a greater chance of achieving self-similarity. In this problem, the initial range of physical length scales in the perturbation was modified to lie between L/32 and L/16 where L is the cross section, and the diffuse layer has an initial thickness L/128. Note that one algorithm (NUT3D) ran a “half-scale” problem, where the initial length-scales were all scaled by a factor of two, instead of four.

To facilitate the use of this test case in future algorithmic development and validation, Fortran and C implementations of the initial conditions may be made available to other research groups by contacting the lead author.

In this subsection, an approach is detailed that can estimate the initial growth rate for a narrowband multimode surface perturbation. This initial growth rate provides a consistent non-dimensionalisation for both problems explored in this paper.

A single mode RMI at low to moderate Mach numbers will have an initial $ȧ=kAt+a0+Δu$⁠, where $a0+=Ca0=(1−Δu/Ui)a0$⁠. The current initial condition consists of a random multimode perturbation represented by a sum of Fourier modes with random phase ϕ,

According to linear theory, after shock passage, the perturbation is given by

where the second expression is valid toward the end of the linear phase, assuming an initial perturbation $a0≪2π/k$⁠. It is assumed that the amplitude compression factor does not depend on the wavenumber $k=ky2+kz2$⁠. The initial variance of the perturbation is given by $σ2(0)=σ02=∑ky,kza02/2=∫0∞P(k)dk$⁠, where P(k) is the power spectrum of the initial perturbation, so that

which implies that $σ̇=k¯At+σ0+Δu$ with $σ0+=Cσ0=C∫0∞P(k)dk$⁠, where $k¯$ is a weighted average wavenumber of the perturbation. For the narrowband perturbation utilised in this paper, P(k) = const. for $kmax/2≤k≤kmax$⁠, giving $k¯=7/12kmax$⁠.

The integral width $W=∫0Lx⟨f1⟩⟨f2⟩dx$ is used to measure the growth of the layer, where $⟨.⟩$ represents a planar averaged in the homogeneous direction and must be related to the variance σ. For the type of perturbation chosen, the height distribution p(A) should be Gaussian,

The dense fluid plane-averaged volume fraction is then

The integral width is therefore

This result is close to that obtained with an assumed linear volume fraction distribution with bubble height equalling spike height h, which gives W = h/3, $σ=h/3$⁠, and thus W = 0.577σ. Thus, the predicted initial rate of increase of the layer width is

For the initial conditions employed within both the standard and quarter scale case, $Ẇ0=12.649$ m/s, from which $W=W0++Ẇ0(t−tshocktime)$ for the linear phase and $t−tshock time≥0$⁠. All non-dimensional quantities within this paper are non-dimensionalised by $Ẇ0$⁠, $λ¯=2π/k¯$⁠, cross-sectional area 4π², and ρ_L = 1. For example, dimensionless time $τ=tẆ0/λ¯$⁠.

This study focuses on the integral properties of the mixing layer including width measures, mixing fractions, total kinetic energy in each direction, and kinetic energy spectra.

The time dependent evolution of integral mixing width W, molecular mixing fraction Θ, and the mixing parameter Ξ are defined as (see, for example, Refs. 8, 11, and 25–27)

where $⟨f1,2⟩$ indicates the y-z plane averaged volume fraction of species 1, 2 where species 1 is the heavy gas. For codes that use mass fractions to track the species, volume fractions are derived from mass fractions assuming local pressure and temperature equilibrium.

Planar averaged kinetic energy in the x, y, and z directions have been computed as the difference of the actual velocities minus the plane averaged velocity, integrated in each cell of volume dV over the entire volume V,

Finally, the two-dimensional variable density spectra taken at the mixing layer centre (defined by $⟨f1⟩=0.5$⁠)^26,28 are

where ν_i = ρ^1/2u_i, $(.)^$ indicates the Fourier transform of a quantity and $(.)^*$ is the complex conjugate of the transform.

Eight independent algorithms have been employed in this study that are briefly described here with references to more complete information.

Ares is an Arbitrary Lagrangian–Eulerian (ALE) code developed at the Lawrence Livermore National Laboratory (LLNL). The Lagrange time step uses a second order predictor-corrector method. The Gauss divergence theorem is used to solve the discrete finite difference compressible multi-component Navier–Stokes equations. Spatial derivatives are approximated using a second-order finite difference method. Artificial viscosity is applied to damp out spurious, high frequency oscillations that are generated near shocks and contact discontinuities. The maximum grid resolution run here is $720×5122$⁠.

Flamenco is a Godunov-type solver that utilises a nominally fifth order reconstruction stencil (in 1D),²⁹ a HLLC (Harten-Lax-van-Leer-Contact) Riemann solver,³⁰ coupled with a low Mach correction,^31–33 and second order Total Variation Diminishing time stepping.³⁴ The governing equations are based on Navier–Stokes plus volume fraction equations.³⁵ The maximum grid resolution run here is $720×5122$⁠.

FLASH is a modular, massively parallel, open source code developed at the University of Chicago to investigate astrophysical flows.³⁶ The “Hydro module” is employed here to solve the compressible Euler equations on a block structured adaptively refined mesh. The Euler equations are solved using the Piecewise Parabolic Method (PPM)³⁷ complemented by the symmetric Monotized Central (MC) limiter.³⁸ The maximum grid resolution run here is $720×5122$⁠.

Hesione is a finite difference Lagrange-remap code. In the Lagrange phase, it uses a second order accurate scheme in space and time as employed in the BBC algorithm.³⁹ For multimaterial flow, the remap phase requires interface reconstruction (IR) based on Young’s method and can operate with mixed meshes.⁴⁰ A key difference in this numerical algorithm is that interface reconstruction inhibits species diffusion, as if the two fluids are immiscible. That means the mixing is heterogeneous (with IR) and homogeneous with all other codes used by the collaboration. The maximum grid resolution run here is $180×1282$⁠.

Miranda uses a tenth-order compact differencing scheme for spatial derivatives and a five-stage, fourth-order Runge–Kutta scheme for temporal integration. Full details of the numerical method are given by Cook.⁴¹ For numerical regularization of the sharp, unresolved gradients in the flow, artificial fluid properties are used to locally damp structures that exist on the length scales of the computational mesh. The maximum grid resolution run here is $720×5122$⁠.

NUT3D solves the Euler equations augmented with mass fraction equations.^42–44 The difference scheme follows the method proposed by Ref. 45. A “sharp” second orderapproximation in space is coupled with the exact solution of the Riemann problem for the fluxes. Time integration is carried out by the predictor-corrector method. The maximum grid resolution run here is $720×5122$⁠.

Triclade solves the Navier–Stokes plus mass fraction equations using a conservative finite difference method and is based on the wave propagation algorithm of LeVeque⁴⁶ with high order accuracy provided by the corrections due to Daru and Tenaud.⁴⁷ A fifth order time-space accuracy is chosen here, referred to as WP5 by Shanmuganathan et al.⁴⁸ and a more detailed description of the method is given in Appendix A.2 of that reference. The maximum grid resolution run here is $1440×10242$⁠.

TURMOIL solves the Euler equations as well as an equation for mass fraction using a Lagrange remap method.^25,49 An un-split, second-order finite difference method is used for the Lagrange phase and the monotonic advection method of van Leer⁵⁰ in the directionally split remap phase. A higher-order artificial viscosity, similar to that proposed by Schultz⁵¹ and Cook, Ulitsky and Miller⁵² is used in these calculations.⁵³ The maximum grid resolution run here is $1440×10242$⁠.

Thus, the eight algorithms include Godunov-type, ALE methods, and compact differences. The algorithms are either Implicit LES (ILES)^25,49,54,55 or Artificial Fluid LES (AFLES).⁴¹ In the case of ILES, it is assumed that there is sufficient separation between the integral length scales and the grid scale such that the growth of the integral length scale is independent of the exact dissipation mechanism. Furthermore, it is assumed that the Reynolds number is sufficiently high that the smallest scales represented on the grid are well mixed. This second assumption is the most restrictive and poses a particular challenge in temporally transitional flows; however, the first assumption has been addressed in the formulation of the problem. Although several of the algorithms have viscous and diffusive terms implemented, the results presented here are from inviscid computations.

Figures 4–6 show the convergence of W, TKX, and Θ for each of the seven codes that provided convergence data. For the standard problem, all algorithms are grid converged with respect to integral width W at the highest grid resolution. Triclade, TURMOIL, and Flamenco are converged with $<$2% variation at $180×1282$ resolution. Miranda, NUT3D, Ares, and FLASH converge at $360×2562$ grid resolution; however, it should be noted that the Miranda and Ares simulations did not use sub-cell sampling to improve the representation of the initial condition on the mesh, which slowed the observed rate of convergence.

As a second measure of the large scales, the kinetic energy converges well for all algorithms; however, most show a moderate increase in early time kinetic energy implying that the solutions are not converged in that region. There is potential for a further increase of up to $≈10$% of the peak kinetic energy at $720×5122$ grid resolution based on Richardson extrapolation of the current parameters. However, it should be noted that Triclade results at $1440×10242$ have $<$2% variation from the $720×5122$ case, lying between the $720×5122$ and $360×2562$ resolutions. This indicates that the kinetic energies are reasonably converged for the very high order algorithms at $720×5122$ resolution.

The mixing measures are challenging to resolve at early times, as the flow is undergoing temporal transition from a sum of linearly growing modes through to a self-similar layer. Thus, at the early linear stages, the measured mixing parameters Θ and Ξ are determined by the ability of the numerical scheme to resolve an extremely stretched contact surface. However, at later times, turbulence increases the amount ofmixing, such that the gradients are again resolved on the given mesh. The target of this problem was to determine the late-time self-similar mixed state of the specified narrowband initial condition, where all modes have saturated and an inertial range established. At the final time, each individual algorithm has converged to within 3%.

For the quarter scale problem, convergence is not as clear. For integral width, based on the standard problem (which is converged at $180×1282$ resolution), Flamenco, FLASH, Triclade, TURMOIL, and Miranda should have converged at the maximum grid resolution. Examining the finest two resolutions (⁠$360×2562$ and $720×5122$⁠), Flamenco and TURMOIL vary by $<$3% for the highest two resolutions; however, for the other algorithms, this difference is larger. This is not possible to demonstrate conclusively without running a higher grid resolution, which is beyond the computational resources available.

For kinetic energies, most algorithms show a worst case variation on the order of 10%–20% between the highest two grid resolutions at the peak post-shock passage value. For mix Θ, each algorithm shows very good convergence, with the greatest difference between the finest two resolutions of 5%.

Visualisations of the standard test case at times τ = 1.23 (t = 0.1 s) and τ = 6.15 (t = 0.5 s) are shown in Figs. 7 and 8 for each algorithm at $720×5122$ and an additional visualisation of the $1440×10242$ Triclade simulation. The visualisations at early time highlight the classical bubble and spike phenomenology associated with RMI at this Atwood number, showing narrow spikes penetrating into the lighter fluid and relatively wide bubbles. The large scale features are qualitatively similar for all algorithms; for example, the prominent spike that is located at $(y,z)≈(5.5,1.3)$⁠. However, there is a clear difference in resolution of the small scales, where the higher order algorithms permit the growth of finer scale fluctuations, as is again highlighted clearly at the aforementioned spike head. At the early time, there are still coherent contact surfaces visible, linking the bubble and spikes in some regions of the flow, e.g., $(y,z)≈(3.5,5.5)$⁠.

The later time visualisations provide an indication of the wider range of length-scales present, implying that the flow is transitioning toward being a fully turbulent mixing layer with a smoother, well mixed, profile. The overall mixing layer width is qualitatively similar between the algorithms, as is the positioning of the key bubbles and spikes. However, there are differences in the fine scale features as highlighted previously, which at later times have caused small differences in the position of these larger features. The computations of the standard problem were terminated at τ = 6.15 (t = 0.5 s) as it is clear that the dominant wavelengths are approaching the domain size, thus impacting the statistical accuracy of the quantities of interest.

Code-averages have been taken of the seven non-interface reconstructing algorithms. Figure 9 plots the integral width for all algorithms, the code-average integral width, and the standard deviation plotted as error bars. The code-averaged results are also compared to three theoretical models, the model developed in Sec. II C for the linear phase, a second inspired by the approach of Mikaelian,⁵⁶ and a third that utilises a straight error function blending of the linear and non-linear stages.

Mikaelian’s model combines growth of the width that is linear in time at early time, matched to a nonlinear growth at late time. This approach is designed to match a post-shock fully turbulent layer width that is linear in time initially with a late time decay through a power law. Here the model is adapted such that the linear phase, given by $W=W0++4αAt+ΔVt$⁠, is calibrated to match the initial linear growth $Ẇ0$ defined in Sec. II C by an appropriate choice of α. The linear phase is then joined to an assumed power law behavior at the temporal location t^* where the growth rates predicted from both models would be equal, giving

where the layer width W = W^b + W^s, θ = θ^b, θ^s, α = α^s, α^b, W^* is the width at the transition between the two power laws, and $W0+$ is the initial layer width. The superscripts (.)^s and (.)^b denote the spike and bubble values that may be allowed to vary. The current curve uses mostly the parameters recommended by Mikaelian,⁵⁶ except for β and α^b,s. Adapting the derivation of Mikaelian,⁵⁶ β governs the point at which the two power laws are matched and is computed as a dimensionless time $β=Δut*σ0≈32$⁠. This is based on an estimation of the matching time t^* being where the code-averaged W varies by 5% from that predicted by the theory in Sec. II C. The second adjustment is in the specification of α^b,s, which has been calibrated to match the initial growth rate provided by the new linear theory. For completeness, the full set of parameters is

Both α^b,s have been set to the same value but could be set to different values as long as the sum remains unchanged. In this way, the coefficients use only quantities that can be estimated a priori. This curve is labelled “Matched Power Laws” in Fig. 9.

The third curve labelled “Blended Power Laws” writes the dimensionless time τ as a function of dimensionless integral width $W̃$ as

where all $(.)̃$ are dimensionless quantities and $W̃0+$ is the dimensionless post-shock integral width from theory. On the right-hand side of this equation, $τshocktime$ takes into account the time taken for the shock to reach the interface (0.0011 s) and a delay in initial growth due to the inversion process (⁠$Δτ≈0.015$⁠). The second term gives a linear growth $W̃=τ$ at early time, where the error function removes this contribution at later time. The third term blends in the zeroth order term τ₀ in the data fit for the late time nonlinear behavior in the form $W̃=Ãτ−τ0θ$⁠, with the linear initial growth rate. Finally, the fourth term is responsible for correctly representing the late time power law behavior.

The dimensionless parameters in the non-linear fit are

These have been determined from the quarter scale problem, as will be discussed in Sec. IV; $B=0.395$ and $C=5.8$ are dimensionless free parameters set “by eye” as non-linear regression in Mathematica® that did not perform well in optimising these coefficients. Assuming that the dimensionless parameters $B$ and $C$ and exponent θ are the same for all narrowband cases, all other quantities can be evaluated a priori.

Examining Fig. 9 and the close up shown in Fig. 10, the new model predicting the linear growth rate of a multimode narrowband perturbation shows excellent agreement with the code-averaged data, and the Flamenco data that are plotted as the sampling rate are much higher than the code-averaged data in the linear phase. Regarding the “blended” and “matched” power laws, both of the models show very good agreement within their regions of applicability and underlying assumptions. The “Matched” power laws overestimate late time behavior as the θ prescribed is higher than that computed from the simulations. The prediction from the “Blended Power Laws” is better; however, it does not match as well as the non-linear regression since this formula is calibrated to match the very late time behavior simulated in the subsequent quarter scale case. For future modeling and code validation, the dimensional data are tabulated in Table I.

There is excellent agreement between all algorithms, which is expected given that this is a measure of the large scales and should be relatively algorithm-independent. The standard deviation of the code-averaged integral width is $σW≤0.7$%. As expected, the interface reconstruction code Hesione differs from the diffuse-interface algorithms as the interface reconstruction method suppresses the dissipation of density fluctuations. The predicted integral width for Hesione is 3.7% higher at the latest time when compared with Ares at the same grid resolution, which uses a similar method but without IR.

Two approaches have been employed to compute the simulated value of θ: the first is to employ non-linear regression to fit the growth curve for $τ>1.23$ (⁠$t>0.1$ s) to a form W = A(t − t₀)^θ, and the second is to utilise the derivatives of the integral width to directly determine $θ−1=1−ẄW/Ẇ2$⁠. This power law behavior of the mixing zone width (⁠$W∝tθ$⁠) is a consequence of the interaction of large scale structures in the mixing zone (low-k part of the spectrum) as proved by Poujade and Peybernes.⁵⁷ Using a simple two buoyancy-drag model and equating fluid inertia with drag (see, for example, Ref. 58), the Richtmyer–Meshkov mixing layer width will grow in the self-similar regime according to

which yields the power law solution for RMI of θ⁻¹ = 1 + C_d.

Using non-linear regression on the code-averaged integral width for $τ>1.23$⁠, assuming W = A(t − t₀)^θ gives

in line with previous results for early time growth (not yet self similar) of the mixing layer.^9,11,12 This is plotted in Fig. 9 as the solid black line that passes through all data points. Note that θ increases at later time, for example, at $t>0.4$⁠, θ = 0.228. Table II tabulates the values of θ for each individual code, illustrating that the range $0.210<θ<0.228$ brackets all diffuse interface algorithms.

The time varying value of θ determined from the derivatives of W is plotted in Fig. 11, where the derivatives were evaluated using second-order central differences with temporal spacing of Δt = 0.04 s. This shows that following the initial linear growth for $τ<0.8$⁠, the value of θ moderately decreases from $τ>1.23$ to $τ<4.8$⁠, followed by an increase in late time. A time varying θ has been observed previously by Thornber et al.,^11,17 Oggian et al.,¹⁴ and Lombardini et al.¹²(although not explicitly discussed therein). The structure of the layer is continually evolving throughout the simulation time from inviscid linear growth through to a developing turbulent layer, where Barenblatt et al.⁵⁹ postulated a narrow-band growth exponent of $θ≈2/3−ν$ with ν being a viscous correction that will naturally evolve in time due to the changing nature of the mixing layer.

The integral mix measures Θ and Ξ are plotted in Fig. 12, including the code-averaged data (not including Hesione) and the standard deviation is plotted as error bars. The data are also tabulated in Table I. For Hesione, interface reconstruction (IR) inhibits numerical diffusion so that the two fluids are treated as immiscible. That means the mixing is heterogeneous, which leads to a predicted $Θ≈0.24$ if mixed cells are homogenised, in contrast to 0.7–0.8 for other algorithms. The α collaboration¹⁹ reported the same difference using HYDRA and ALEGRA in their interface reconstruction (IR) version.The computation of mixing properties for the IR algorithms is a useful limit, effectively describing the variation of the interface area with time when fluid mixing is inhibited. Clearly, the IR algorithm does not represent well the growth of miscible layers; however, it is more representative of the immiscible limit.

Referring specifically to the diffuse interface algorithms, there is substantially greater variation in integral mix measures compared with the integral width. The early time quasi-linear phase shows good agreement between all codes, where the initialised diffuse interface is thinned by the extension of the contact surface as the fluids inter-penetrate. As time progresses, the energy containing scales begin to break down, transferring energy to smaller scales, which then increase the efficiency of the mixing by stirring the fluid, thus increasing the cross-sectional area.

Employing the time-to-transition criteria of Zhou⁶⁰ and taking δ = L/4 and $u=2/3(TKX+TKY+TKZ)$ evaluated at t = 0.03 s implies that the transition time is $τtr≈0.18$ (⁠$ttr≈0.015$ s). However, at such an early stage, the components of kinetic energy are strongly biased toward the shock-direction (x); using only the TKY and TKZ components gives $τtr≈0.42$ (⁠$ttr≈0.034$ s). This is approximately the time of minimum mix parameters, i.e., where mixing begins to accelerate. The time of minimum mixing measures is not very sensitive to algorithm choice. As time increases, the layer tends toward a fully established mixing layer in which the mix measures increase. Should perfect self-similarity be achieved, these integral mix measures should tend toward a constant value.

From the mix parameters, it is clear that (i) the layer has not achieved self-similarity in the prescribed period and (ii) uncertainty in mix is greatest at the transition period. This can be seen clearly in the error-bars plotted for the code-averaged results, where the maximum standard deviation is σ_Θ = 12% at τ = 0.37 (t = 0.03 s), moves below 10% at τ = 1.1 (t = 0.09 s), and is 2.8% at τ = 6.15 (t = 0.5 s). This uncertainty in the results from differing algorithms is determined by two numerical factors. First, at an early time, the contact surface is stretched and becomes thinner than the local grid scale. Thus, the computed mix measures are dependent on the minimum thickness contact surface that each algorithm is able to represent, since once the contact surface is smaller than that, the layer can no longer thin. The second factor is that each algorithm has a different minimum number of points required to resolve the growth of a single vortex, which again is dependent on the dissipative properties of the algorithm (strongly linked to the order of accuracy). The mixing layer must be resolved with a sufficient number of points to be able to simulate the physically present range of the length scales in the problem. It should be noted that in the high Re limit and given sufficient grid resolution, the species variance should be sufficiently resolved using a Large-Eddy Simulation, and thus all algorithms should converge.

Combining the complicated physical and numerical factors, it is encouraging that the standard deviation of the code-averaged results reduces considerably at late time (see Table I), indicating that although self-similarity has not yet been achieved, the simulations appear to be tending to a much narrower uncertainty band at late times. At the latest time, $Θ=0.769±0.0225$ and $Ξ=0.767±0.0229$⁠, which is consistent with previous studies of narrowband and other closely related perturbations,^{11,12,14,17,22} and close to the experimentally measured results of Krivets et al.⁶¹ The predicted mixing is lower than those by Orlicz et al.⁶² that employed a substantially different, strongly two-dimensional (at the large scale) perturbation of a gas curtain. It is expected that mixing would be higher in that experiment since, aside from the large differences in initial perturbation, the gas curtain can mix on both sides of the interface. Weber et al.⁶³ measured Ξ = 0.85 at the latest time; however, they noted that the mixing parameter was decreasing at late times.

Table II gives the values of Ξ and Θ for each algorithm at the latest time. Where convergence is uniform, asymptotic values have been computed using Richardson extrapolation of the three provided grid levels. Note that the Richardson extrapolated values are not expected to be realistic (particularly for Ares and Flamenco), even though their convergence was relatively uniform, hence their latest time value is also given in parentheses. The time of minimum mixing measures is not very sensitive to algorithm choice.

Figure 13 plots the fluctuating kinetic energy in each component direction and a simple measure of anisotropy 2TKX/(TKY + TKZ) as a function of time, where in the plots all kinetic energies (and spectra) are non-dimensionalised by $ρLẆ02λ¯A/2$⁠. Figure 14 plots the average of all of the algorithms except the two outliers, TURMOIL and Hesione. A sample of the same data is also tabulated in Table III.

As kinetic energy is primarily a measure of the large scales, it should converge well with increasing grid resolution. At early time, there is a strong peak in fluctuating kinetic energy in all directions, as the measure also includes contributions from the non-planar shock. As the shock clears the mixing layer and straightens, the kinetic energy deposition by the shock wave dominates and peaks at τ = 0.12 (t = 10⁻² s) in the shock direction, and τ = 0.62 (⁠$t≈5×10−2$ s) in the in-plane directions. This delay is due to the transfer of energy from the shock-direction component into the in-plane components and is visible in the anisotropy measure. This time scale matches well with the time-to-transition estimate of $τtr≈0.42$ (⁠$ttr≈0.034$ s) in Sec. III B. The kinetic energy components then decay with an approximate power-law form, where an additional line is plotted in Fig. 13 representing the predicted self-similar decay rate based on $W∝tθ$⁠, which gives fluctuating kinetic energy $∝t3θ−2$ (see, for example, Ref. 11).

Overall, six of the algorithms collapse nearly perfectly, which is an excellent result given that each algorithm has quite a different order of accuracy and dissipation mechanism. This strongly suggests that a reasonable separation of energetic scales from dissipative scales has been achieved. There are two algorithms that give different solutions, Hesione and TURMOIL, which have consistently higher kinetic energy throughout the layer, the reason for which has not yet been identified but may be linked to the treatment of acoustic waves in the ALE method. The code averages have been taken of the six algorithms, neglecting these two outliers, as otherwise the standard deviations are effectively dominated by just the single algorithm. Taking this approach, the standard deviations are remarkably small, $≪$1% at early time and $≈5$% at the latest time.

Plotting both the anisotropy for each algorithm and the code averaged kinetic energy anisotropy [2TKX/(TKY + TKZ)], it can be seen that for this initial condition and time period, the results are not tending toward isotropy. At the final time, this anisotropy measure is 1.43 for the code-averaged data, where the anisotropy predicted by each individual algorithm ranges from 1.25 (TURMOIL) to 1.55 (Miranda). This lends substantial additional support to prior observations that although the initial anisotropy reduces due to the transfer from the shock direction to the in-plane directions, this transfer does not result in isotropy prior to the commencement of the power-law decay.^11,12,14,64

The variable density kinetic energy spectra taken at a single plane at the centre of the mixing layer defined as $⟨f1⟩=0.5$ are shown in Fig. 15 at $τ>1.23$ (⁠$t>0.1$ s) and τ = 6.15 (t = 0.5 s) for Ares, Flamenco, FLASH, Miranda, and Triclade at $720×5122$ and Triclade at $1440×10242$⁠. The wavenumbers have been normalised by $k¯=6.11$⁠. The spectra at τ = 1.23 peak at k = 3–5, then for most algorithms show the early stages of a constant power-law region to $k/k¯≈16$⁠, dependent on algorithm. This behavior is mirrored in the τ = 6.15 results, where the peak is now located at $k/k¯=0.33$–0.5. The constant power-law region has a slope of $≈k−1.3$ for $k/k¯=3$–5, which is shallower than a classical Kolmogorov spectrum or that proposed by Zhou.⁶⁵ 

Comparing the algorithms, the agreement at the large scales is remarkably good considering that this is a spectrumfrom a single plane located at the centre of the mixing layer. There are differences in the dissipative properties of the schemes, as expected from the different orders of accuracy and algorithmic constructions. The spectra from FLASH depart from the bulk of the results at $k/k¯≈1.3$–1.6, most likely due to numerical dissipation. Ares has a higher kinetic energy in the range k = 0.5–8.2, which may be a numerical manifestation of the physically observed “bottleneck” phenomenon and then damps the high wavenumber energy to zero.

The second-order TURMOIL scheme also shows significant fine-scale structure. The ability of TURMOIL to resolve fine scale turbulent features has been documented elsewhere^{11,48,53,66,67} and is a result of a number of factors. In particular, the second-order Lagrange phase of the TURMOIL algorithm is close to symplectic in flows where artificial viscosity is small; the remap phase is third order in each swept direction and can use less conservative assumptions when reconstructing gradients than a full Godunov scheme, as it is performing a purely passive advection. Both phases are independent of the Mach number in the low-Mach limit. For flows subject to non-linear perturbations, the overall fidelity of a numerical scheme may be more important than its formal order of accuracy for linear perturbations.

Triclade, Flamenco, and Miranda agree very well throughout the power-law region, with the key difference between the three in the approach to the grid cutoff. Relative to the Triclade simulation at $1440×10242$⁠, the Miranda and Triclade simulations have a small “bump” at the end of the power-law region, then uniformly damp fluctuations to zero as the cutoff is approached. Miranda has the sharpest damping at high wavenumbers. The spectrum in the Flamenco simulation does not damp to zero at the high wavenumbers. This high wavenumber behavior is most likely responsible for the observed differences in the integral mix measures in Sec. III B, as the schemes with the highest energy at high wavenumbers have the highest values of Θ and Ξ (see Fig. 12).

Based on the observation that the key integral quantities were grid converged at $180×1282$ resolution for several algorithms, a second simulation was undertaken at $720×5122$ resolution where all initial physical length scales were divided by four, with the computational domain size fixed. This permitted computations to more than eight times the dimensionless time (⁠$tȧ/λmin≈380$–760). Due to the computational expense, each simulation was run for the longest time possible. Most of the simulations terminated before τ = 123 (t = 2.5 s), which is when parts of the leading spikes exit the domain. Some simulations were continued as far as τ = 246 (t = 5 s). Table IV details the expected maximum error in the integral widthestimated from the standard problem, the final time each algorithm was run to, and key quantities recorded or computedfrom the results. In this section, all results are presented together.

A visualisation of the resultant flow field is shown at four times from the Flamenco code in Fig. 16. The earliest time is at the same dimensionless time as the first time instant of the standard problem shown in Fig. 7. Here the contact surface between the two fluids is relatively well defined, although individual spikes and bubbles are breaking down. At the later times, the mixing layer develops through transition to the start of an established mixing layer.

Figure 17 (left) plots the integral width for all algorithms, and Table IV summarises the data fits for the integral width for each individual algorithm. The initial conditions were designed to be at the limit of expected convergence. As such, it is important to clarify the expected confidence in the solutions. Column one of Table IV details the expected error based on the error between the $180×1282$ and $720×5122$ simulations on the standard problem (for NUT3D this was based on the error from $360×2562$ to $720×5122$⁠). The results in Fig. 17 show that there is a grouping of algorithms whose time-dependent integral width agrees relatively closely, including Flamenco, FLASH, NUT3D, Triclade, and TURMOIL. This is consistent with the predicted convergence error. The standard case demonstrates that the Ares simulation should converge down toward the other algorithms and that the Miranda simulation should converge upwards if a higher grid resolution was run. This behavior has been confirmed by shorter simulations at higher resolution at early time, however, not for the full simulation for which resources were not available.

The integral mixing measures vary by $<$5% following τ = 24.6 (t = 0.5 s), thus non-linear regression was employed again to determine θ, A, and t₀ for each algorithm. The second column of Table IV details the latest time to which the problem was simulated, which is the latest time of the data fit if the algorithm was not run past τ = 123 (t = 2.5 s). If it was run to later time, the computed data fit used only data up to τ = 123 (t = 2.5 s). The code-average was not employed due to the larger scatter relative to the standard problem. Table IV plots these values and their means. The mean values were $θ=0.291±0.025$ where all algorithms were considered; however, if Miranda and Ares results are excluded, then $θ=0.292±0.009$⁠. The computed values of t₀ and A are equally consistent. Note that this value of θ is quite close to the value consistent with the self-similar solution of the multicomponent k–epsilon turbulence model, which gives θ = 0.3 using standard coefficient values based on shear turbulence and predicts very good agreement with many experimentally measured pre-reshock mixing layer widths in numerically modeled Richtmyer–Meshkov instability.^68,69

Figure 18 plots the integral width for the quarter scale problem with the standard problem and an additional simulation run with a sharp interface with different random numbers.¹⁷ The integral width for the quarter scale problem varies by $<$0.1% from the case with a different set of random numbers (and random number generator) and sharp interface, thus indicating that there is (i) little residual dependence on the phase relation of the initial conditions, (ii) little residual dependence on statistical variations in the individual modal amplitudes, and (iii) a small impact of the diffuse layer on the vorticity deposited in the layer. The results from the standard case collapse on top of the quarter scale case under the chosen non-dimensionalisations. The “blended” and “matched” power-law models are plotted, and the agreement with the dataset is good, with the matched power laws showing a 13.7% overestimation at the latest time and blended within 1.2% as expected, as it is based on the non-linear regression data fit.

The Flamenco simulations were continued to τ = 246 to explore the late time behavior of the mixing layer; however, note that past τ = 123, spikes are observed exiting the domain. This extended data set permitted the plotting of θ(t), as shown in Fig. 17. At late time, there is a substantial reduction in θ; however, this is also mirrored in the non-linear regression, which as a function of the starting time gives θ = 0.285 $(τ>50)$⁠, θ = 0.274 $(τ>100)$⁠, θ = 0.263 $(τ>147)$⁠, θ = 0.251 $(τ>197)$⁠, for example.

Jumps in θ can be seen due to small oscillations in the mixing layer width (not visible to the naked eye) that the computed value of θ is very sensitive to (visible in variations of $Ẅ$ but not of $Ẇ$⁠). This strong sensitivity has been noted before, for example, Fig. 4 in Oggian et al.¹⁴ where the authors further filtered the width data to smooth these small oscillations. In the Flamenco computations, the spikes in θ are caused by very small pressure waves generated at the computational domain boundary traversing. The first at $τ≈8$ is caused by a weak reflected pressure wave generated by the shock wave exiting the computational domain (amplitude 0.1% of the post-shock pressure, $≈400$ Pa). The latter spikes in θ, for example, at $τ≈1$⁠,are caused by weak pressure waves (less than 1 Pa) generated as particularly strong (in the context of the decaying layer) vortical structures approach and exit the computational domain.

Next, the values of θ have been computed using two alternative mixing measures, the mixed mass $M$⁠,⁷⁰ and the layer width h:^16,71

where $⟨.⟩$ indicates a plane-averaged quantity and Y_i denotes the mass fraction of gas i. They both give a power-law dependence similar in form to the integral width. Using the Flamenco data and non-linear regression, the $θ(τ>24.6)$ based on these quantities gives 0.292 for h and 0.287 for $M$ indicating a relative insensitivity of the computation of θ using these different mixing measures (⁠$θ=0.291±0.005$⁠).

Figure 19 plots the integral mixing measures, Θ and Ξ, and shows a comparison of the results of Ξ and Θ against a recently proposed normalised mixed mass measure:⁷⁰ 

The integral mix measures Θ and Ξ at the final time, or τ = 123, are detailed in Table IV for each algorithm. The mean of both are

showing excellent inter-algorithm agreement of these late time, quasi-self-similar statistics. As with the standard case, these are consistent with previous results^{11,12,14,17,22,63} and close to experimentally measured results of Krivets et al.⁶¹ but lower than those of Orlicz et al.⁶² The maximum difference between the two at the final time for all algorithms is 1.6%. All algorithms but Ares show decreasing mix measures as time progresses; however, the maximum variation for $149≤τ≤123$ is $<$1% and for $374≤t≤748$ Miranda and Flamenco show variations of $<$0.4%. Normalised mixed mass Ψ mirrors the behavior of Θ, lying between 0.987Θ and 0.982Θ for the entire simulation and asymptoting to 0.795 and 0.769 for Flamenco and Triclade, respectively, at the latest time.

The turbulent kinetic energy components are plotted in Fig. 20, along with the previously defined anisotropy measure (also detailed at the final time in Table IV). There is a greater scatter between codes as a result of the marginal grid convergence in this test case. All algorithms show power-law behavior in time for the kinetic energy components (TKZ is omitted since it is negligibly different from TKY); however, the exponents differ by small amounts. Plotted on Fig. 20 is the power-law t^3θ−2 with the code-averaged θ = 0.291 that shows reasonable agreement with Flamenco, Miranda, and Triclade at late time $τ>24.6$⁠.

The observed anisotropy in this case increases following the minimum observed for the standard case at the latest time (Fig. 13), then ranges from 1.49 to 1.66 over all algorithms. From this result, it can be concluded that the decrease in anisotropy at “late” times in the standard case was a transient phenomena, as a minimum is reached in the quarter scale problem at $τ≈6$ that is close to the final time of the standard problem. This recovery can also be seen in the results of Oggian et al. albeit at a shorter dimensionless time.¹⁴ There is no return to isotropy, and the anisotropy is maintained at anearly constant value after the initial transients within the simulation time scale.

The ratio $⟨λb⟩/⟨hb⟩$ of the dominant bubble wavelength to the bubble front amplitude is a measure of self-similarity in the flow. For self-similar flows, $⟨λb⟩/⟨hb⟩$ saturates to a constant value but may still depend weakly on other factors such as the initial conditions or the Atwood number of the flow. The self-similarity ratio is also an input parameter in buoyancy-drag models,^72,73 where it is used to evolve the lateral length scales beyond the point of nonlinear transition in the flow. Figure 21(a) is a late-time image (τ = 6.2) of the bubble front from FLASH simulations of the standard problem. The bubble front is defined as the streamwise location of the 1% volume fraction iso-surface at any given (y,z) point. The autocorrelation procedure developed in Refs. 74 and 75 is adopted to extract the dominant bubble wavelengths from such images. Briefly, the autocorrelation function η(y, z) of the bubble front Z_b(y, z) is computed according to

The dominant bubble diameter $⟨Db⟩$ is then obtained as the radial length, where the azimuthally averaged function $⟨η⟩b$ drops below a threshold value of 0.3 (calibrated from severaltest images). The corresponding bubble wavelength is then obtained from the relation $⟨λb⟩≈⟨Db⟩(ρh+ρL)/ρh$⁠, where the post-shock densities are used.⁷⁶ 

The results of this analysis are shown in Figs. 21(b) and 21(c) for the standard problem and the quarter-scale problem, respectively. For the standard problem, the analysis was performed for FLASH simulation data and shows a late-time τ = 3.69 (⁠$t>0.3$ s) saturation to a nearly constant value of $0.85±0.01$⁠. This time-to-saturation may reasonably be taken as the time required by the dominant wavelength in the initial wavepacket to reach nonlinearity, thus triggering the onset of self-similarity for the flow. For the quarter-scale problem in Fig. 21(c), the analysis was performed on FLASH and Triclade data, but they both obtain very similar saturated values of $0.85±0.07$ and $0.86±0.02$⁠, respectively. Furthermore, since the dominant wavelength from the initial condition is smaller for the quarter-scale problem, the time required by that mode to reach nonlinearity is also shorter as expected from the $W∝1/λmin$ scaling.

Finally, $⟨f1⟩$ and $⟨f1f2⟩$ are plotted against scaled distance x/W for several times for Flamenco and Triclade in Fig. 22. By scaling the spatial dimension by the layer width, this figure illustrates the validity of the assumption that the mixing layer profiles can be collapsed via a single length scale. Should the bubbles and spikes evolve with a differing power law (different θ^b and θ^s), there should be a noticeable lack of collapse under a single scaling. Examining the profiles of $⟨f1⟩$ shows that the agreement between Flamenco and Triclade is very good and that both are showing a slight lack of self-similarity, most notably a smoothing of the bubble side (close to $⟨f1⟩=1$⁠) at later τ, as shown in the inset of Fig. 22 (left). Examining the $⟨f1f2⟩$ profiles shows that the extremes of the spike side (⁠$x>0$⁠), and to a lesser extent, the extremes of the bubble side (⁠$x<0$⁠) have not reached self-similarity. Despite this, the core of the mixing layer $−4≤x/W≤4$ is reasonably collapsed at all times indicating that a single length scale and, equivalently a single θ, is a reasonable approximation.

A carefully designed narrowband multimode Richtmyer–Meshkov diffuse interface test case was proposed, and results were obtained for eight independent algorithms representing a broad range of numerical approaches. A new analytical estimation of the initial linear growth of a narrowband multimode perturbation was also presented.

A detailed study of the integral measures for a standard case for which grid convergence can be demonstrated has been presented, showing excellent agreement for all algorithms at the highest grid resolution. Code-averages and standard deviations were then computed, enabling an uncertainty estimate not possible without multiple independent simulations. This standard case demonstrated that the transitional RMI problem has width W = 0.807(t − 0.0309)^0.219 (⁠$0.210<θ<0.228$⁠) for diffuse interface algorithms, mix measures $Θ=0.769±0.0225$⁠, $Ξ=0.767±0.0229$⁠, and a power-law decay of fluctuating kinetic energy approximately proportional to t^3θ−2. It is highly anisotropic, where the shock direction kinetic energy is 43% greater than each in-plane component. Spectra demonstrated a short power-law range of the form $≈k−1.3$⁠, shallower than expected from previous results, although the mixing layer had not achieved self-similarity.

Based on these results, a second problem was run by a subset of the codes with one quarter the physical length scales but the same computational domain size, enabling the simulations to reach more than eight times the dimensionless time. From the converged algorithms, the growth exponent was $θ=0.292±0.009$ computed using non-linear regression; however, the exponent was decreasing with time. This value of θ was insensitive to the choice of width measure when comparing against h or normalised mixed mass $M$⁠. The integral mix measures were $Θ=0.792±0.014$ and $Ξ=0.800±0.014$ and, whilst not constant during the simulation, vary by $<$2.5% in the period during which the growth exponent is fit. The fluctuating kinetic energy showed greater uncertainty in this case; however, all algorithms clearly showed a minimum anisotropy at early times, followed by a marginal increase in anisotropy at later times to $1.49≤2TKX/(TKY+TKZ)≤1.66$⁠.

The analytical estimation of the linear growth rate of a multimode narrowband perturbation was found to be in excellent agreement with numerical data, and two models of the growth rate are presented and validated, which merge linear and non-linear growth relations into a single expression.

Thus, this study has quantified the value of θ for a narrowband mixing layer, including uncertainty bands based on an ensemble of simulations with state-of-the-art algorithms. The RMI layer is predicted to have late time anisotropy for a wide range of algorithms for the time scales simulated. The value of θ reflects a likely lower-bound on the growth rate due to a multimode pure RM instability in the planar case, whereas the mixing parameters reflect an upper bound. Plotting the plane averaged volume fraction profiles against x/W indicates that the bulk of the layer collapses self-similarly through a single length-scale, implying that there is a single θ or that $θb≈θs$ for this problem. Furthermore, the proposed case is now sufficiently documented and analysed to benchmark future algorithms for shock-induced turbulent mixing.

An important conclusion of this study is that with current computational power, it is very difficult to converge the integral properties over a number of algorithms. Thus, when examining higher order moments such as probability distribution functions and spectra, great care must be taken to ensure that they are well converged such that the underlying physics can be correctly elucidated.

The authors would like to strongly acknowledge Dr. Karnig Mikaelian who suggested this collaboration at the International Workshop on the Physics of Compressible Turbulent Mixing in 2014. This research was supported under Australian Research Council’s Discovery Projects funding scheme (Project No. DP150101059). The author would like to acknowledge the computational resources at the National Computational Infrastructure through the National Computational Merit Allocation Scheme that were employed for all cases presented here. FLASH was developed by the DOE-sponsored ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. P. Ramaprabhu was partially supported by the U.S. Department of Energy (DOE) under Contract No. DE-AC52-06NA2-5396. Results from TURMOIL are ©British Crown Owned Copyright 2017/AWE. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA2734.

The perturbed interface A(y, z) is determined according to a specified power spectrum and standard deviation following previous studies where full details may be found;^9,11 here an outline will be given. A perturbation power spectrum P(k) is assumed

where $k=ky2+kz2$ is the wave vector of the perturbation. First, the power spectrum is rewritten as a two-dimensional power spectrum in wave space