Turbulent Richtmyer–Meshkov instability (RMI) is investigated through a series of high resolution three-dimensional simulations of two initial conditions with eight independent codes. The simulations are initialised with a narrowband perturbation such that instability growth is due to non-linear coupling/backscatter from the energetic modes, thus generating the lowest expected growth rate from a pure RMI. By independently assessing the results from each algorithm and computing ensemble averages of multiple algorithms, the results allow a quantification of key flow properties as well as the uncertainty due to differing numerical approaches. A new analytical model predicting the initial layer growth for a multimode narrowband perturbation is presented, along with two models for the linear and non-linear regimes combined. Overall, the growth rate exponent is determined as θ=0.292±0.009, in good agreement with prior studies; however, the exponent is decaying slowly in time. Also, θ is shown to be relatively insensitive to the choice of mixing layer width measurements. The asymptotic integral molecular mixing measures Θ=0.792±0.014, Ξ=0.800±0.014, and Ψ=0.782±0.013 are lower than some experimental measurements but within the range of prior numerical studies. The flow field is shown to be persistently anisotropic for all algorithms, at the latest time having between 49% and 66% higher kinetic energy in the shock parallel direction compared to perpendicular and does not show any return to isotropy. The plane averaged volume fraction profiles at different time instants collapse reasonably well when scaled by the integral width, implying that the layer can be described by a single length scale and thus a single θ. Quantitative data given for both ensemble averages and individual algorithms provide useful benchmark results for future research.

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,1 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.3 

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,6 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,3 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 collaboration19 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 studies11,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.

FIG. 1.

Schematic of the initial condition (not to scale).

FIG. 1.

Schematic of the initial condition (not to scale).

Close modal

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

0.0<x<3.0(ρ,u,p)=6.37497kg/m3,61.48754m/s,399995.9Pa,
(1)
initial heavy(ρ,u,p)=(3.0kg/m3,291.575m/s,100000Pa),
(2)
initial light(ρ,u,p)=(1.0kg/m3,291.575m/s,100000Pa).
(3)

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/m3, 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/m3 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 Ui 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 f1 and f2 = 1 − f1 within the mixed cells. Here, a diffuse interface of gradient thickness of δ = L/32 is employed, thus f1 is computed as

f1(y,z)=12  erfcπ[xS(y,z)]δ,
(4)

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 f1 = 1 − f2 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 f1 = 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 t0.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.17 As an additional check, the normalised two-point, in-plane longitudinal correlation,

Rvv+ww(s,t)=12v(x,y,z,t)v(x,y+s,z,t)v(x,y,z,t)v(x,y,z,t)+w(x,y,z,t)w(x,y,z+s,t)w(x,y,z,t)w(x,y,z,t),
(5)

is computed for 0sL/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.

FIG. 2.

Iso-surface of volume fraction f1 = 0.5 at 180×1282 at τ = 0.

FIG. 2.

Iso-surface of volume fraction f1 = 0.5 at 180×1282 at τ = 0.

Close modal
FIG. 3.

Two-point longitudinal velocity correlations plotted for various τ from the Flamenco computation at 720×5122 resolution.

FIG. 3.

Two-point longitudinal velocity correlations plotted for various τ from the Flamenco computation at 720×5122 resolution.

Close modal

The computational domain is Cartesian as shown in Fig. 1 and measures x×y×z=2.8π×2π×2π m3. 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,24 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 ϕ,

A(y,z)=ky,kza0cos(kyy+kzz+ϕ).
(6)

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

A(y,z)=ky,kzCa0(1+kAt+Δut)cos(kyy+kzz+ϕ)
(7)
ky,kzCa0kAt+Δutcos(kyy+kzz+ϕ),
(8)

where the second expression is valid toward the end of the linear phase, assuming an initial perturbation a02π/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=0P(k)dk, where P(k) is the power spectrum of the initial perturbation, so that

σ2(t)k12Ca0kAt+Δut2CAt+Δut20k2P(k)dk,
(9)
k¯=0k2P(k)dk0P(k)dk,
(10)

which implies that σ̇=k¯At+σ0+Δu with σ0+=Cσ0=C0P(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/2kkmax, giving k¯=7/12kmax.

The integral width W=0Lxf1f2dx 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,

p(A)=1σ2πexpA22σ2.
(11)

The dense fluid plane-averaged volume fraction is then

f¯1(x)=prob(A>x)=xp(A)dA=0p(A)dA0xp(A)dA=121erfx2σ.
(12)

The integral width is therefore

W=141erf2x2σdx=24σ1erf2(t)dt,
(13)
0.564σ.
(14)

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

Ẇ0=0.564712kmaxAt+σ+Δu.
(15)

For the initial conditions employed within both the standard and quarter scale case, Ẇ0=12.649 m/s, from which W=W0++Ẇ0(ttshocktime) for the linear phase and ttshock time0. All non-dimensional quantities within this paper are non-dimensionalised by Ẇ0, λ¯=2π/k¯, cross-sectional area 4π2, 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)

W(t)=0Lxf1f2dx,
(16)
Θ(t)=f1f2dxf1f2dx,Ξ(t)=min(f1,f2)dxmin(f1,f2)dx,
(17)

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,

TKX=xyz12ρ(uũ)2dV,   ũ=yzρudVyzρdV,
(18)
TKY=xyz12ρv2dV,
(19)
TKZ=xyz12ρw2dV.
(20)

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

E(k)=12ν^i*ν^i,
(21)

where νi = ρ1/2ui, (.)^ 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),29 a HLLC (Harten-Lax-van-Leer-Contact) Riemann solver,30 coupled with a low Mach correction,31–33 and second order Total Variation Diminishing time stepping.34 The governing equations are based on Navier–Stokes plus volume fraction equations.35 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.36 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)37 complemented by the symmetric Monotized Central (MC) limiter.38 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.39 For multimaterial flow, the remap phase requires interface reconstruction (IR) based on Young’s method and can operate with mixed meshes.40 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.41 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 LeVeque46 with high order accuracy provided by the corrections due to Daru and Tenaud.47 A fifth order time-space accuracy is chosen here, referred to as WP5 by Shanmuganathan et al.48 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 Leer50 in the directionally split remap phase. A higher-order artificial viscosity, similar to that proposed by Schultz51 and Cook, Ulitsky and Miller52 is used in these calculations.53 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).41 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.

FIG. 4.

Convergence of W.

FIG. 4.

Convergence of W.

Close modal
FIG. 5.

Convergence of TKX.

FIG. 5.

Convergence of TKX.

Close modal
FIG. 6.

Convergence of mix measure Θ.

FIG. 6.

Convergence of mix measure Θ.

Close modal

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).

FIG. 7.

Visualisation of heavy gas volume fraction f1 at τ = 1.23 (t = 0.1 s) for the standard problem. (a) Flamenco, (b) NUT3D, (c) Triclade, (d) TURMOIL, (e) Ares, (f) Miranda, (g) FLASH, (h) Triclade 1440 × 10242.

FIG. 7.

Visualisation of heavy gas volume fraction f1 at τ = 1.23 (t = 0.1 s) for the standard problem. (a) Flamenco, (b) NUT3D, (c) Triclade, (d) TURMOIL, (e) Ares, (f) Miranda, (g) FLASH, (h) Triclade 1440 × 10242.

Close modal
FIG. 8.

Visualisation of heavy gas volume fraction f1 at τ = 6.15 (t = 0.5 s) for the standard problem. (a) Flamenco, (b) NUT3D, (c) Triclade, (d) TURMOIL, (e) Ares, (f) Miranda, (g) FLASH, (h) Triclade 1440 × 10242.

FIG. 8.

Visualisation of heavy gas volume fraction f1 at τ = 6.15 (t = 0.5 s) for the standard problem. (a) Flamenco, (b) NUT3D, (c) Triclade, (d) TURMOIL, (e) Ares, (f) Miranda, (g) FLASH, (h) Triclade 1440 × 10242.

Close modal

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,56 and a third that utilises a straight error function blending of the linear and non-linear stages.

FIG. 9.

Integral width W for all algorithms at the highest grid resolution and the code-averaged W compared to the new linear model [Eq. (15)], the matched power law [Eq. (22)], the blended power law [Eq. (25)], and the non-linear regression data fit for the standard problem (solid line passing through the symbols, starting at τ = 1.23).

FIG. 9.

Integral width W for all algorithms at the highest grid resolution and the code-averaged W compared to the new linear model [Eq. (15)], the matched power law [Eq. (22)], the blended power law [Eq. (25)], and the non-linear regression data fit for the standard problem (solid line passing through the symbols, starting at τ = 1.23).

Close modal

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

W=W*1+Ẇ*θW*(tt*)θ,
(22)

where the layer width W = Wb + Ws, θ = θ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,56 except for β and αb,s. Adapting the derivation of Mikaelian,56 β governs the point at which the two power laws are matched and is computed as a dimensionless time β=Δut*σ032. 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

αb=αs=0.564712kmaxσ+4=0.022,θb=0.25,θs=0.31,
(23)
β=32,   At+=0.487,   Ui=434  m/s,   t*=3.41×103  s.
(24)

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

τ=τshock time+W̃W̃0+1Z+τ0Z+W̃/Ã1/θ,
(25)
Z=121+erfCW̃W̃0+B,
(26)

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 τ0 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

Ã=Aλ¯1θẆ0θ=0.3438,   τ0=2.075,   θ=0.296,τshocktime=0.0011Ẇ0λ¯+0.015.
(27)

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.

FIG. 10.

The predicted integral width from the new linear model [Eq. (15)] and the two power law fits [Eqs. (22) and (25)] compared to Flamenco data.

FIG. 10.

The predicted integral width from the new linear model [Eq. (15)] and the two power law fits [Eqs. (22) and (25)] compared to Flamenco data.

Close modal
TABLE I.

Mixing layer width and mix measures as a function of time for the code-averaged data for the standard problem.

t (s)τWσWΘσΘΞσΞ
0.05 0.62 0.351 0.0049 0.253 0.0334 0.247 0.0330 
0.10 1.23 0.450 0.0039 0.425 0.0404 0.403 0.0375 
0.15 1.85 0.507 0.0029 0.541 0.0435 0.512 0.0415 
0.20 2.46 0.547 0.0028 0.616 0.0426 0.589 0.0424 
0.25 3.08 0.579 0.0028 0.667 0.0395 0.644 0.0408 
0.30 3.69 0.606 0.0028 0.702 0.0357 0.683 0.0379 
0.35 4.31 0.629 0.0029 0.726 0.0317 0.712 0.0341 
0.40 4.92 0.649 0.0031 0.745 0.0280 0.735 0.0300 
0.45 5.54 0.668 0.0035 0.759 0.0249 0.753 0.0261 
0.50 6.15 0.684 0.0040 0.769 0.0225 0.767 0.0229 
t (s)τWσWΘσΘΞσΞ
0.05 0.62 0.351 0.0049 0.253 0.0334 0.247 0.0330 
0.10 1.23 0.450 0.0039 0.425 0.0404 0.403 0.0375 
0.15 1.85 0.507 0.0029 0.541 0.0435 0.512 0.0415 
0.20 2.46 0.547 0.0028 0.616 0.0426 0.589 0.0424 
0.25 3.08 0.579 0.0028 0.667 0.0395 0.644 0.0408 
0.30 3.69 0.606 0.0028 0.702 0.0357 0.683 0.0379 
0.35 4.31 0.629 0.0029 0.726 0.0317 0.712 0.0341 
0.40 4.92 0.649 0.0031 0.745 0.0280 0.735 0.0300 
0.45 5.54 0.668 0.0035 0.759 0.0249 0.753 0.0261 
0.50 6.15 0.684 0.0040 0.769 0.0225 0.767 0.0229 

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 σW0.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(tt0)θ, and the second is to utilise the derivatives of the integral width to directly determine θ1=1W/Ẇ2. This power law behavior of the mixing zone width (Wtθ) 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.57 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

=CdẆ2W,
(28)

which yields the power law solution for RMI of θ−1 = 1 + Cd.

Using non-linear regression on the code-averaged integral width for τ>1.23, assuming W = A(tt0)θ gives

W=0.807(t0.0309)0.219,
(29)

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.

TABLE II.

θ for each algorithm, along with Richardson extrapolated values of Θ and Ξ at τ = 6.15 for the standard problem. Note that Richardson extrapolation is not necessarily robust, for example, the extrapolated values for Flamenco and Ares would not be considered to be a realistic bound, hence the actual value at the latest time is also shown in parentheses.

Codeθ(τ>1.23)ΘRΞR
Ares 0.215 0.7827 (0.7442) 0.9291 (0.7351) 
Flamenco 0.220 0.8717 (0.7943) 0.6872 (0.7912) 
FLASH 0.221 0.8446 0.8027 
Hesione 0.242 0.2740a 0.2708a 
Miranda 0.228 0.8010a 0.7922a 
NUT3D 0.210 0.7677 0.7655 
Triclade 0.218 0.7685 0.7842 
TURMOIL 0.222 0.8312 0.7867a 
Codeθ(τ>1.23)ΘRΞR
Ares 0.215 0.7827 (0.7442) 0.9291 (0.7351) 
Flamenco 0.220 0.8717 (0.7943) 0.6872 (0.7912) 
FLASH 0.221 0.8446 0.8027 
Hesione 0.242 0.2740a 0.2708a 
Miranda 0.228 0.8010a 0.7922a 
NUT3D 0.210 0.7677 0.7655 
Triclade 0.218 0.7685 0.7842 
TURMOIL 0.222 0.8312 0.7867a 
a

Indicates that the measure is not uniformly converging, thus the Richardson extrapolated value is replaced by the value at the highest grid resolution.

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.,14 and Lombardini et al.12(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.59 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.

FIG. 11.

The time-dependent value of θ computed from W and its derivatives (solid line) compared to non-linear regression (dashed line) from Flamenco for the standard problem.

FIG. 11.

The time-dependent value of θ computed from W and its derivatives (solid line) compared to non-linear regression (dashed line) from Flamenco for the standard problem.

Close modal

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 α collaboration19 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.

FIG. 12.

Mix parameters Θ and Ξ for all algorithms for the standard problem (top) and the code averaged values with standard deviations (bottom).

FIG. 12.

Mix parameters Θ and Ξ for all algorithms for the standard problem (top) and the code averaged values with standard deviations (bottom).

Close modal

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 Zhou60 and taking δ = L/4 and u=2/3(TKX+TKY+TKZ) evaluated at t = 0.03 s implies that the transition time is τtr0.18 (ttr0.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 τtr0.42 (ttr0.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.61 The predicted mixing is lower than those by Orlicz et al.62 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.63 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.

FIG. 13.

Turbulent kinetic energy components and anisotropy ratio 2TKX/(TKY + TKZ) for the standard problem.

FIG. 13.

Turbulent kinetic energy components and anisotropy ratio 2TKX/(TKY + TKZ) for the standard problem.

Close modal
FIG. 14.

Code-averaged dimensional kinetic energies and anisotropy ratio for the standard problem.

FIG. 14.

Code-averaged dimensional kinetic energies and anisotropy ratio for the standard problem.

Close modal
TABLE III.

Dimensional kinetic energies and the anisotropy measure as a function of time and standard deviations for the code-averaged data for the standard problem.

t (s)τTKXσTKXTKYσTKYTKZσTKZAnisotropy
0.05 0.62 10983.94 318.25 4930.47 305.09 4936.73 316.48 2.23 
0.10 1.23 5978.18 55.41 3541.15 144.92 3664.06 132.97 1.66 
0.15 1.85 3766.78 45.97 2396.90 101.51 2502.81 90.99 1.54 
0.20 2.46 2644.78 47.58 1724.41 67.77 1797.25 76.27 1.50 
0.25 3.08 1987.92 41.08 1323.75 47.25 1369.86 54.28 1.48 
0.30 3.69 1574.77 43.72 1058.62 36.75 1092.80 41.19 1.46 
0.35 4.31 1289.17 43.06 875.09 32.92 901.20 33.67 1.45 
0.40 4.92 1084.62 41.88 736.59 28.71 764.62 26.72 1.44 
0.45 5.54 930.11 41.10 632.30 26.58 662.20 23.22 1.44 
0.50 6.15 812.86 40.63 552.11 23.37 583.82 22.72 1.43 
t (s)τTKXσTKXTKYσTKYTKZσTKZAnisotropy
0.05 0.62 10983.94 318.25 4930.47 305.09 4936.73 316.48 2.23 
0.10 1.23 5978.18 55.41 3541.15 144.92 3664.06 132.97 1.66 
0.15 1.85 3766.78 45.97 2396.90 101.51 2502.81 90.99 1.54 
0.20 2.46 2644.78 47.58 1724.41 67.77 1797.25 76.27 1.50 
0.25 3.08 1987.92 41.08 1323.75 47.25 1369.86 54.28 1.48 
0.30 3.69 1574.77 43.72 1058.62 36.75 1092.80 41.19 1.46 
0.35 4.31 1289.17 43.06 875.09 32.92 901.20 33.67 1.45 
0.40 4.92 1084.62 41.88 736.59 28.71 764.62 26.72 1.44 
0.45 5.54 930.11 41.10 632.30 26.58 662.20 23.22 1.44 
0.50 6.15 812.86 40.63 552.11 23.37 583.82 22.72 1.43 

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−2 s) in the shock direction, and τ = 0.62 (t5×102 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 τtr0.42 (ttr0.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 Wtθ, 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 k1.3 for k/k¯=3–5, which is shallower than a classical Kolmogorov spectrum or that proposed by Zhou.65 

FIG. 15.

Variable density kinetic ener gy spectra at τ = 1.23 (t = 0.1 s) (left) and τ = 6.15 (t = 0.5 s) (right) for the standard problem.

FIG. 15.

Variable density kinetic ener gy spectra at τ = 1.23 (t = 0.1 s) (left) and τ = 6.15 (t = 0.5 s) (right) for the standard problem.

Close modal

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 elsewhere11,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ȧ/λmin380–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.

TABLE IV.

θ for each algorithm in the quarter scale problem defined using non-linear regression [W = A(tt0)θ, tabulated data are dimensional], Θ, Ξ, and anisotropy measure at the final time for the quarter scale problem.

CodeExperimental error (%)Final τ [t (s)]θ(τ>24.6)At0ΘΞAnisotropy
Ares 5.8 123 (2.5) 0.248 0.2972 -0.010 42 0.804 0.8138 1.53 
Flamenco 0.5 246 (5.0) 0.296 0.2800 -0.043 27 0.807 0.8189 1.49 
FLASH 5.1 64 (1.3) 0.297 0.2787 -0.086 64 0.777 0.7882 N/A 
Miranda 3.8 246 (5.0) 0.331 0.2625 -0.103 66 0.7985 0.7963 1.55 
NUT3D 1.8 60 (0.65) 0.282 0.2819 -0.070 83 0.772 0.7794 N/A 
Triclade 0.2 123 (2.5) 0.283 0.2852 -0.040 59 0.782 0.7946 1.55 
TURMOIL 1.6 108 (2.2) 0.301 0.2825 -0.099 09 0.800 0.8095 1.66 
Mean   0.291 0.2820 -0.064 9 0.792 0.8000 1.55 
σ   0.025 0.0110 -0.034 65 0.0141 0.0144 0.06 
CodeExperimental error (%)Final τ [t (s)]θ(τ>24.6)At0ΘΞAnisotropy
Ares 5.8 123 (2.5) 0.248 0.2972 -0.010 42 0.804 0.8138 1.53 
Flamenco 0.5 246 (5.0) 0.296 0.2800 -0.043 27 0.807 0.8189 1.49 
FLASH 5.1 64 (1.3) 0.297 0.2787 -0.086 64 0.777 0.7882 N/A 
Miranda 3.8 246 (5.0) 0.331 0.2625 -0.103 66 0.7985 0.7963 1.55 
NUT3D 1.8 60 (0.65) 0.282 0.2819 -0.070 83 0.772 0.7794 N/A 
Triclade 0.2 123 (2.5) 0.283 0.2852 -0.040 59 0.782 0.7946 1.55 
TURMOIL 1.6 108 (2.2) 0.301 0.2825 -0.099 09 0.800 0.8095 1.66 
Mean   0.291 0.2820 -0.064 9 0.792 0.8000 1.55 
σ   0.025 0.0110 -0.034 65 0.0141 0.0144 0.06 

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.

FIG. 16.

Visualisations of the quarter scale problem at τ = 1.23, 25.6124, and 246 (t = 0.025, 0.52, 2.52, and 5.00 s) using Flamenco.

FIG. 16.

Visualisations of the quarter scale problem at τ = 1.23, 25.6124, and 246 (t = 0.025, 0.52, 2.52, and 5.00 s) using Flamenco.

Close modal

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.

FIG. 17.

Integral width, the time-dependent value of θ computed from W and its derivatives (solid line) compared to non-linear regression (dashed line) from Flamenco for the quarter scale problem.

FIG. 17.

Integral width, the time-dependent value of θ computed from W and its derivatives (solid line) compared to non-linear regression (dashed line) from Flamenco for the quarter scale problem.

Close modal

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 t0 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 t0 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.17 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.

FIG. 18.

Integral width from Flamenco for the standard and quarter scale problem and an additional simulation run with a sharp interface and different random numbers, plotted along with the linear multimode theory and two models [“matched” power law Eq. (22) and “blended” power law Eq. (25)] and the result of the non-linear regression.

FIG. 18.

Integral width from Flamenco for the standard and quarter scale problem and an additional simulation run with a sharp interface and different random numbers, plotted along with the linear multimode theory and two models [“matched” power law Eq. (22) and “blended” power law Eq. (25)] and the result of the non-linear regression.

Close modal

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.14 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,70 and the layer width h:16,71

M(t)=ρY1Y2dV,      h(t)=2min(f1,1f1)dx,
(30)

where . indicates a plane-averaged quantity and Yi 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:70 

Ψ(t)=ρY1Y2dxρY1Y2dx.
(31)
FIG. 19.

Integral mixing measures for the quarter scale problem and a comparison of three mix measures computed using the results from Flamenco.

FIG. 19.

Integral mixing measures for the quarter scale problem and a comparison of three mix measures computed using the results from Flamenco.

Close modal

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

Θ=0.792±0.01   Ξ=0.800±0.014,
(32)

showing excellent inter-algorithm agreement of these late time, quasi-self-similar statistics. As with the standard case, these are consistent with previous results11,12,14,17,22,63 and close to experimentally measured results of Krivets et al.61 but lower than those of Orlicz et al.62 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 374t748 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 t3θ−2 with the code-averaged θ = 0.291 that shows reasonable agreement with Flamenco, Miranda, and Triclade at late time τ>24.6.

FIG. 20.

Turbulent kinetic energy components and ratio 2TKX/(TKY + TKZ) for the quarter scale problem.

FIG. 20.

Turbulent kinetic energy components and ratio 2TKX/(TKY + TKZ) for the quarter scale problem.

Close modal

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.14 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 Zb(y, z) is computed according to

η(y,z)=y,zZb(y,z)ZbZb(y+y,z+z)Zby,zZb(y,z)Zb2.
(33)

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 λbDb(ρh+ρL)/ρh, where the post-shock densities are used.76 

FIG. 21.

(a) Visualization of the bubble front at late time τ = 6.2 (t = 0.5 s) from the standard scale problem using FLASH. Time evolution of the self-similarity ratio λb/hb from (b) the standard problem (FLASH) and (c) the quarter scale problem (FLASH and Triclade).

FIG. 21.

(a) Visualization of the bubble front at late time τ = 6.2 (t = 0.5 s) from the standard scale problem using FLASH. Time evolution of the self-similarity ratio λb/hb from (b) the standard problem (FLASH) and (c) the quarter scale problem (FLASH and Triclade).

Close modal

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 W1/λ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 4x/W4 is reasonably collapsed at all times indicating that a single length scale and, equivalently a single θ, is a reasonable approximation.

FIG. 22.

Plane-averaged volume fraction f1 plotted against scaled distance x/W for several times for the Flamenco and Triclade simulations (left) and the equivalent for f1f2 for Flamenco only (right).

FIG. 22.

Plane-averaged volume fraction f1 plotted against scaled distance x/W for several times for the Flamenco and Triclade simulations (left) and the equivalent for f1f2 for Flamenco only (right).

Close modal

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 t3θ−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 k1.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.492TKX/(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

P(k)=Cwhenkmin<k<kmax0otherwise,
(A1)

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

σ2=0P(k)dk=12πP(k)kdkydkz.
(A2)

The power in the spectrum is represented as a series of discrete modes of wave vector k = (ky, kz) with amplitudes a(k). The specified power spectrum is the sum of the squares of the discrete mode amplitudes within an annular band. As the number of modes within each band increases with k=ky2+kz2, then, the power in each individual mode must be reduced, on the mean, by a factor of k as k increases. Thus, for a given power spectrum P(k), the amplitude of each individual mode will be

a(k)=P(k)k.
(A3)

To initialise modes within a certain band, the inverse Fourier transform of this relation can be taken, noting that the amplitude is a real function,

a(y,z)=m,n=NNRecm,nexpik0my+nz,
(A4)

where kmax=2Nπ/L, the cross section is L×L, and cm,n are the amplitudes with the mode number m in the y direction and n in the z direction. To satisfy the given power spectrum, Eq. (A4) can be simplified considerably by expanding using the Euler formula and considering only the real part,

A(y,z)=m,n=kmaxkmaxrm,ncosk0my+nz+ϕm,n,
(A5)

where cm,n = c1 + ic2, thus rm,n=c12+c22 and ϕm,n=tanc2/c1. Next, the cosine term is expanded using trigonometric relations giving

A(y,z)=m,n=0Nam,ncosk0mycosk0nz+bm,ncosk0mysink0nz+cm,nsink0mycosk0nz+dm,nsink0mysink0nz.
(A6)

To initialise a random field, am,n, bm,n,…. must be chosen from a distribution randomly so that the standard deviation is proportional to the Fourier coefficients in Eq. (A3). The random variables are selected from a Gaussian distribution to give the desired mean amplitude at a given wavenumber km,n, in this case

14(ām,n2+b¯m,n2+c¯m,n2+d¯m,n2)=12πP(km,n)km,nΔkyΔkz,
(A7)

where km,n=kym2+kzn2, the wavenumber in the y direction is kym = 2πm/L, the wavenumber in the z-direction kzn = 2πn/L, and L is the width of domain, which in this case is assumed to have a square cross section.

The random seeds are generated using a Mersenne Twister routine, which is deterministic. The initial coefficients am,n, bm,n,…are rescaled such that the initial perturbation has an rms amplitude of 0.1λmin when perfectly resolved.

1.
R. D.
Richtmyer
, “
Taylor instability in shock acceleration of compressible fluids
,”
Commun. Pure Appl. Math.
13
,
297
319
(
1960
).
2.
E. E.
Meshkov
, “
Instability of the interface of two gases accelerated by a shock wave
,”
SSSR Mekh. Zhidk. Gaza
4
,
151
157
(
1969
).
3.
Y.
Zhou
, “
Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. I
,”
Phys. Rep.
(published online);
Y.
Zhou
, “
Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence and mixing. II
,”
ibid.
(published online).
4.
D. S.
Clark
 et al, “
Three-dimensional simulations of low foot and high foot implosion experiments on the National Ignition Facility
,”
Phys. Plasmas
23
(
5
),
056302
(
2016
).
5.
A. L.
Kuhl
,
J. B.
Bell
,
V. E.
Beckner
,
K.
Balakrishnan
, and
A. J.
Aspden
, “
Spherical combustion clouds in explosions
,”
Shock Waves
23
(
3
),
233
249
(
2013
).
6.
J.
Yang
,
T.
Kubota
, and
E. E.
Zukoski
, “
Applications of shock-induced mixing to supersonic combustion
,”
AIAA J.
31
(
5
),
854
862
(
1993
).
7.
A.
Burrows
, “
Supernova explosions in the universe
,”
Nature
403
,
723
733
(
2000
).
8.
D. L.
Youngs
, “
Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities
,”
Laser Part. Beams
12
,
725
750
(
1994
).
9.
D. L.
Youngs
, “
Effect of initial conditions on self-similar turbulent mixing
,” in
Proceedings of the 9th International Workshop on the Physics of Compressible Turbulent Mixing
,
2004
.
10.
D. J.
Hill
,
C.
Pantano
, and
D. I.
Pullin
, “
Large-eddy simulation and multi-scale modelling of a Richtmyer-Meshkov instability with reshock
,”
J. Fluid Mech.
557
,
29
61
(
2006
).
11.
B.
Thornber
,
D.
Youngs
,
D.
Drikakis
, and
R. J. R.
Williams
, “
The influence of initial conditions on turbulent mixing due to Richtmyer-Meshkov instability
,”
J. Fluid Mech.
654
,
99
139
(
2010
).
12.
M.
Lombardini
,
D. I.
Pullin
, and
D. I.
Meiron
, “
Transition to turbulence in shock-driven mixing: A Mach number study
,”
J. Fluid Mech.
690
,
203
226
(
2012
).
13.
V. K.
Tritschler
,
B. J.
Olson
,
S. K.
Lele
,
S.
Hickel
,
X. Y.
Hu
, and
N. A.
Adams
, “
On the Richtmyer–Meshkov instability evolving from a deterministic multimode planar interface
,”
J. Fluid Mech.
755
,
429
462
(
2014
).
14.
T.
Oggian
,
D.
Drikakis
,
D. L.
Youngs
, and
R. J. R.
Williams
, “
Computing multi-mode shock-induced compressible turbulent mixing at late times
,”
J. Fluid Mech.
779
,
411
431
(
2015
).
15.
H.
Liu
and
Z.
Xiao
, “
Scale-to-scale energy transfer in mixing flow induced by the Richtmyer-Meshkov instability
,”
Phys. Rev. E
93
,
053112
(
2016
).
16.
C. R.
Weber
,
A. W.
Cook
, and
R.
Bonazza
, “
Growth rate of a shocked mixing layer with known initial perturbations
,”
J. Fluid Mech.
725
,
372
401
(
2013
).
17.
B.
Thornber
, “
Impact of domain size and statistical errors in simulations of homogeneous decaying turbulence and the Richtmyer-Meshkov instability
,”
Phys. Fluids
28
(
4
),
045106
(
2016
).
18.
B. J.
Olson
and
J. A.
Greenough
, “
Large eddy simulation requirements for the Richtmyer-Meshkov instability
,”
Phys. Fluids
26
(
4
),
044103
(
2014
).
19.
G.
Dimonte
,
D. L.
Youngs
,
A.
Dimits
,
S.
Wunsch
,
C.
Garasi
,
M. J.
Andrews
,
A. C.
Calder
,
P.
MacNeice
,
P.
Ricker
,
S.
Weber
,
M.
Marinak
,
A.
Robinson
,
P.
Ramaprabhu
,
B.
Fryxell
,
K.
Olson
,
R.
Rosner
,
J.
Biello
,
L.
Dursi
,
F.
Timmes
,
H.
Tufo
,
Y.-N.
Young
, and
M.
Zingale
, “
A comparative study of the turbulent Rayleigh-Taylor instability using high resolution three dimensional numerical simulations: The alpha-group collaboration
,”
Phys. Fluids
16
(
5
),
1668
1693
(
2004
).
20.
R. V.
Morgan
,
R.
Aure
,
J. D.
Stockero
,
J. A.
Greenough
,
W.
Cabot
,
O. A.
Likhachev
, and
J. W.
Jacobs
, “
On the late-time growth of the two-dimensional Richtmyer-Meshkov instability in shock tube experiments
,”
J. Fluid Mech.
712
,
354
383
(
2012
).
21.
B. E.
Morgan
, “
The 2D ‘shock-jet’problem
,” Technical Report,
Lawrence Livermore National Laboratory (LLNL)
,
Livermore, CA
,
2013
.
22.
Y.
Zhou
and
B.
Thornber
, “
A comparison of three approaches to compute the effective Reynolds number of an implicit LES
,”
ASME J. Fluids Eng.
138
(
7
),
1
7
(
2016
).
23.
Y.
Zhou
,
F. F.
Grinstein
,
A. J.
Wachtor
, and
B. M.
Haines
, “
Estimating the effective Reynolds number in implicit large-eddy simulation
,”
Phys. Rev. E
89
(
1
),
013303
(
2014
).
24.
P. E.
Dimotakis
, “
The mixing transition in turbulent flows
,”
J. Fluid Mech.
409
,
69
98
(
2000
).
25.
D. L.
Youngs
, “
Three-dimensional numerical simulation of turbulent mixing by Rayleigh-Taylor instability
,”
Phys. Fluids A
3
(
5
),
1312
1320
(
1991
).
26.
A. W.
Cook
and
Y.
Zhou
, “
Energy transfer in Rayleigh-Taylor instability
,”
Phys. Rev. E
66
,
026312
(
2002
).
27.
M.
Latini
,
O.
Schilling
, and
W. S.
Don
, “
Effects of WENO flux reconstruction order and spatial resolution on reshocked two-dimensional Richtmyer-Meshkov instability
,”
J. Comput. Phys.
221
,
805
836
(
2007
).
28.
B.
Thornber
and
Y.
Zhou
, “
Energy transfer in the Richtmyer-Meshkov instability
,”
Phys. Rev. E
86
,
056302
(
2012
).
29.
K. H.
Kim
and
C.
Kim
, “
Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows part II: Multi-dimensional limiting process
,”
J. Comput. Phys.
208
,
570
615
(
2005
).
30.
E. F.
Toro
,
M.
Spruce
, and
W.
Speares
, “
Restoration of the contact surface in the HLL-Riemann solver
,”
Shock Waves
4
(
1
),
25
34
(
1994
).
31.
B.
Thornber
,
D.
Drikakis
,
R.
Williams
, and
D.
Youngs
, “
On entropy generation and dissipation of kinetic energy in high-resolution shock-capturing schemes
,”
J. Comput. Phys.
227
,
4853
4872
(
2008
).
32.
B.
Thornber
,
A.
Mosedale
,
D.
Drikakis
,
D.
Youngs
, and
R.
Williams
, “
An improved reconstruction method for compressible flows with low Mach number features
,”
J. Comput. Phys.
227
,
4873
4894
(
2008
).
33.
B.
Thornber
and
D.
Drikakis
, “
Numerical dissipation of upwind schemes in low Mach flow
,”
Int. J. Numer. Methods Fluids
56
,
1535
1541
(
2008
).
34.
R. J.
Spiteri
and
S. J.
Ruuth
, “
A class of optimal high-order strong-stability preserving time discretization methods
,”
SIAM J. Numer. Anal.
40
(
2
),
469
491
(
2002
).
35.
G.
Allaire
,
S.
Clerc
, and
S.
Kokh
, “
A five-equation model for the simulation of interfaces between compressible fluids
,”
J. Comput. Phys.
181
,
577
616
(
2002
).
36.
B.
Fryxell
,
K.
Olson
,
P.
Ricker
,
F. X.
Timmes
,
M.
Zingale
,
D. Q.
Lamb
,
P.
MacNeice
,
R.
Rosner
,
J. W.
Truran
, and
H.
Tufo
, “
FLASH: An adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes
,”
Astrophys. J., Suppl. Ser.
131
(
1
),
273
(
2000
).
37.
P.
Colella
and
P. R.
Woodward
, “
The piecewise parabolic method (PPM) for gas-dynamical simulations
,”
J. Comput. Phys.
54
(
1
),
174
201
(
1984
).
38.
B.
van Leer
, “
Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow
,”
J. Comput. Phys.
23
(
3
),
263
275
(
1977
).
39.
P.
Woodward
and
P.
Colella
, “
The numerical simulation of two-dimensional fluid flow with strong shocks
,”
J. Comput. Phys.
54
(
1
),
115
173
(
1984
).
40.
D. L.
Youngs
, “
An interface tracking method for a 3D Eulerian hydrodynamics code
,” Technical Report 44/92/35,
Atomic Weapons Research Establishment (AWRE)
,
1984
.
41.
A. W.
Cook
, “
Artificial fluid properties for large-eddy simulation of compressible turbulent mixing
,”
Phys. Fluids
19
(
5
),
055103
(
2007
).
42.
V. F.
Tishkin
,
V. V.
Nikishin
,
I. V.
Popov
, and
A. P.
Favorski
, “
Finite difference schemes of three-dimensional gas dynamics for the study of Richtmyer–Meshkov instability
,”
Matem. Mod.
7
(
5
),
15
25
(
1995
).
43.
I. G.
Lebo
and
V. F.
Tishkin
,
The Study of Hydrodynamic Instability in Problems of Laser Fusion by Methods of Mathematical Modeling
(
FIZMATLIT
,
Moscow
,
2006
).
44.
P. A.
Kuchugov
, “
Dynamics of turbulent mixing processes in laser targets
,” Ph.D. thesis,
Keldysh Institute of Applied Mathematics of RAS
,
2014
.
45.
K. V.
Vyaznikov
,
V. F.
Tishkin
, and
A. P.
Favorski
, “
Construction of monotone high resolution difference schemes for hyperbolic systems
,”
Matem. Mod.
1
(
5
),
95
120
(
1989
).
46.
R. J.
LeVeque
,
Finite Volume Methods for Hyperbolic Problems
, Cambridge Texts in Applied Mathematics (
Cambridge University Press
,
2002
), Vol. 31.
47.
V.
Daru
and
C.
Tenaud
, “
High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations
,”
J. Comput. Phys.
193
,
563
594
(
2004
).
48.
S.
Shanmuganathan
,
D. L.
Youngs
,
J.
Griffond
,
B.
Thornber
, and
R. J. R.
Williams
, “
Accuracy of high-order density-based compressible methods in low Mach vortical flows
,”
Int. J. Numer. Methods Fluids
74
(
5
),
335
358
(
2014
).
49.
Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics
, edited by
F. F.
Grinstein
,
L. G.
Margolin
, and
W. J.
Rider
(
Cambridge University Press
,
Cambridge
,
2007
).
50.
B.
van Leer
, “
Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection
,”
J. Comput. Phys.
23
,
276
299
(
1977
).
51.
W. D.
Schultz
, “
Two-dimensional Lagrangian hydrodynamic difference equations
,” in
Methods in Computational Physics
(
Academic Press
,
1964
).
52.
A. W.
Cook
,
M. S.
Ulitsky
, and
D. S.
Miller
, “
Hyperviscosity for unstructured ALE meshes
,”
Int. J. Comput. Fluid Dyn.
27
(
1
),
32
50
(
2013
).
53.
R. J. R.
Williams
, “
Sub-grid properties and artificial viscous stresses in staggered-mesh schemes
” (unpublished).
54.
J. P.
Boris
,
F. F.
Grinstein
,
E. S.
Oran
, and
R. L.
Kolbe
, “
New insights into large eddy simulation
,”
Fluid Dyn. Res.
10
,
199
228
(
1992
).
55.
D.
Drikakis
,
M.
Hahn
,
A.
Mosedale
, and
B.
Thornber
, “
Large eddy simulation using high-resolution and high-order methods
,”
Philos. Trans. R. Soc., A
367
(
1899
),
2985
2997
(
2009
).
56.
K. O.
Mikaelian
, “
Testing an analytic model for Richtmyer–Meshkov turbulent mixing widths
,”
Shock Waves
25
(
1
),
35
45
(
2015
).
57.
O.
Poujade
and
M.
Peybernes
, “
Growth rate of Rayleigh-Taylor turbulent mixing layers with the foliation approach
,”
Phys. Rev. E
81
,
016316
(
2010
).
58.
G.
Dimonte
and
M.
Schneider
, “
Turbulent Richtmyer-Meshkov instability experiments with strong radiatively driven shocks
,”
Phys. Plasmas
4
(
12
),
4347
4357
(
1997
).
59.
G. I.
Barenblatt
,
G.
Looss
, and
D. D.
Joseph
,
Nonlinear Dynamics and Turbulence
(
Pitman Publishing
,
1983
).
60.
Y.
Zhou
, “
Unification and extension of the similarity scaling criteria and mixing transition for studying astrophysics using high energy density laboratory experiments or numerical simulations
,”
Phys. Plasmas
14
(
8
),
082701
(
2007
).
61.
V. V.
Krivets
,
K. J.
Ferguson
, and
J. W.
Jacobs
, “
Turbulent mixing induced by Richtmyer-Meshkov instability
,”
AIP Conf. Proc.
1793
,
150003
(
2017
).
62.
G. C.
Orlicz
,
S.
Balasubramanian
, and
K. P.
Prestridge
, “
Incident shock Mach number effects on Richtmyer-Meshkov mixing in a heavy gas layer
,”
Phys. Fluids
25
(
11
),
114101
(
2013
).
63.
C.
Weber
,
N.
Haehn
,
J.
Oakley
,
D.
Rothamer
, and
R.
Bonazza
, “
Turbulent mixing measurements in the Richtmyer-Meshkov instability
,”
Phys. Fluids
24
(
7
),
074105
(
2012
).
64.
J. R.
Ristorcelli
,
A. A.
Gowardhan
, and
F. F.
Grinstein
, “
Two classes of Richtmyer-Meshkov instabilities: A detailed statistical look
,”
Phys. Fluids
25
(
4
),
044106
(
2013
).
65.
Y.
Zhou
, “
A scaling analysis of turbulent flows driven by Rayleigh-Taylor and Richtmyer-Meshkov instabilities
,”
Phys. Fluids
13
(
2
),
538
543
(
2001
).
66.
D. L.
Youngs
and
R. J. R.
Williams
, “
Turbulent mixing in spherical implosions
,”
Int. J. Numer. Methods Fluids
56
(
8
),
1597
1603
(
2008
).
67.
P.
Tsoutsanis
,
I. W.
Kokkinakis
,
L.
Könözsy
,
D.
Drikakis
,
R. J. R.
Williams
, and
D. L.
Youngs
, “
Comparison of structured- and unstructured-grid, compressible and incompressible methods using the vortex pairing problem
,”
Comput. Methods Appl. Mech. Eng.
293
,
207
231
(
2015
).
68.
J. T.
Morán-López
and
O.
Schilling
, “
Multicomponent Reynolds–averaged Navier–Stokes simulations of reshocked Richtmyer–Meshkov instability-induced mixing
,”
High Energy Density Phys.
9
(
1
),
112
121
(
2013
).
69.
J. T.
Morán-López
and
O.
Schilling
, “
Multi-component Reynolds-averaged Navier–Stokes simulations of Richtmyer–Meshkov instability and mixing induced by reshock at different times
,”
Shock Waves
24
(
3
),
325
343
(
2014
).
70.
Y.
Zhou
,
W. H.
Cabot
, and
B.
Thornber
, “
Asymptotic behavior of the mixed mass in Rayleigh–Taylor and Richtmyer–Meshkov instability induced flows
,”
Phys. Plasmas
23
(
5
),
052712
(
2016
).
71.
A. W.
Cook
and
P. E.
Dimotakis
, “
Transition stages of Rayleigh–Taylor instability between miscible fluids
,”
J. Fluid Mech.
443
,
69
99
(
2001
).
72.
U.
Alon
,
J.
Hecht
,
D.
Ofer
, and
D.
Shvarts
, “
Power laws and similarity of Rayleigh-Taylor and Richtmyer-Meshkov mixing fronts at all density ratios
,”
Phys. Rev. Lett.
74
(
4
),
534
537
(
1995
).
73.
Y.
Srebro
,
Y.
Elbaz
,
O.
Sadot
,
L.
Arazi
, and
D.
Shvarts
, “
A general buoyancy-drag model for the evolution of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities
,”
Laser Part. Beams
21
(
3
),
347
353
(
2003
).
74.
G.
Dimonte
and
M.
Schneider
, “
Density ratio dependence of Rayleigh-Taylor mixing for sustained and impulsive acceleration histories
,”
Phys. Fluids
12
,
304
321
(
2000
).
75.
P.
Ramaprabhu
,
G.
Dimonte
, and
M. J.
Andrews
, “
A numerical study of the influence of initial perturbations on the turbulent Rayleigh-Taylor instability
,”
J. Fluid Mech.
536
,
285
319
(
2005
).
76.
B. J.
Daly
, “
Numerical study of two fluid Rayleigh-Taylor instability
,”
Phys. Fluids
10
(
2
),
297
307
(
1967
).