# Numerical Simulation of an Idealised Richtmyer-Meshkov Instability Shock
Tube Experiment

Michael Groom1<EMAIL_ADDRESS>Ben Thornber1  1School of
Aerospace, Mechanical and Mechatronic Engineering, University of Sydney,
Sydney, NSW 2006, Australia

###### Abstract

The effects of initial conditions on the evolution of the Richtmyer–Meshkov
instability (RMI) at early to intermediate time are analysed, using numerical
simulations of an idealised version of recent shock tube experiments performed
at the University of Arizona (Sewell et al., _J. Fluid Mech._ (2021), 917,
A41). The experimental results are bracketed by performing both implicit
large-eddy simulations (ILES) of the high-Reynolds number limit as well as
direct numerical simulations (DNS) at Reynolds numbers lower than those
observed in the experiments, both using the Flamenco finite-volume code.
Various measures of the mixing layer width $h$, known to scale as $\sim
t^{\theta}$ at late time, based on both the plane-averaged turbulent kinetic
energy (TKE) and volume fraction (VF) profiles are used to explore the effects
of initial conditions on $\theta$ and are compared with the experimental
results. The decay rate $n$ of the total fluctuating kinetic energy is also
used to estimate $\theta$ based on a relationship that assumes self-similar
growth of the mixing layer. The estimates for $\theta$ range between 0.44 and
0.52 for each of the broadband perturbations considered and are in good
agreement with the experimental results. Decomposing the mixing layer width
into separate bubble and spike heights $h_{b}$ and $h_{s}$ shows that, while
the bubbles and spikes initially grow at different rates, their growth rates
$\theta_{b}$ and $\theta_{s}$ have equalised by the end of the simulations
indicating that the mixing layer is approaching self-similarity. Anisotropy of
the Reynolds stresses is also analysed and is shown to persist throughout each
of the simulations. Outer-scale Reynolds numbers and various key length scales
are calculated for the DNS cases, showing that fully developed turbulence is
not obtained due to the challenges associated with performing DNS for
broadband initial conditions. Overall the results demonstrate important
differences between broadband and narrowband surface perturbations, as well as
persistent effects of finite bandwidth on the growth rate of mixing layers
evolving from broadband perturbations. Good agreement is obtained with the
experiments for the different quantities considered, however the results also
show that care must be taken when using measurements based on the velocity
field to infer properties of the concentration field, as well as when it is
appropriate to assume the mixing layer is growing self-similarly with a single
growth rate $\theta$.

###### keywords:

shock waves, turbulent mixing, transition to turbulence

##  1 Introduction

This paper analyses the effects of initial conditions on the evolution of the
Richtmyer–Meshkov instability (RMI), which occurs when an interface separating
two materials of differing densities is accelerated impulsively, typically by
an incident shock wave (Richtmyer, 1960; Meshkov, 1969). The instability
evolves due to the deposition of baroclinic vorticity at the interface, caused
by a misalignment of density and pressure gradients during the shock–interface
interaction. This occurs either from surface perturbations on the interface,
or when the shock wave is non-uniform or inclined relative to the interface.
The baroclinic vorticity that is deposited on the interface leads to the
growth of surface perturbations and the development of secondary shear layer
instabilities, which drive the transition to a turbulent mixing layer. Unlike
the closely related Rayleigh–Taylor instability (RTI), the RMI is induced for
both light to heavy and heavy to light configurations. In both cases the
initial growth of the interface is linear in time and can be described by
analytical expressions (Richtmyer, 1960; Meyer & Blewett, 1972;
Vandenboomgaerde et al., 1998). However, as the amplitudes of modes in the
perturbation become large with respect to their wavelengths the growth becomes
nonlinear, whereby numerical simulation is required to calculate the
subsequent evolution of the mixing layer. Another key difference between RTI
and RMI is that, for the RMI, baroclinic vorticity is only deposited initially
and not continuously generated, compared to the (classical) RTI where the
interface is continuously accelerated. For a comprehensive and up-to-date
review of the literature on both RTI, RMI and the Kelvin-Helmholtz instability
(KHI), the reader is referred to Zhou (2017a, b); Zhou et al. (2021), as well
as Livescu (2020) for an excellent review on variable-density turbulence more
generally.

The understanding of mixing due to RMI is of great importance in areas such as
inertial confinement fusion (ICF) (Lindl et al., 2014), where a spherical
capsule containing thermonuclear fuel is imploded using powerful lasers with
the aim of compressing the contents to sufficient pressures and temperatures
so as to initiate nuclear fusion. The compression is performed using a series
of strong shocks, which trigger hydrodynamic instabilities at the ablation
front due to capsule defects and drive asymmetries (Clark et al., 2016). The
subsequent mixing of ablator material and fuel that ensues can dilute and cool
the hotspot, which reduces the overall efficiency of the implosion. As a
contrast to ICF, in high-speed combustion such as in a scramjet or rotating
detonation engine, RMI due to weak shocks improves the mixing of fuel and
oxidiser leading to more efficient combustion (Yang et al., 1993, 2014). An
understanding of mixing due to RMI is also important for many astrophysical
phenomena such as supernovae and the dynamics of interstellar media (Arnett,
2000). Note that in such applications RTI usually occurs alongside RMI and in
general it is impossible to separate the effects of both instabilities.
However, there is still great value in studying RMI independently,
particularly when comparing with shock tube experiments that have been
designed to isolate its effects using an RT-stable configuration.

In the applications mentioned above, the most important statistical quantity
one would like to know is typically the mixing layer width, denoted by $h$. At
late time $h$ scales as $\sim t^{2}$ for RTI and $\sim t^{\theta}$ for RMI
where the exponent $\theta\leq 1$ has been shown to depend on initial
conditions (Youngs, 2004; Thornber et al., 2010). Various approaches have been
taken to define $h$, which fall into one of two categories. The first is to
consider the distance between two cutoff locations based on a particular
threshold of some spatially-averaged profile in the direction normal to the
mixing layer (i.e. the direction of the shock-induced acceleration). Examples
include the visual width (Cook & Dimotakis, 2001) based on the 1% and 99%
locations of the mean volume fraction profile (the choice of a 1% threshold is
somewhat arbitrary; see Zhou & Cabot (2019) for a comparison of different
thresholds in the context of RTI). Such measures have the advantage of being
easily interpretable but can be sensitive to statistical fluctuations. The
second approach is to define an integral measure by integrating a particular
spatially-averaged profile in the normal direction, for example the integral
width (Andrews & Spalding, 1990). Integral measures are less susceptible to
statistical fluctuations but are also less interpretable, as different
profiles can give the same integrated value. The recently proposed mixed mass
(Zhou et al., 2016) and integral bubble and spike heights (Youngs & Thornber,
2020a) are attempts to combine the best aspects of both approaches.

Over the last few decades, both shock tube experiments and numerical
simulations have been performed in order to better understand the fundamentals
of RMI, such as the value of $\theta$ at late time. Previous numerical studies
have typically used large-eddy simulation (LES) or implicit LES (ILES) to
predict mixing at late time in the high Reynolds number limit (Youngs, 1994;
Hill et al., 2006; Thornber et al., 2010; Lombardini et al., 2012; Tritschler
et al., 2014a; Thornber et al., 2017; Soulard et al., 2018). Key findings
include the dependence of $\theta$ on the type of surface perturbation used to
initiate the instability (Youngs, 2004; Thornber et al., 2010). Narrowband
perturbations, which include only a small, annular band of modes in wavenumber
space, have been found to give values of $\theta$ at late-time between 0.25
(Soulard & Griffond, 2022) and 0.33 (Youngs & Thornber, 2020b) whereas
perturbations including additional long wavelength modes, known as broadband
perturbations, have been found to give values of $\theta$ as high as 0.75
(Groom & Thornber, 2020). Studies of the effects of initial conditions in RTI
have found similar results for the growth rate $\alpha$ when additional long
wavelength modes were included in the initial perturbation (Ramaprabhu et al.,
2005; Banerjee & Andrews, 2009). When only short wavelength perturbations are
present the growth rate of RTI is limited by the nonlinear coupling of
saturated short wavelength modes (bubble merger), while additional long
wavelength perturbations cause the growth rate to become limited by the
amplification and saturation of long wavelength modes (bubble competition).
Futhermore, Aslangil et al. (2020) considered the case of RTI where the
applied acceleration is completely withdrawn after initial development. The
resulting mixing layer is closely related to an RMI-induced mixing layer,
differing only by the mechanism of the initial acceleration, with the growth
rate exponent for narrowband initial conditions shown to be within the bounds
of 0.2 to 0.28 suggested by Weber et al. (2013).

Early shock tube experiments made use of membranes to form the initial
perturbation between the two gases (Vetter & Sturtevant, 1995), however these
tended to leave fragments that dampened the subsequent instability growth,
inhibited mixing and interfered with diagnostics. In order to circumvent this,
modern shock tube experiments use membraneless interfaces, for example by
forming by a shear layer between counter-flowing gases (Weber et al., 2012,
2014; Reese et al., 2018; Mohaghar et al., 2017, 2019), using a gas curtain
(Balakumar et al., 2008; Balasubramanian et al., 2012) or by using
loudspeakers to generate Faraday waves at the interface (Jacobs et al., 2013;
Krivets et al., 2017; Sewell et al., 2021).

These methods of interface generation typically result in the formation of a
broadband surface perturbation and as such these experiments have obtained
values of $\theta$ that are higher than the 0.25–0.33 expected for narrowband
initial conditions. For example Weber et al. (2012, 2014) measured $\theta$ in
the range 0.43–0.58, while later experiments on the same facility by Reese et
al. (2018) obtained $\theta=0.34\pm 0.01$ once the concentration field was
adjusted to remove larger-scale structures from the mixing layer prior to
averaging in the spanwise direction. Jacobs et al. (2013) found that their
measurements of mixing layer width prior to reshock could be partitioned into
two groups with different power law exponents. The particular diagnostic used
was the mixing layer half width, found by taking the distance between the 10%
and 90% average concentration locations and halving this. Prior to reshock,
both groups initially had growth rates close to 0.5 ($\theta=0.51$ and
$\theta=0.54$), while at later times the growth rates were smaller but also
more different ($\theta=0.38$ and $\theta=0.29$ respectively). Krivets et al.
(2017) also found a wide range of $\theta$ for the integral width prior to
reshock, ranging from $\theta=0.18$ to $\theta=0.57$, using a similar
experimental setup. During these experiments the timing of the arrival of the
shock wave relative to the phase of the forcing cycle was not controlled,
which resulted in large variations in the initial amplitudes of the
perturbation. More recent experiments by Sewell et al. (2021) took this into
account and divided the results into a low-amplitude and high-amplitude group.
Using a measure for the mixing layer width based on 5% threshold locations of
the turbulent kinetic energy profile, they found $\theta=0.45\pm 0.08$ and
$\theta=0.51\pm 0.04$ for the low- and high-amplitude groups prior to reshock.

In this paper, both ILES and direct numerical simulations (DNS) are performed
of 3D RMI with narrowband and broadband perturbations, using a setup that
represents an idealised version of the shock tube experiments performed at the
University of Arizona (Jacobs et al., 2013; Krivets et al., 2017; Sewell et
al., 2021) to investigate the effects of long wavelength modes in the initial
perturbation. A similar study was performed in Groom & Thornber (2020) but the
main aim in that paper was to approximate the regime where there are always
longer and longer wavelength modes in the initial condition that are yet to
saturate (referred to as the infinite bandwidth limit). Of primary interest
here is to explore the impacts of finite bandwidth broadband perturbations on
the mixing layer growth over the length and time scales of a typical shock
tube experiment and compare the results with those of both narrowband
perturbations and broadband perturbations in the infinite bandwidth limit.
While the main aim is not to match the experiments as closely as possible, it
is anticipated that the results generated in this study could in principle be
verified experimentally. Direct comparisons are also still able to be made
through appropriate non-dimensionalisations, which has previously been
difficult to do when comparing results between simulations and experiments. An
assessment will also be made as to the validity of using measurements based on
the velocity field to draw conclusions about the concentration field (and vice
versa).

The paper is organised as follows. In §2, an overview of the governing
equations and numerical methods employed to solve these equations is given, as
well as a description of the computational setup and initial conditions. This
section also gives a brief discussion on some of the challenges associated
with performing DNS with broadband surface perturbations. §3 details an
analysis of many of the same quantities presented in Sewell et al. (2021),
including turbulent kinetic energy profiles and spectra as well as various
measures of the mixing layer width that are used to estimate the growth rate
$\theta$. The evolution of key length scales and Reynolds numbers is also
given for the DNS cases. Finally, §4 gives a summary of the main findings, as
well as directions for future work on this problem.

Figure 1: A schematic of the problem setup. The major ticks correspond to a
grid spacing of $\Delta x=1.0$ m. The interface is initially located at
$x=3.0$ m and the shock is initially located at $x=2.5$ m in the light fluid
and travels from light to heavy.

##  2 Computational Setup

###  2.1 Governing Equations

The computations presented in this paper all solve the compressible
Navier–Stokes equations extended to a five-equation, quasi-conservative system
of equations based on volume fractions rather than the conventional four-
equation, fully-conservative model based on mass fractions for multicomponent
flows. This ensures that pressure and temperature equilibrium is maintained
across material interfaces when upwind discretisations are used and the ratio
of specific heats varies across the interface, as is the case for air and SF6,
which greatly improves the accuracy and efficiency of the computation (Allaire
et al., 2002; Massoni et al., 2002). This is a well-established approach for
inviscid computations and was recently extended to include the effects of
species diffusion, viscosity and thermal conductivity by Thornber et al.
(2018), enabling accurate and efficient DNS to be performed for this class of
problems. The full set of equations for binary mixtures is

∂ρ∂t

+∇(ρu) & = 0
∂ρu∂t+∇(ρuu^t+pδ) = ∇σ
∂ρe∂t+∇([ρe+p]u) = ∇(σu-q)
∂ρ1f1∂t+∇(ρ_1 f_1u) = ∇(ρD_12W1f1W)
∂f1∂t+u f_1 = (D_12 f_1)-MD_12 f_1 f_1+D_12 f_1NN.

In (2.1), $\rho$ is the mass density, $\boldsymbol{u}=[u,v,w]^{t}$ is the
mass-weighted velocity vector, $p$ is the pressure, $f_{n}$ is the volume
fraction of species $n$ and $e=e_{i}+e_{k}$ is the total energy per unit mass,
where $e_{k}=\frac{1}{2}\boldsymbol{u\cdot u}$ is the kinetic energy and the
internal energy $e_{i}$ is given by the equation of state. Note that only
(2.1$e$) is in non-conservative form, hence the term quasi-conservative as
conservation errors are negligible (only species internal energies are not
conserved). All computations are performed using the ideal gas equation of
state

$e_{i}=\frac{p}{\rho(\overline{\gamma}-1)}$ (1)

where $\overline{\gamma}$ is the ratio of specific heats of the mixture. For
the five-equation model this is given by

$\frac{1}{\overline{\gamma}-1}=\sum_{n}\frac{f_{n}}{\gamma_{n}-1}$ (2)

which is an isobaric closure (individual species temperatures are retained in
the mixture). The viscous stress tensor $\boldsymbol{\sigma}$ for a Newtonian
fluid is

$\boldsymbol{\sigma}=-\overline{\mu}\big{[}\boldsymbol{\nabla
u}+(\boldsymbol{\nabla
u})^{t}\big{]}+\frac{2}{3}\overline{\mu}(\boldsymbol{\nabla\cdot
u})\boldsymbol{\delta}$ (3)

where $\overline{\mu}$ is the dynamic viscosity of the mixture. Note that in
(3) the bulk viscosity is assumed to be zero according to Stokes’ hypothesis.
The heat flux $\boldsymbol{q}=\boldsymbol{q}_{c}+\boldsymbol{q}_{d}$, with the
conductive heat flux $\boldsymbol{q_{c}}$ given by Fourier’s law

$\boldsymbol{q}_{c}=-\overline{\kappa}\boldsymbol{\nabla}T$ (4)

where $\overline{\kappa}$ is the thermal conductivity of the mixture, and $T$
is the temperature. The thermal conductivity of species $n$ is calculated
using kinetic theory as
$\kappa_{n}=\mu_{n}\left(\frac{5}{4}\frac{\mathcal{R}}{W_{n}}+c_{p,n}\right)$,
while the thermal conductivity of the mixture (as well as the mixture
viscosity) is calculated using Wilke’s rule. The enthalpy flux
$\boldsymbol{q}_{d}$, arising from changes in internal energy due to mass
diffusion, is given by

$\boldsymbol{q}_{d}=\sum_{n}h_{n}\boldsymbol{J}_{n}$ (5)

where $h_{n}=c_{p,n}T$ is the enthalpy of species $n$ and $c_{p,n}$ the
specific heat at constant pressure. The diffusion flux on the RHS of (2.1$d$)
invokes Fick’s law of binary diffusion, written in terms of volume fraction.
$W_{n}$ is the molecular weight of species $n$, $W$ is the molecular weight of
the mixture and the binary diffusion coefficient $D_{12}$ is calculated by
assuming both species have the same Lewis number
($\mbox{{Le}}_{1}=\mbox{{Le}}_{2}=\mbox{{Le}}$), such that

$D_{12}=\frac{\overline{\kappa}}{\mbox{{Le}}\rho\bar{c}_{p}}$ (6)

with $\bar{c}_{p}$ the specific heat at constant pressure for the mixture.
Finally in (2.1$e$), $\mathcal{M}=\frac{W_{1}-W_{2}}{W_{1}f_{1}+W_{2}f_{2}}$
and $N=p/k_{b}T$ is the number density.

###  2.2 Numerical method

The governing equations presented in §2.1 are solved using the University of
Sydney code Flamenco, which employs a method of lines discretisation approach
in a structured, multiblock framework. Spatial discretisation is performed
using a Godunov-type finite-volume method, which is integrated in time via a
second-order TVD Runge-Kutta method (Spiteri & Ruuth, 2002). The spatial
reconstruction of the inviscid terms uses a fifth-order MUSCL scheme (Kim &
Kim, 2005), which is augmented by a modification to the reconstruction
procedure to ensure the correct scaling of pressure, density and velocity
fluctuations in the low Mach number limit (Thornber et al., 2008). The
inviscid flux component is calculated using the HLLC Riemann solver (Toro et
al., 1994), while the viscous and diffusive fluxes are calculated using
second-order central differences. Following Abgrall (1996), the non-
conservative volume fraction equation is written as a conservative equation
minus a correction term

$\frac{\partial f_{1}}{\partial
t}+\boldsymbol{\nabla}\bcdot(\mathcal{U}f_{1})-f_{1}(\boldsymbol{\nabla}\bcdot\mathcal{U})=\bnabla\bcdot(D_{12}\bnabla{f_{1}})$
(7)

with $\mathcal{U}=\boldsymbol{u}+\mathcal{M}D_{12}\bnabla
f_{1}-D_{12}\frac{\bnabla N}{N}$. The additional terms in $\mathcal{U}$ that
arise from species diffusion must be included in the calculation of the
inviscid flux component, as even though they are viscous in nature they modify
the upwind direction of the advection of volume fraction in the solution to
the Riemann problem at each cell interface. In the HLLC Riemann solver used in
Flamenco this is achieved by modifying the wave speeds to incorporate the
additional diffusion velocity, see Thornber et al. (2018) for further details.
In the absence of viscosity and thermal conductivity the governing equations
reduce to the inviscid five-equation model of Allaire et al. (2002), which has
been used in previous studies of RMI (Thornber, 2016; Thornber et al., 2017).
The numerical algorithm described above has been extensively demonstrated to
be an effective approach for both ILES and DNS of shock-induced turbulent
mixing problems (see Thornber et al., 2010, 2011; Groom & Thornber, 2019,
2021).

###  2.3 Problem Description and Initial Conditions

The computational setup is similar to previous studies of narrowband and
broadband RMI by Groom & Thornber (2019, 2020) but with a few key differences
that will be described here. A Cartesian domain of dimensions $x\times y\times
z=L_{x}\times L\times L$ where $L=2\upi$ m is used for all simulations. The
extent of the domain in the $x$-direction is either $L_{x}=1.5\upi$ for the
ILES cases or $L_{x}=0.75\upi$ for the DNS cases. Periodic boundary conditions
are used in the $y$\- and $z$-directions, while in the $x$-direction outflow
boundary conditions are imposed very far away from the test section so as to
minimise spurious reflections from outgoing waves impacting the flow field.
The initial mean positions of the shock wave and the interface are $x_{s}=2.5$
m and $x_{0}=3.0$ m respectively and the initial pressure and temperature of
both (unshocked) fluids is $p=0.915$ atm and $T=298$ K, equal to that in the
experiments of Jacobs et al. (2013). All computations employ the ideal gas
equation of state with a fixed value of $\gamma$ for each species. A schematic
of the initial condition is shown in Figure 1.

The shock Mach number is $M=1.5$, which is higher than the $M=1.2$ shock used
in Jacobs et al. (2013); Krivets et al. (2017) and the $M=1.17$ shock used in
Sewell et al. (2021). This is so that the initial velocity jump is larger,
which makes more efficient use of the explicit time stepping algorithm, but
not so large that it introduces significant post-shock compressibilty effects.
Therefore the post-shock evolution of the mixing layer is still approximately
incompressible in both the present simulations and the experiments in (Jacobs
et al., 2013; Krivets et al., 2017; Sewell et al., 2021). The initial
densities of air and SF6 are $\rho_{1}=1.083$ kg/m3 and $\rho_{2}=5.465$ kg/m3
and the post-shock densities are $\rho_{1}^{+}=2.469$ kg/m3 and
$\rho_{2}^{+}=15.66$ kg/m3 respectively. This gives a post-shock Atwood number
of $A^{+}=0.72$, which is essentially the same as the value of $0.71$ given in
Jacobs et al. (2013), indicating that the effects of compressibilty are
minimal. The variation in $\rho$ and $f_{1}$ across the interface are computed
based on the surface perturbation described in (7) below. The evolution of the
interface is solved in the post-shock frame of reference by applying a shift
of $\Delta u=-158.08$ m/s to the initial velocities of the shocked and
unshocked fluids. The initial velocity field is also modified to include an
initial diffusion velocity at the interface, which is calculated as in
previous DNS studies of RMI (Groom & Thornber, 2019, 2021). To improve the
quality of the initial condition, three-point Gaussian quadrature is used in
each direction to accurately compute the cell averages required by the finite-
volume algorithm.

Property | Air | SF6
---|---|---
$W_{l}$ | 28.964 | 146.057
$\gamma_{l}$ | 1.4 | 1.1
$\mu_{l}$ | 1.836 | 1.535
$\Pran_{l}$ | 0.71 | 0.90
$\mbox{{Sc}}_{l}$ | 0.71 | 0.90
Table 1: The molecular weight $W_{l}$ (g/mol), ratio of specific heats
$\gamma$, dynamic viscosities ($\times 10^{5}$ Pa-s) and Prandtl and Schmidt
numbers of air and SF6.

Table 1 gives the thermodynamic properties of each fluid. The dynamic
viscosities of both fluids are calculated using the Chapman–Enskog viscosity
model at a temperature of $T=298$ K, while the diffusivities are calculated
under the assumption of Lewis number equal to unity (hence
$\Pran_{l}=\mbox{{Sc}}_{l}$). In the DNS calculations, the actual values of
viscosity used are much higher, so as to give a Reynolds number that is able
to be fully resolved, but are kept in the same proportion to each other. This
is so that the same domain width $L$ can be used for each calculation.

Based on the interface characterisation of the low-amplitude set of
experiments performed in Sewell et al. (2021), four different initial surface
perturbations of a planar interface are considered which follow an idealised
power spectrum of the form

$P(k)=Ck^{m}.$ (8)

Three broadband initial conditions are simulated, containing length scales in
the range $\lambda_{max}=L/2$ to $\lambda_{min}=L/32$ and with a spectral
exponent $m=-1$, $-2$ and $-3$ respectively. The choice of bandwidth
$R=\lambda_{max}/\lambda_{min}=16$ is based on estimates of the minimum
initial wavelength performed in Jacobs et al. (2013) of $\lambda_{min}=2.9$ to
$3.2$ mm, relative to a test section width of $L=8.9\times 10^{-2}$ m. When
scaled to the dimensions of the experiment, the perturbations in this study
all have a minimum wavelength of $\lambda_{min}=2.8$ mm. Note also that the
diagnostic spatial resolution of the PIV method used in Sewell et al. (2021)
is $1.98$ mm, resulting in attenuation of the measured scales that are smaller
than this. The constant $C$ dictates the overall standard deviation of the
perturbations and is set such that all initial amplitudes are linear and each
perturbation has the same amplitude in the band between $k_{max}/2$ and
$k_{max}$, specifically $a_{k_{max}}k_{max}=1$. See Groom & Thornber (2020)
for further details, noting that unlike the broadband perturbations analysed
in that study the perturbations considered here have different total standard
deviations for the same bandwidth.

The power spectra for these three perturbations are shown in Figure 2, along
with the mean power spectrum of the low-amplitude experiments from Sewell et
al. (2021).

Figure 2: Power spectra of the broadband perturbations as well as the mean
power spectrum of the low-amplitude experiments from Sewell et al. (2021).
Note that the spectra are scaled to match the dimensions of the experiment.

In Figure 2 it can be seen that the $m=-3$ initial condition is the closest
match to the experiments (with an estimated slope of $m=-2.99$ over the same
range of modes), with the other perturbations included to study the effects of
varying $m$. A fourth perturbation (not shown) is also considered; a
narrowband perturbation with a constant power spectrum (i.e. $m=0$) and length
scales in the range $\lambda_{min}=L/16$ to $\lambda_{max}=L/32$. This is used
to study the effects of additional long wavelength modes in the initial
condition and is essentially the same perturbation as the quarter-scale scale
case in Thornber et al. (2017), however the initial amplitudes are larger and
are defined such that $a_{k_{max}}k_{max}=1$, which is at the limit of the
linear regime. Note that in the experiments of Jacobs et al. (2013),
$a_{k_{max}}k_{max}$ ranged between 2.82 and 3.14, which is much more
nonlinear. The choice of restricting the mode amplitudes such that all modes
are initially linear is made so that the results may be easily scaled by the
initial growth rate and compared with the results of the previous studies.

The amplitudes and phases of each mode are defined using a set of random
numbers that are constant across all grid resolutions and cases, thus allowing
for a grid convergence study to be performed for each case. The interface is
also initially diffuse for this same reason, with the profile given by an
error function with characteristic initial thickness
$\updelta=\lambda_{min}/4$. The volume fractions $f_{1}$ and $f_{2}=1-f_{1}$
are computed as

$f_{1}(x,y,z)=\frac{1}{2}\textrm{erfc}\left\\{\frac{\sqrt{\upi}\left[x-S(y,z)\right]}{\updelta}\right\\}$
(9)

where $S(y,z)=x_{0}+A(y,z)$, with $A(y,z)$ being the amplitude perturbation
satisfying the specified power spectrum and $x_{0}$ the mean position of the
interface. The amplitude perturbation $A(y,z)$ is given by

$\displaystyle
A(y,z)=\sum_{m,n=0}^{{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}N_{max}}}}\big{[}$
$\displaystyle a_{mn}$
$\displaystyle\cos(mk_{0}y)\cos(nk_{0}z)+b_{mn}\cos(mk_{0}y)\sin(nk_{0}z)$
$\displaystyle+$ $\displaystyle c_{mn}$
$\displaystyle\sin(mk_{0}y)\cos(nk_{0}z)+d_{mn}\sin(mk_{0}y)\sin(nk_{0}z)\big{]}$
(10)

where
${{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}N_{max}}}=k_{max}L/(2\upi)$,
$k_{0}=2\upi/L$ and $a_{mn}\ldots d_{mn}$ are selected from a Gaussian
distribution. Crucially, the Mersenne Twister pseudorandom number generator is
employed which allows for the same random numbers to be used across all
perturbations. This facilitates grid convergence studies for DNS and ensures
that the phases of each mode are identical when comparing across perturbations
with different values of $m$; only the amplitudes are varied. For full details
on the derivation of the surface perturbation see Thornber et al. (2010, 2017)
and Groom & Thornber (2020). A visualisation of each initial perturbation is
shown in figure 3. Whilst there is a noticeable difference between the
narrowband and broadband surface perturbations, the differences between the
$m=-1$ and $m=-2$ perturbations in particular are quite subtle. Nevertheless
these subtle differences in the amplitudes of the additional, longer
wavelengths are responsible for quite noticeable differences in the subsequent
evolution of the mixing layer, as will be shown in the following sections.
This highlights the importance of understanding the sensitivity to initial
conditions in RMI-induced flows.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 3: Contours of volume fraction $f_{1}$ for the ILES cases at $t=0$ and
$z=0$. The major ticks on both axes correspond to a grid spacing of $\Delta
x={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
y=}}1$ m.

For each perturbation, the weighted-average wavelength can be defined as
$\bar{\lambda}=2\upi/\bar{k}$, where

$\bar{k}=\frac{\sqrt{\displaystyle\int_{k_{min}}^{k_{max}}k^{2}P(k)\>\mathrm{d}k}}{\sqrt{\displaystyle\int_{k_{min}}^{k_{max}}P(k)\>\mathrm{d}k}}.$
(11)

Similarly, the initial growth rate of the perturbation variance is given by

$\dot{\sigma_{0}}=\sigma_{0}^{+}A^{+}\Delta u\bar{k}/\psi$ (12)

where $\sigma_{0}^{+}=C_{V}(1-\Delta u/U_{s})\sigma_{0}$ is the post-shock
standard deviation, $\sigma_{0}$ is the initial standard deviation and $\psi$
is a correction factor to account for the diffuse interface (Duff et al.,
1962; Youngs & Thornber, 2020b). Here $C_{V}=(A^{-}+C_{R}A^{+})/(2C_{R}A^{+})$
is an additional correction factor that is applied to the Richtmyer
compression factor $C_{R}=(1-\Delta u/U_{s})$ to give the impulsive model of
Vandenboomgaerde et al. (1998). For the present gas combination and
configuration, $C_{V}=1.16$ and is used to account for deficiencies in the
original impulsive model of Richtmyer (1960) for certain cases. Thornber et
al. (2017) showed that for a Gaussian height distribution, the integral width
$W=\int\langle f_{1}\rangle\langle f_{2}\rangle\>\mathrm{d}x$ is equal to
$0.564\sigma$ and therefore $\dot{W_{0}}=0.564\dot{\sigma_{0}}$. For the DNS
cases, the initial Reynolds number is calculated in line with previous studies
as

$\Rey_{0}=\frac{\bar{\lambda}\dot{W}_{0}\overline{\rho^{+}}}{\overline{\mu}}$
(13)

$\overline{\rho^{+}}=9.065$ kg/m3 is the mean post-shock density. Table 2
gives the initial growth rate and weighted-average wavelength for each
perturbation.

Quantity | $m=0$ | $m=-1$ | $m=-2$ | $m=-3$
---|---|---|---|---
$R$ | 2 | 16 | 16 | 16
$\bar{\lambda}$ | 0.278 | 0.463 | 0.785 | 1.33
$\dot{W_{0}}$ | 16.74 | 20.03 | 23.84 | 34.32
Table 2: The bandwidth, weighted-average wavelength (m) and initial growth
rate of integral width (m/s) for each of the four perturbations.

###  2.4 Direct Numerical Simulations

Prior to presenting results for each perturbation, it is important to discuss
some of the challenges present when performing DNS of RMI with broadband
perturbations. Previous DNS studies of 3D multi-mode RMI have focussed
exclusively on narrowband perturbations (Olson & Greenough, 2014; Groom &
Thornber, 2019; Wong et al., 2019; Groom & Thornber, 2021) or perturbations
with a dominant single mode (Tritschler et al., 2014b). The present set of
broadband DNS use a perturbation with $8\times$ the bandwidth of initial modes
compared to the narrowband perturbation analysed in Groom & Thornber (2019,
2021), but still require the same number of cells per initial minimum
wavelength for a given Reynolds number in order to fully resolve the
calculation. To be considered fully resolved and thus qualify as ”strict” DNS,
grid convergence must be demonstrated for statistics that depend on the
smallest scales in the flow, such as enstrophy and scalar dissipation rate. Of
the previously cited studies, only Groom & Thornber (2019, 2021) fully resolve
these gradient-dependent quantities and none of the studies mentioned (as well
as the present study) resolve the internal structure of the shock wave.
Demonstration of grid convergence for enstrophy and scalar dissipation rate in
the present set of DNS cases is given in Appendix A, however this comes at the
cost of limiting the Reynolds number that can be achieved, as discussed below.

Regarding the Reynolds number, using the standard width-based definition
$\Rey_{h}=\dot{h}h/\nu$ where the width $h\propto t^{\theta}$ then the
Reynolds number, and hence the grid resolution requirements, can either
increase or decrease in time depending on the value of $\theta$ since

$\Rey_{h}\propto\frac{\theta t^{\theta-1}t^{\theta}}{\nu}\propto
t^{2\theta-1}.$ (14)

Therefore for $\theta<1/2$ the Reynolds number is decreasing and vice versa
for $\theta>1/2$. Youngs (2004); Thornber et al. (2010) showed that the value
of $\theta$ depends on both the bandwidth and spectral slope $m$ of the
initial condition, which was recently demonstrated in Groom & Thornber (2020)
using ILES for perturbations of the form given by (8) with $m=-1$, $-2$ and
$-3$. For the largest bandwidths simulated, these perturbations gave values of
$\theta=0.5$, $0.63$ and $0.75$ respectively, which for the $m=-1$ and $-2$
cases are quite close to the theoretical values of $\theta=1/2$ and
$\theta=2/3$. What these results imply is that the Reynolds number of a
broadband perturbation with $m\leq-1$ will either be constant or increase with
time as the layer develops, which make performing fully grid-resolved DNS more
challenging than for a narrowband layer where $\theta\leq 1/3$ (Elbaz &
Shvarts, 2018; Soulard et al., 2018).

For DNS of narrowband RMI the number of cells per $\lambda_{min}$ can be
maximised, which sets the smallest scale that can be grid resolved and
therefore the maximum Reynolds number that can be obtained on a given grid.
For fully developed isotropic turbulence, it is well known that grid
resolution requirements scale as $\Rey^{9/4}$ and the total number of floating
point operations required to perform a simulation to a given time scales as
$Re^{3}$ (Pope, 2000). For transitional RMI, empirically the scaling appears
to be less severe (closer to $\Rey^{2}$), but available computing power still
quickly limits the maximum Reynolds number that can be obtained. The
simulations presented in Groom & Thornber (2021) represent the current state
of the art in terms of maximum Reynolds number that can be achieved using the
Flamenco algorithm. Even then, the highest Reynolds number simulation in that
study was still short of meeting the mixing transition requirement for fully
developed turbulence in unsteady flows (Zhou et al., 2003).

For DNS of broadband RMI, assuming the same grid resolution is used, the
larger bandwidth necessitates a smaller Reynolds number since the number of
cells per $\lambda_{min}$ required to resolve the shock-interface interaction
and subsequent evolution is the same. This is before any considerations about
whether additional grid resolution is required at later time due to increasing
Reynolds number. The requirement that all initial amplitudes be linear also
limits the initial velocity jump (and hence the Reynolds number) that can be
obtained, and the diffuse profile across the interface that is required to
properly resolve the shock-interface interaction in DNS also dampens the
initial velocity jump (relative to if a sharp interface was used). All of this
results in the fact that for the current maximum grid sizes simulated in this
and previous studies (e.g. $2048^{2}$ cross-sectional resolution), DNS can be
performed at either a moderate Reynolds number but small bandwidth (i.e. too
narrow to be indicative of real surface perturbations) as in Groom & Thornber
(2021) or a moderate bandwidth but low Reynolds number (i.e. too diffuse to be
indicative of fully-developed turbulence) as in the present study. These
observations are not exclusive to DNS of RMI but also apply to RTI, Kelvin-
Helmholtz instability and other flows where the effects of initial conditions
are important and realistic initial perturbations need to be considered.

In spite of all this, DNS is still a useful tool in the context of this study
as it provides results that may be considered a plausible lower bound to the
experimental results in a similar manner to which ILES results may be
considered a plausible upper bound. It is also necessary for computing
statistical quantities that depend on the smallest scales of motion being
sufficiently resolved, such as the turbulent length scales and Reynolds
numbers presented in §3.6 as well as many other quantities that are important
for informing modelling of these types of flows (see Groom & Thornber (2021);
Wong et al. (2022) for some examples). Comments on how some of the limitations
mentioned above might be resolved are given in §4.

##  3 Results

Using the initial conditions and computational setup described in §2, six
simulations are performed with Flamenco. These consist of four ILES
corresponding to the four different initial conditions as well as two DNS; one
for the $m=-1$ initial condition and one for the $m=-2$ initial condition. The
viscosity used in these DNS is $\overline{\mu}=0.3228$ Pa$\cdot$s, which
corresponds to initial Reynolds numbers of $\Rey_{0}=261$ and $\Rey_{0}=526$
for the $m=-1$ and $m=-2$ cases respectively. While this viscosity is much
higher than would occur experimentally, it is equivalent to using a much
smaller value of $\overline{\lambda}$ to obtain the same Reynolds number due
to the various simplifications employed in the governing equations, such as no
variation in viscosity with temperature. For each simulation, grid convergence
is assessed using the methodology outlined in Thornber et al. (2017) for ILES
and Groom & Thornber (2019) for DNS. The simulations were run up to a physical
time of $t=0.1$ s, at which point some of the spikes were observed to have
reached the domain boundaries in the $m=-3$ ILES case. The complete set of
simulations is summarised in table 3.

Case | $m$ | $\Rey_{0}$ | Simulation time (s) | Domain size (m3) | Grid resolution
---|---|---|---|---|---
1 | 0 | - | 0.1 | $1.5\upi\times 2\upi\times 2\upi$ | $384\times 512^{2}$
2 | -1 | - | 0.1 | $1.5\upi\times 2\upi\times 2\upi$ | $384\times 512^{2}$
3 | -2 | - | 0.1 | $1.5\upi\times 2\upi\times 2\upi$ | $384\times 512^{2}$
4 | -3 | - | 0.1 | $1.5\upi\times 2\upi\times 2\upi$ | $384\times 512^{2}$
5 | -1 | 261 | 0.1 | $0.75\upi\times 2\upi\times 2\upi$ | $384\times 1024^{2}$
6 | -2 | 526 | 0.1 | $0.75\upi\times 2\upi\times 2\upi$ | $384\times 1024^{2}$
Table 3: The initial power spectrum slope, initial Reynolds number (DNS only),
total simulation time, domain size and maximum grid resolution employed for
each case.

Figure 4 shows visualisations of the solution at the latest time of $t=0.1$ s
for the four ILES cases. Bubbles of light fluid can be seen flowing into the
heavy fluid on the lower side of the mixing layer, while heavy spikes are
penetrating into the light fluid on the upper side. In the narrowband case the
mixing layer has remained relatively uniform over the span of the domain,
whereas in the broadband cases, particularly the $m=-2$ and $m=-3$ cases,
large-scale entrainment is starting to occur at scales on the order of the
domain width. Another noticeable phenomenon at this time is that in the
narrowband case some spikes have penetrated much further away from the main
mixing layer than in the broadband cases. This is shown in greater detail in
figure 6 where isosurfaces of volume fraction $f_{1}=0.001$ and $f_{1}=0.999$
are plotted for both the $m=0$ narrowband case and the $m=-2$ broadband case
to highlight the differences in spike behaviour. Note that in the narrowband
case there are taller structures on the spike side that in some instances have
been ejected from the main layer. See also Figure 5 from Youngs & Thornber
(2020a) for a similar visualisation at a lower Atwood number. A plausible
explanation for this is that the slower but more persistent growth of the low
wavenumber modes in the broadband cases cause the main mixing layer to
eventually disrupt the trajectory of any spikes that were initially ejected
from high wavenumber modes. Future work will study this comparison of spike
behaviour between narrowband and broadband mixing perturbations at higher
Atwood numbers that are more relevant to ICF.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 4: Contours of volume fraction $f_{1}$ for the ILES cases at $t=0.1$ s
and $z=0$. The major ticks on both axes correspond to a grid spacing of
$\Delta
x={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
y=}}1$ m.

Figure 5 shows visualisations at the same physical time for the two DNS cases.
As discussed in §2.4, these DNS are at quite low Reynolds number so as to be
able to fully resolve the wide range of initial length scales. They are
therefore quite diffuse, however good agreement can still be observed in the
largest scales of motion with the corresponding ILES cases. The fluctuating
kinetic energy spectra presented in §3.5 also corroborate this observation.

(a) $m=-1$, $\Rey_{0}=261$.

(b) $m=-2$, $\Rey_{0}=526$.

Figure 5: Contours of volume fraction $f_{1}$ for the DNS cases at $t=0.1$ s
and $z=0$. The major ticks on both axes correspond to a grid spacing of
$\Delta
x={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta
y=}}1$ m.

Figure 6: Isosurfaces of volume fraction $f_{1}$ for the $m=0$ (left) and
$m=-2$ (right) ILES cases at $t=0.1$ s.

###  3.1 Non-dimensionalisation

The results in the following sections are appropriately non-dimensionalised to
allow for direct comparisons with the experiments in Jacobs et al. (2013) and
Sewell et al. (2021). All length scales are normalised by $\lambda_{min}$,
which is equal to $0.196$ m in the simulations and is estimated to lie between
$2.9$ mm and $3.2$ mm in the experiments. As the effects of different initial
impulses are of primary interest, it does not makes sense to use $\dot{W_{0}}$
as the normalising velocity scale, therefore all velocities are normalised by
$A^{+}\Delta u$ instead. In the simulations $A^{+}=0.72$ and $\Delta u=158.08$
m/s, while in the experiments $A^{+}=0.71$ and $\Delta u=74$ m/s. Therefore
the non-dimensional time is given by

$\tau=\frac{(t-t_{0})A^{+}\Delta u}{\lambda_{min}}$ (15)

where $t_{0}=0.0011$ s is the shock arrival time. This equates to a
dimensionless time of $\tau=57.4$ at the latest time considered in the
simulations ($t=0.1$ s), $107\leq\tau\leq 118$ at the latest time prior to
reshock in the experiments of Jacobs et al. (2013) ($t-t_{0}=6.5$ ms) and
$73.9\leq\tau\leq 81.5$ at the latest time prior to reshock in the experiments
of Sewell et al. (2021) ($t-t_{0}=4.5$ ms), assuming the same range of values
for $\lambda_{min}$ of $2.9$ to $3.2$ mm. Figure 7 shows a subset of the image
sequence taken from a typical vertical shock tube experiment in Jacobs et al.
(2013) using the Mie diagnostic . For comparison with the present simulations,
a dimensionless time of $\tau=57.4$ corresponds to a physical time in the
range of $t=3.17$ ms to $t=3.50$ ms, which may be compared with the images
shown for times $t=3.00$ ms and $t=3.50$ ms in figure 7.

Figure 7: Image sequence taken from a typical vertical shock tube experiment
using the Mie diagnostic. Times relative to shock impact are shown in each
image. Reshock occurs at $t=6.50$ ms. Source: Figure 3 of Jacobs et al.
(2013).

###  3.2 Turbulent Kinetic Energy and Mix Width

In this section comparisons are made both between the present simulation
results and those of the experiments, as well as between the methods for
calculating those results in the experiments with methods that have been
commonly employed in previous simulation studies of RMI. To measure the mixing
layer width, Jacobs et al. (2013) used Mie scattering over a single plane,
with each image then row-averaged to obtain the mean smoke concentration in
the streamwise direction. For each concentration profile, the mixing layer
width is defined as the distance between the 10% and 90% threshold locations.
This is similar to the definition of visual width used in simulation studies
of both RMI and RTI (see Cook & Dimotakis, 2001; Cook & Zhou, 2002; Zhou &
Cabot, 2019), where the plane-averaged mole fraction or volume fraction
profile is used along with a typical threshold cuttoff of 1% and 99%, e.g.

$h=x\left(\langle f_{1}\rangle=0.01\right)-x\left(\langle
f_{1}\rangle=0.99\right).$ (16)

This is a useful definition of the outer length scale of the mixing layer,
however the choice of cutoff location is somewhat arbitrary and when used to
estimate growth rates the results are influenced by both the the choice of
cutoff location as well as statistical fluctuations (Zhou & Cabot, 2019). For
that purpose, an integral definition is typically used such as the integral
width (Andrews & Spalding, 1990)

$W=\int\langle f_{1}\rangle\langle f_{2}\rangle\>\mathrm{d}x.$ (17)

If $f_{1}$ varies linearly with $x$ then $h=6W$ (Youngs, 1994). See also the
recent paper by Youngs & Thornber (2020a) where integral definitions of the
bubble and spike heights are proposed that are of similar magnitude to the
visual width. These are presented in Appendix B and are discussed in §3.3
below.

In the experiments of Sewell et al. (2021), PIV was used as the main
diagnostic and therefore an alternate definition of the mixing layer width was
required. In that study, the row-averaged turbulent kinetic energy was used
and a mixing layer width defined as the distance between the $x$-locations at
which the TKE is 5% of its peak value. This definition assumes that the
turbulent velocity field spreads at the same rate as the mixing layer. Figure
8 shows streamwise profiles of mean turbulent kinetic energy for each of the
four initial conditions, defined as

$\mathrm{TKE}=\frac{1}{2}\overline{u_{i}^{\prime}u_{i}^{\prime}}$ (18)

where $\psi^{\prime}=\psi-\overline{\psi}$ indicates a fluctuating quantity
and the ensemble average $\overline{\psi}=\langle\psi\rangle$ is calculated as
a plane average taken over the statistically homogeneous directions (in this
case $y$ and $z$). The volume fraction profile $\langle f_{1}\rangle\langle
f_{2}\rangle$ is also shown on the right axis of each plot, as well as the
(outermost) $x$-locations at which the TKE is 5% of its peak value. An
important feature worth noting when comparing the narrowband case with the
other broadband cases is that the 5% cutoff on the spike side ($x<x_{c}$) is
further from the mixing layer centre $x_{c}$ than in the $m=-1$ and $m=-2$
cases, despite these cases having a greater overall amplitude in the initial
perturbation. There is also a greater amount of mixed material, as measured by
the product $\langle f_{1}\rangle\langle f_{2}\rangle$, at this location than
in those two broadband cases, which is in line with the observations made in
figure 4 about the greater penetration distances of spikes from the main layer
in the narrowband case. In all cases the TKE profile is asymmetric, with the
5% cutoff on the spike side being located further away from the mixing layer
centre than the corresponding 5% cutoff on the bubble side. This asymmetry,
along with the implications it has for the growth rate exponent $\theta$, is
discussed in further detail §3.3.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 8: Plane-averaged profiles of $\mathrm{TKE}$ for each initial condition
at time $\tau=57.4$. Solid lines indicate ILES results and dotted lines
indicate DNS results. Also shown are the volume fraction profiles (grey solid
lines), along with the 5% cutoff locations (crosses) and the TKE centroid
(black dashed lines).

In Sewell et al. (2021) a definition for the mixing layer centre is given as
the centroid of the mean turbulent kinetic energy profile, i.e.

$x_{c}=\frac{\displaystyle\int xf(x)\>\mathrm{d}x}{\displaystyle\int
f(x)\>\mathrm{d}x}$ (19)

where $f(x)$ is the mean turbulent kinetic energy profile. This centroid is
also shown in figure 8. This definition is compared with an alternate
definition in terms of the $x$-location of equal mixed volumes,

$\int_{-\infty}^{x_{c}}\langle
f_{2}\rangle\>\mathrm{d}x=\int_{x_{c}}^{\infty}\langle
f_{1}\rangle\>\mathrm{d}x$ (20)

which has been used previously in both computational (Walchli & Thornber,
2017; Groom & Thornber, 2021) and experimental (Krivets et al., 2017) studies
of RMI. Figure 9 plots the temporal evolution of both of these definitions for
$x_{c}$ for each initial condition, showing that the TKE centroid consistently
drifts towards the spike side of the layer as time progresses. The definition
in terms of position of equal mixed volumes is much more robust and remains
virtually constant throughout the simulation. There is also little variation
between cases for this definition, unlike the TKE centroid which is more
biased towards the spike side in the $m=-3$ and $m=0$ cases. The choice of
definition for the mixing layer centre is important as it will influence the
bubble and spike heights that are based off it (as well as their ratio), along
with any quantities that are plotted at the mixing layer centre over various
points in time.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 9: Temporal evolution of the mixing layer centre $x_{c}$, comparing the
definition based on the centroid of the mean turbulent kinetic energy profile
with the definition based on the $x$-location of equal mixed volumes.

Figure 10 shows the temporal evolution of the mixing layer width, using both
the visual width definition based on the mean volume fraction profile
(referred to as the VF-based width) as well as the definition from Sewell et
al. (2021) based on the distance between the 5% cutoff locations in the mean
turbulent kinetic energy profile (referred to as the TKE-based width). The
mean volume fraction $f_{1}$ at these 5% cutoff locations is $\geq 0.997$ on
the spike side ($x<x_{c}$) and $\leq 0.003$ on the bubble side ($x>x_{c}$) in
all cases, hence why the TKE-based width is larger than the VF-based width in
each of the plots as the VF-based width is defined using a 1% and 99% cutoff
in the volume fraction profile. Using nonlinear regression to fit a function
of the form $h=\beta(\tau-\tau_{0})^{\theta}$, the growth rate exponent
$\theta$ can be obtained for the TKE-based width, VF-based width and the
integral width (not shown in figure 10) for each case. Following Sewell et al.
(2021), the fit is performed only for times satisfying
$\bar{k}\dot{\sigma_{0}}t>1$ so that the flow is sufficiently developed. The
estimated value of $\theta$ for each case is given in table 4. Note that the
uncertainties reported are merely taken from the variance of the curve-fit and
do not represent uncertainties in the true value of $\theta$.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 10: Temporal evolution of mixing layer width $h$ based on the distance
between cutoff locations using either the mean turbulent kinetic energy or
mean volume fraction profiles. Solid lines indicate ILES results and dotted
lines indicate DNS results. Curve-fits to the data are also shown, with the
relevant data points used given by the symbols in each plot.

Analysing the results in table 4, there is good agreement between the values
of $\theta$ obtained from the visual and integral widths for all cases. This
is mainly a verification that the results are not severely impacted by a lack
of statistical resolution at the lowest wavenumbers, which would result in the
visual width measurements being dependent on the specific realisation. The
small differences in the values of $\theta$ reported indicate that there is
still some influence of statistical fluctuations, therefore the estimates made
using the integral width should be regarded as the most accurate. When
comparing the TKE-based and VF-based threshold widths, there is good agreement
for the broadband ILES cases and in particular for the $m=-3$ ILES case. For
the narrowband ILES case however, the VF-based (and integral) width is growing
at close to the theoretical value of $\theta=1/3$ for self-similar decay
proposed by Elbaz & Shvarts (2018), whereas the TKE-based width is growing at
a much fast rate of $\theta=0.589$. This is even faster than any of the
broadband cases and is due to the sensitivity of the TKE-based width to spikes
located far from the mixing layer centre in the narrowband case, which contain
very little material but are quite energetic and which grow at a faster rate
than the rest of the mixing layer. For the broadband DNS, the growth rate of
the TKE-based width is slightly lower than that of the VF-based width for both
cases, indicating that turbulent fluctuations are more confined to the core of
the mixing layer. In the $m=-1$ case, the value of $\theta$ obtained from the
integral width is Reynolds number independent, while for $m=-2$ the value of
$\theta$ obtained from the integral width in the DNS case is converging
towards the high Reynolds number limit given by the ILES case. Given that the
broadband perturbations, specifically the $m=-3$ perturbation, are the most
relevant to the experiments in Jacobs et al. (2013) and Sewell et al. (2021),
it is reassuring to note that estimates of $\theta$ made using TKE-based
widths measured with PIV correspond well with estimates based off the
concentration field.

Case | $m$ | $\Rey_{0}$ | TKE-based width $\theta$ | VF-based width $\theta$ | Integral width $\theta$ | TKE decay rate $\theta$
---|---|---|---|---|---|---
1 | 0 | - | $0.589\pm 1.20\times 10^{-2}$ | $0.323\pm 4.89\times 10^{-3}$ | $0.330\pm 1.27\times 10^{-3}$ | $0.253\pm 7.00\times 10^{-3}$
2 | -1 | - | $0.460\pm 1.03\times 10^{-2}$ | $0.450\pm 1.54\times 10^{-3}$ | $0.442\pm 1.10\times 10^{-4}$ | $0.429\pm 5.65\times 10^{-3}$
3 | -2 | - | $0.479\pm 3.92\times 10^{-3}$ | $0.522\pm 3.59\times 10^{-3}$ | $0.514\pm 3.60\times 10^{-4}$ | $0.512\pm 3.47\times 10^{-3}$
4 | -3 | - | $0.493\pm 6.25\times 10^{-3}$ | $0.492\pm 1.25\times 10^{-3}$ | $0.510\pm 1.91\times 10^{-3}$ | $0.562\pm 2.22\times 10^{-3}$
5 | -1 | 261 | $0.444\pm 1.41\times 10^{-2}$ | $0.501\pm 8.40\times 10^{-4}$ | $0.441\pm 1.00\times 10^{-4}$ | $0.492\pm 8.08\times 10^{-3}$
6 | -2 | 526 | $0.456\pm 3.69\times 10^{-3}$ | $0.556\pm 2.27\times 10^{-3}$ | $0.549\pm 1.52\times 10^{-3}$ | $0.576\pm 4.69\times 10^{-3}$
Table 4: Estimates of the growth rate exponent $\theta$ from curve-fits to the
TKE-based, VF-based and integral widths, as well as from the decay rate of
total turbulent kinetic energy.

An alternative method for estimating $\theta$ is also given in Sewell et al.
(2021), which makes use of the decay rate of total fluctuating kinetic energy
and a relationship between this decay rate $n$ and the mixing layer growth
rate $\theta$ originally derived by Thornber et al. (2010). Assuming that
$h\propto t^{\theta}$ and the mean fluctuating kinetic energy
$q_{k}\propto\dot{h}^{2}$ gives the relation $q_{k}\propto t^{2\theta-2}$.
Since the total fluctuating kinetic energy is proportional to the width of the
mixing layer multiplied by the mean fluctuating kinetic energy, this gives
$\mathrm{TKE}\propto t^{3\theta-2}\propto t^{n}$. Directly measuring the decay
rate $n$ therefore gives an alternative method for estimating $\theta$, which
is particularly useful in experimental settings where only velocity field data
is available. This predicted value of $\theta=(n+2)/3$ has been found to be in
good agreement with the measured growth rate from the integral width in
multiple studies of narrowband RMI Thornber et al. (2010, 2017). However,
Groom & Thornber (2020) showed that for RMI evolving from broadband
perturbations with bandwidths as large as $R=128$ the measured values of
$\theta$ do not agree with this theoretical prediction, indicating that longer
periods of growth dominated by just-saturating modes are required than can
currently be obtained in simulations. Figure 11 shows the temporal evolution
of TKE, where the integration has been performed between the 5% cutoff
locations used to define the TKE-based width. Nonlinear least squares
regression is again used to estimate $n$ for each case, with the fit performed
for times greater than the point at which the curvature becomes convex. The
corresponding value of $\theta$ for each $n$ using the relation $n=3\theta-2$
is given in table 4.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 11: Temporal evolution of total fluctuating kinetic energy, integrated
between the 5% cutoff locations. Solid lines indicate ILES results and dotted
lines indicate DNS results. Curve-fits to the data are also shown, with the
relevant data points used given by the symbols in each plot.

For the narrowband case the estimate of $\theta$ from the TKE decay rate does
not agree with the other estimates, indicating that the mixing layer growth is
not sufficiently self-similar (a key assumption in the derivation) and lags
the decay in TKE. This is still true even when the range of times used in the
curve-fitting procedure is restricted to be the same as for the curve-fit to
the decay rate (not shown). For the broadband cases there is better agreement
however, particularly in the $m=-1$ and $m=-2$ ILES cases. In all broadband
cases the bandwidth of the initial perturbation is relatively small compared
to the perturbations analysed in Groom & Thornber (2020) and the longest
initial wavelength saturates early on in the overall simulation, therefore the
conclusions made in that study regarding the $n=3\theta-2$ relation do not
necessarily apply here as the current broadband cases are not in the self-
similar growth regime. They are also likely not in full self-similar decay
however, especially if the narrowband case is not, yet the values of $\theta$
are in better agreement than in the narrowband case. Further work is required
to determine why this is indeed the case.

Comparing the estimates of $\theta$ with those in Sewell et al. (2021) using
both the TKE-based width and TKE decay rate, the $m=-3$ simulation results are
in between the results of the low-amplitude and high-amplitude experiments.
For the low-amplitude experiments (prior to reshock), the TKE-based width
measurements gave $\theta=0.45$ and the TKE decay rate measurements gave
$\theta=0.68$ (which would correspond to no decay of TKE if the layer was
homogeneous (Barenblatt et al., 1983)). The equivalent results in the $m=-3$
simulation were $\theta=0.493$ and $\theta=0.562$, i.e. larger and smaller
than the respective experimental results but both within the experimental
margins of error. Similarly for the high-amplitude experiments, both the TKE-
based width measurements and the TKE decay rate measurements gave
$\theta=0.51$, indicating that the turbulence in the mixing layer is more
developed and closer to self-similar prior to reshock. The $m=-3$ simulation
results are also within the experimental margins of error for these results.
Overall, the combination of experimental and computational evidence indicates
that there are persistent effects of initial conditions when broadband surface
perturbations are present for a much greater period of time than just the time
to saturation of the longest initial wavelength (as considered in previous
simulation studies of broadband RMI) and last for the duration of the first-
shock growth in a typical shock tube experiment. Furthermore, a consideration
of the impact of finite bandwidth in the initial power spectrum (also referred
to as confinement) is required when adapting theoretical results for infinite
bandwidth (unconfined, see Youngs (2004); Thornber et al. (2010); Soulard et
al. (2018); Soulard & Griffond (2022)) to a specific application.

###  3.3 Bubble and Spike Heights

In order to help better explain the estimates for $\theta$ given in table 4,
it is useful to decompose the TKE-based and VF-based widths into separate
bubble and spike heights, $h_{b}$ and $h_{s}$, defined as the distance from
the mixing layer centre $x_{c}$ to the relevant cutoff location on the bubble
and spike side of the layer respectively. Given the drift in time for the
centroid of the TKE profile shown in figure 9, the $x$-location of equal mixed
volumes is used as the definition of the mixing layer centre for both the VF-
based and TKE-based bubble and spike heights. Figures 12 and 13 show the
evolution in time of $h_{b}$ and $h_{s}$ respectively for heights based off
both the 5% TKE cutoff (referred to as TKE-based heights) and the 1% and 99%
volume fraction cutoff (referred to as VF-based heights).

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 12: Temporal evolution of the bubble height $h_{b}$ based on the
distance between cutoff locations using either the mean turbulent kinetic
energy or mean volume fraction profiles. Solid lines indicate ILES results and
dotted lines indicate DNS results. Curve-fits to the data are also shown, with
the relevant data points used given by the symbols in each plot.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 13: Temporal evolution of the spike height $h_{s}$ based on the
distance between cutoff locations using either the mean turbulent kinetic
energy or mean volume fraction profiles. Solid lines indicate ILES results and
dotted lines indicate DNS results. Curve-fits to the data are also shown, with
the relevant data points used given by the symbols in each plot.

Some important trends can be observed. Firstly, the VF-based heights are
smoother than the corresponding TKE-based heights indicating that they are
less sensitive to statistical fluctuations. Secondly, the TKE-based $h_{b}$
and $h_{s}$ are greater than the corresponding VF-based heights in all cases
and for both measures the spike height is greater than the bubble height. This
can also be seen in figure 14, which plots the ratio $h_{s}/h_{b}$ vs. time
and shows that $h_{s}/h_{b}>1$ for all cases. The same trend was observed in
Youngs & Thornber (2020a) for both $At=0.5$ and $At=0.9$ but in a heavy-light
configuration where the heavy spikes are being driven into the lighter fluid
in the same direction as the shock wave. Appendix B plots the same integral
definitions of the bubble and spike heights used in Youngs & Thornber (2020a),
verifying that the behaviour is very similar to the VF-based heights presented
here. The ratio of spike to bubble heights using both threshold measures is
also very similar at late time in all cases with the exception of the
narrowband case. The ratio $h_{s}/h_{b}$ also appears to be converging to the
same value at late-time in all cases except for the TKE-based heights in the
narrowband case, suggesting it is only dependent on the Atwood number.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 14: Temporal evolution of the ratio of spike to bubble height. Solid lines indicate ILES results and dotted lines indicate DNS results. Curve-fits to the data are also shown, with the relevant data points used given by the symbols in each plot. Case | $m$ | $\Rey_{0}$ | TKE-based $\theta_{b}$ | VF-based $\theta_{b}$ | TKE-based $\theta_{s}$ | VF-based $\theta_{s}$
---|---|---|---|---|---|---
1 | 0 | - | $0.493\pm 2.43\times 10^{-2}$ | $0.441\pm 2.43\times 10^{-2}$ | $0.615\pm 1.72\times 10^{-2}$ | $0.277\pm 5.38\times 10^{-3}$
2 | -1 | - | $0.350\pm 8.94\times 10^{-3}$ | $0.514\pm 1.04\times 10^{-3}$ | $0.509\pm 1.57\times 10^{-2}$ | $0.425\pm 2.46\times 10^{-3}$
3 | -2 | - | $0.355\pm 1.31\times 10^{-2}$ | $0.466\pm 4.29\times 10^{-3}$ | $0.543\pm 8.58\times 10^{-3}$ | $0.550\pm 3.31\times 10^{-3}$
4 | -3 | - | $0.280\pm 2.95\times 10^{-2}$ | $0.282\pm 1.79\times 10^{-3}$ | $0.586\pm 1.07\times 10^{-2}$ | $0.606\pm 1.49\times 10^{-3}$
5 | -1 | 261 | $0.338\pm 1.36\times 10^{-2}$ | $0.461\pm 3.20\times 10^{-4}$ | $0.509\pm 1.61\times 10^{-2}$ | $0.523\pm 1.46\times 10^{-3}$
6 | -2 | 526 | $0.284\pm 8.39\times 10^{-3}$ | $0.458\pm 2.57\times 10^{-3}$ | $0.561\pm 4.89\times 10^{-3}$ | $0.613\pm 2.23\times 10^{-3}$
Table 5: Estimates of the growth rate exponents $\theta_{b}$ and $\theta_{s}$
from curve-fits to the TKE-based and VF-based bubble and spike heights.

Figure 14 shows that the ratio of $h_{s}/h_{b}$ is approximately constant by
the end of the simulations. This indicates that a single $\theta$ is
appropriate for describing the growth of the mixing layer beyond this point.
However, prior to that $h_{b}$ and $h_{s}$ do grow at different rates as shown
in table 5, where the bubble growth rate exponent is denoted by $\theta_{b}$
and the spike growth rate exponent is denoted by $\theta_{s}$. Two key trends
can be observed; the VF-based $\theta_{b}$ is greater than the TKE-based
$\theta_{b}$ in all cases other than the narrowband ($m=0$) case, while the
VF-based $\theta_{s}$ is greater than the TKE-based $\theta_{s}$ in all cases
other than the $m=-1$ ILES case and the narrowband case. The $m=-3$ case also
has the smallest difference in $\theta_{b}$ and $\theta_{s}$ for both
threshold measures. Comparing the DNS cases with their respective ILES cases,
the VF-based $h_{b}$ is almost independent of the Reynolds number in both the
$m=-1$ and $m=-2$ cases. This is also true for the TKE-based $h_{s}$ in the
$m=-2$ cases. A higher degree of Reynolds number dependence is observed for
both definitions of $h_{s}$, which is consistent with previous observations
made about turbulence developing preferentially on the spike side of the
mixing layer Groom & Thornber (2021). This can also be observed for the
integral definitions of $h_{b}$ and $h_{s}$ given in Appendix B.

This analysis provides evidence that, prior to reshock, $h_{b}$ and $h_{s}$ do
grow at different rates in a typical shock tube experiment. However, their
growth rate exponents have equalised by the time reshock arrives. This is a
complicating factor when estimating a single value for $\theta$ at early times
and points to the difficulties in obtaining self-similar growth for RMI in
both experiments and simulations. This also suggests that the ratio of spike
to bubble heights could be used to determine when it is appropriate to start
curve-fitting for estimating a single value of $\theta$, and that measurements
based on the concentration field are likely more accurate in this regard than
those made using the velocity field.

###  3.4 Anisotropy

The anisotropy of the fluctuating velocity field is explored using the same
two measures presented in Sewell et al. (2021). The first is a global measure
of anisotropy, defined as

$\mathrm{TKR}=\frac{2\times\mathrm{TKX}}{\mathrm{TKY}+\mathrm{TKZ}}$ (21)

where $\mathrm{TKX}=\frac{1}{2}\overline{u^{\prime}u^{\prime}}$,
$\mathrm{TKY}=\frac{1}{2}\overline{v^{\prime}v^{\prime}}$ and
$\mathrm{TKZ}=\frac{1}{2}\overline{w^{\prime}w^{\prime}}$, with each quantity
integrated between the cutoff locations based on 5% of the maximum TKE. The
second measure is the Reynolds stress anisotropy tensor, whose components are
defined by

$b_{ij}=\frac{\overline{u_{i}^{\prime}u_{j}^{\prime}}}{\frac{1}{2}\overline{u_{i}^{\prime}u_{i}^{\prime}}}-\frac{1}{3}\delta_{ij}.$
(22)

This tensor, specifically the $x$-direction principal component
$\boldsymbol{b}_{11}$ for this particular flow, is a measure of anisotropy in
the energy-containing scales of the fluctuating velocity field with a value of
0 indicating isotropy in the direction of that component. The local version of
TKR (i.e. with TKX, TKY and TKZ not integrated in the $x$-direction) can be
written in terms of $\boldsymbol{b}_{11}$ as

$\frac{2\overline{u^{\prime}u^{\prime}}}{\overline{v^{\prime}v^{\prime}}+\overline{w^{\prime}w^{\prime}}}=\frac{2\boldsymbol{b}_{11}+2/3}{2/3-\boldsymbol{b}_{11}}$
(23)

allowing the two measures to be related to one another.

Figure 15 shows the temporal evolution of the global anisotropy measure TKR
for each case. Compared to the equivalent figure 13 in Sewell et al. (2021)
the peak in anisotropy at early time is less pronounced, however this is due
to only integrating TKX, TKY and TKZ between the 5% cutoff locations. Figure
10 in Groom & Thornber (2019) shows the same measure without this limit on the
integration for a similar case, with the peak in anisotropy much closer to
that observed in Sewell et al. (2021). This indicates that much of the
anisotropy observed at very early times is due to the shock wave. At an
equivalent dimensionless time to the latest time simulated here, the
anisotropy ratio presented in Sewell et al. (2021) is approximately 2 for the
high-amplitude experiments and 3 for the low-amplitude experiments. For the
$m=-3$ perturbation that most closely matches those experiments the TKR at the
latest time is 2.46, while for the other ILES cases the late-time TKR
decreases as $m$ increases. For the $m=0$ narrowband case the late-time value
is 1.55, which is within the range of 1.49-1.66 observed across codes on the
$\theta$-group quarter-scale case (Thornber et al., 2017); a case which is
essentially the same perturbation but at a lower Atwood number. For the DNS
cases a very different trend is observed where the anisotropy continually
grows as time progresses. This is due to the very low Reynolds numbers of
these simulations, with the lack of turbulence preventing energy from being
transferred to the transverse directions.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 15: Temporal evolution of the global anisotropy measure, with each
component integrated between the 5% cutoff locations. Solid lines indicate
ILES results and dotted lines indicate DNS results.

The spatial variation in anisotropy is shown in figure 16, plotted between the
5% cutoff locations for each case. For the broadband cases the anisotropy is
slightly higher on the spike side of the layer, with the greatest increase in
the $m=-3$ case. This mirrors the results shown in Sewell et al. (2021) for
$\boldsymbol{b}_{11}$, with quite good agreement observed between the $m=-3$
case at the latest time and the low-amplitude experiments just prior to
reshock. In the narrowband case the increase in anisotropy from the mixing
layer centre to the spike side is greater but the overall magnitude of
$\boldsymbol{b}_{11}$ is lower, consistent with what was observed for TKR. The
DNS results show that the biggest increase in anisotropy at low Reynolds
numbers is in the centre of the mixing layer; there is a smaller difference in
anisotropy between the DNS and ILES cases at either edge.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 16: Spatial distribution of the $x$-direction principal component of
the Reynolds stress anisotropy tensor at time $\tau=57.4$. Solid lines
indicate ILES results and dotted lines indicate DNS results. Also shown is the
mixing layer centre defined by the TKE centroid (black dashed lines).

Figure 17 shows the temporal evolution of $\boldsymbol{b}_{11}$ at the mixing
layer centre, both for the definition of $x_{c}$ in terms of the TKE centroid
(shown in figure 16) as well as the alternate definition in terms of the
position of equal mixed volumes. The results for both definitions are similar
across all cases, with the anisotropy at the position of equal mix being
slightly lower in all cases. In the DNS cases $\boldsymbol{b}_{11}$ is
approximately constant in time, indicating that the growth in anisotropy that
was observed for TKR in figure 15 is occurring on either side of the mixing
layer centre. The range of values are also comparable to those given in Wong
et al. (2019) prior to reshock.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 17: Temporal evolution of $x$-direction principal component of the
Reynolds stress anisotropy tensor at the mixing layer centre plane. Solid
lines indicate ILES results and dotted lines indicate DNS results.

###  3.5 Spectra

The distribution of fluctuating kinetic energy per unit mass across the
different scales of motion is examined using radial power spectra of the
transverse and normal components, calculated as

$E_{i}(\kappa)=\widehat{u^{\prime}_{i}}^{\dagger}\widehat{u^{\prime}_{i}}$
(24)

where $\kappa=\sqrt{\kappa_{y}+\kappa_{z}}$ is the (dimensionless) radial
wavenumber in the $y$-$z$ plane at $x=x_{c}$ (given by the $x$-location of
equal mixed volumes), $\widehat{(\ldots)}$ denotes the 2D Fourier transform
taken over this plane and $\widehat{(\ldots)}^{\dagger}$ is the complex
conjugate of this transform. As isotropy is expected in the transverse
directions, the $E_{y}(\kappa)$ and $E_{z}(\kappa)$ spectra are averaged to
give a single transverse spectrum $E_{yz}(\kappa)$.

The normal and transverse spectra are shown in figure 18 for each of the ILES
and DNS cases at the latest simulated time. Curve-fits are made to the data to
determine the scaling of each spectrum, with some interesting trends observed.
For broadband cases evolving from perturbations of the form given in (8), a
scaling of $E(\kappa)\sim\kappa^{(m+2)/2}$ is expected for the low wavenumbers
at early time while the growth of the mixing layer is being dominated by the
just-saturating mode (Groom & Thornber, 2020). This is not observed in figure
18 since saturation of the longest wavelength occurs quite early relative to
the end time of the simulations, however some lingering effects can still be
seen at the lowest wavenumbers. For all three broadband ILES cases there are
two distinct ranges in both the normal and transverse spectra, which
approximately correspond to wavenumbers lower and higher than
$\kappa_{max}=k_{max}(L/2\upi)=32$. Thornber et al. (2010) modified the
analysis of Zhou (2001) to take into account the effects of the initial
perturbation spectrum, resulting in an expected scaling for broadband
perturbations of the form $E(\kappa)\sim\kappa^{(m-6)/4}$. This scaling is
observed for the transverse spectra at wavenumbers greater than
$\kappa_{max}$, while for the normal spectra a scaling of
$E(\kappa)\sim\kappa^{(m-5)/4}$ is observed, the reason for which is currently
unclear.

For wavenumbers less than $\kappa_{max}$ the normal spectra scale as
$\kappa^{-3/2}$ in the $m=-2$ and $m=-3$ cases, which is in good agreement
with previous calculations for narrowband perturbations (Thornber, 2016; Groom
& Thornber, 2019). The narrowband case presented here has a slightly less
steep scaling for both the normal and transverse spectra, although it has not
been run to as late of a dimensionless time as in previous studies such as
Thornber et al. (2017). The normal spectrum in the $m=-1$ case also has a
scaling that is less steep than $\kappa^{-3/2}$. A possible explanation for
this is that saturation occurs a lot later in this case than the other
broadband cases and therefore it may still be transitioning between an
$E(\kappa)\sim\kappa^{(m+2)/2}$ and a $\kappa^{-3/2}$ scaling. For the
transverse spectra in each of the broadband cases at wavenumbers less than
$\kappa_{max}$ a similar trend is observed, with each spectrum having a
scaling that is shallower than $\kappa^{-3/2}$. The same argument of
transition between an $E(\kappa)\sim\kappa^{(m+2)/2}$ and a $\kappa^{-3/2}$
scaling may also be applied here, however simulations to later time would be
required to confirm this.

Finally, for the DNS cases no inertial range is observed due to the low
Reynolds numbers that are simulated. For the normal spectra there is quite
good agreement between the DNS and ILES data in the energy-containing scales
at low wavenumbers. The transverse spectra contain less energy at these
wavenumbers in the DNS cases due to suppression of secondary instabilities
that transfer energy from the normal to transverse directions. Sewell et al.
(2021) did not observe an inertial range in their TKE spectra prior to
reshock, however they noted that there is likely some attenuation of the
spectra at scales smaller than the effective window size of their PIV method,
which is equivalent to a dimensionless wavenumber of $\kappa=47$. This makes
it difficult to compare and verify the current findings with their existing
experimental setup.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$.

(d) $m=0$ ($R=2$).

Figure 18: Transverse and normal components of fluctuating kinetic energy per
unit mass at the mixing layer centre plane at time $\tau=57.4$. Solid lines
indicate ILES results and dotted lines indicate DNS results.

###  3.6 Turbulent Length Scales and Reynolds Numbers

In order to give a better indication of how the present set of results compare
with the experiments of Jacobs et al. (2013) and Sewell et al. (2021), the
outer-scale Reynolds numbers and key turbulent length scales used to evaluate
whether a flow has transitioned to turbulence are computed using the DNS data.
For the purposes of comparison, both the TKE-based and VF-based threshold
widths are used as the outer length scale $h$ from which to compute the outer-
scale Reynolds number as

$\Rey_{h}=\frac{\overline{\rho^{+}}h\dot{h}}{\overline{\mu}}.$ (25)

Figure 19 shows the temporal variation for both definitions of the outer-scale
Reynolds number. The outer-scale Reynolds numbers using the TKE-based
definition for $h$ are roughly a factor of 2 larger, mostly due to the TKE-
based width being a lot larger than the VF-based width in all cases, with
neither definition close to reaching the critical value of $\Rey_{h}\gtrsim
1$-$2\times 10^{4}$ for fully developed turbulence (Dimotakis, 2000). For both
the $m=-1$ and $m=-2$ perturbations the VF-based Reynolds number is
approximately constant in time, consistent with the measured values of
$\theta$ given in table 4.

(a) $m=-1$.

(b) $m=-2$.

Figure 19: Outer-scale Reynolds numbers vs. time.

Dimotakis (2000) showed that for stationary flows, fully developed turbulence
is obtained when $\lambda_{L}/\lambda_{V}\geq 1$ where
$\lambda_{L}=5\lambda_{T}$ is the Liepmann–Taylor length scale and
$\lambda_{V}=50\lambda_{K}$ is the inner-viscous length scale, with
$\lambda_{T}$ and $\lambda_{K}$ the Taylor and Kolmogorov length scales
respectively. These length scales may be related to the outer-scale Reynolds
number by

$\displaystyle\lambda_{L}$ $\displaystyle=5\Rey_{h}^{-1/2}h$ (26)
$\displaystyle\lambda_{V}$ $\displaystyle=50\Rey_{h}^{-3/4}h$ (27)

from which it can be shown that $\Rey_{h}\geq 10^{4}$ for fully developed
turbulence. For a time-dependent flow, Zhou et al. (2003) showed that an
additional length scale $\lambda_{D}=5(\nu t)^{1/2}$ that characterises the
growth rate of shear-generated vorticity must be considered, referred to as
the diffusion layer scale. The condition for fully developed turbulence then
becomes

$\textrm{min}(\lambda_{L},\lambda_{D})>\lambda_{V}.$ (28)

Figure 20 shows the temporal variation of each length scale in (28), with
$\lambda_{L}$ and $\lambda_{V}$ calculated from the outer-scale Reynolds
number using both definitions for $h$. In both cases there is good agreement
between the length scales calculated from either definition of $\Rey_{h}$. The
inner-viscous length scale is greater than the Liepmann–Taylor scale at all
times in both cases, consistent with other observations in this paper on the
lack of fully developed turbulence in the DNS cases at the Reynolds numbers
capable of being simulated currently.

(a) $m=-1$.

(b) $m=-2$.

Figure 20: The Liepmann–Taylor (circles), inner-viscous (squares) and
diffusion length scales vs. time for both definitions of the outer-scale
Reynolds number.

Sewell et al. (2021) also observed $\lambda_{L}<\lambda_{V}$ at all times
prior to reshock in their low-amplitude experiments. The authors note that,
because of the different dependence of each length scale on $\Rey_{h}$, for
$\theta\leq 0.5$ the flow can never transition to turbulence as $\lambda_{V}$
will grow faster than $\lambda_{D}$. Furthermore, the definition for
$\lambda_{D}$ implies that it will be 0 at time $t=0$, which would seem to
imply that an RMI-induced flow with $\theta\leq 0.5$ can never become
turbulent. However, the virtual time origin is neglected in the original
definition for $\lambda_{D}$; if it is included then this allows for the
possibility that $\lambda_{V}<\lambda_{D}$ at early time. In that situation,
transition to turbulence will occur provided the initial velocity jump is
strong enough to produce $\lambda_{L}>\lambda_{V}$ for some period of time.
The turbulence will still be decaying over time if $\theta\leq 0.5$ though and
will eventually no longer be fully developed, reflecting a fundamental
difficulty to obtaining universal behaviour in experiments or numerical
simulations of RMI.

##  4 Conclusions

This paper has presented simulations of an idealised shock tube experiment
between air and sulphur hexafluoride that builds upon the previous results and
analysis presented in Groom & Thornber (2020, 2021). In particular, the
effects of additional long wavelength modes in the initial perturbation were
explored by comparing the results obtained using a narrowband surface
perturbation (similar to the one presented in Groom & Thornber (2021)) and
three broadband perturbations (similar to those presented in Groom & Thornber
(2020)). Both implicit large-eddy simulations (ILES) of the high-Reynolds
number limit as well as direct numerical simulations (DNS) at Reynolds numbers
lower than those observed in the experiments were performed with the Flamenco
finite-volume code.

Various measures of the mixing layer width, based on both the plane-averaged
turbulent kinetic energy and volume fraction profiles, were compared in order
to explore the effects of initial conditions as well as the validity of using
measurements based on the velocity field to draw conclusions about the
concentration field (and vice versa) as is commonly done in experiments due to
the difficulties of using diagnostics for both fields simultaneously. The
effects of initial conditions on the growth rate exponent $\theta$ were
analysed by curve-fitting the expected power law behaviour for the mixing
layer width $h$ to two different definitions of $h$; one based on a threshold
of 5% of the peak turbulent kinetic energy (TKE) and the other based on 1% and
99% of the mean volume fraction (VF). A third method for estimating $\theta$
was also considered, based on the relationship between the total fluctuating
kinetic energy decay rate $n$ and $\theta$ that is derived under the
assumption that the mixing layer growth is self-similar.

In general, estimates of $\theta$ using either definition for $h$ were found
to be in good agreement with one another, particularly for the $m=-3$
broadband perturbation that is the most representative of the initial
conditions used in the experiments of Sewell et al. (2021). The estimates of
$\theta$ based on $h$ for all three broadband cases were between 0.44 and
0.52, which is in very good agreement with the experimental estimates in
Sewell et al. (2021), who found $\theta=0.45\pm 0.08$ for their low-amplitude
cases and $\theta=0.51\pm 0.04$ for their high-amplitude cases prior to
reshock. When the TKE decay rate was used to estimate $\theta$ the results
were generally close to the estimates based on $h$, indicating that the mixing
layer growth is close to self-similar by the end of the simulation. Comparing
the ILES and DNS results also shows that there is only a small Reynolds number
dependence, which is consistent with previous observations in Groom & Thornber
(2019) that the integral quantities are mostly determined by the largest
scales of motion. When the mixing widths were decomposed into individual
bubble and spike heights $h_{b}$ and $h_{s}$, it was found that $h_{b}\sim
t^{\theta_{b}}$ and $h_{s}\sim t^{\theta_{s}}$ with $\theta_{b}\neq\theta_{s}$
at early time. However, it was shown that $\theta_{b}\approx\theta_{s}$ by the
end of each simulation by examining the ratio of $h_{s}/h_{b}$ and showing
this to be tending towards a constant at late time.

The particular regime being analysed here is different to the self-similar
growth regime analysed in Groom & Thornber (2020) as the current set of
broadband perturbations have a much smaller bandwidth and therefore saturate
quite early relative to the total simulation time. The present findings, which
are supported by the experiments, are that while the growth rate in the
saturated regime is less sensitive to the specific power spectrum of the
initial conditions, the effects of additional long wavelength modes are quite
persistent over the duration of a typical shock tube experiment and give rise
to growth rates much higher than for narrowband perturbations.

Comparing $\theta$ for the two definitions of $h$ in the narrowband case also
leads to some interesting observations. For the TKE-based mixing layer width
the value of $\theta$ that is measured is almost a factor of two higher than
the value that is measured for the VF-based width. This is due to spikes that
penetrate further into the lighter fluid and in some cases are ejected from
the main layer. These spikes have been observed in previous studies of similar
cases, such as Thornber & Zhou (2012); Youngs & Thornber (2020a), and are
quite energetic but contain very little heavy material. Therefore they affect
the TKE-based width much more than the VF-based width, which can be seen in
the greater relative difference between the two measures for the spike height
$h_{s}$ than the bubble height $h_{b}$. Presumably if such spikes are ejected
at early time in the broadband cases then they get overtaken by the linear
growth of the low wavelength modes; future work will investigate this in
further detail as it is potentially quite an important phenomenon for
applications where multiple interfaces are located in close proximity to one
another. Future work will also aim to further quantify the effects of finite
bandwidth on $\theta$ and other important integral quantities, see Soulard &
Griffond (2022) for an initial discussion in this direction.

Analysing the anisotropy of the fluctuating velocity field showed that the
mixing layer is persistently anisotropic in the direction of the shock wave in
all cases, in good agreement with previous experiments (prior to reshock) as
well as numerical studies. For the broadband ILES cases, the energy spectra in
both the normal and transverse directions showed two distinct scalings either
side of the highest wavenumber $k_{max}$ in the initial perturbation and which
were dependent on the specific initial condition. These scalings were also
different for the normal vs. transverse energy spectrum in each case. This was
also observed in the narrowband case but only for wavenumbers higher than
$k_{max}$. Finally, calculations of outer-scale Reynolds numbers and turbulent
length scales in the DNS cases showed that the outer-scale Reynolds numbers
are approximately constant throughout the simulations, as expected from the
estimates of $\theta\approx 0.5$, and that good agreement was obtained between
the turbulent length scales calculated using either the TKE-based or VF-based
width as the outer length scale.

Overall the results of this study show that, in general, care needs to be
taken when using measurements based on the velocity field to infer properties
of the concentration field such as the growth rate $\theta$. This is
particularly true when using thresholds rather than integral quantities to
represent the mixing layer width. At early times (i.e. prior to reshock in a
typical shock tube experiment) the mixing layer is not growing self-similarly,
which makes it difficult to determine the value for the growth rate exponent
$\theta$ as a single value may not even be appropriate. However, at the latest
time simulated here (just prior to reshock in the experiments of Jacobs et al.
(2013); Sewell et al. (2021)) the mixing layer is tending toward self-
similarity and good agreement was able to be obtained with the experimental
results across a wide range of quantities, providing additional insight on how
to correctly interpret such results and when it is valid to use a single
growth rate to describe the mixing layer.

[Acknowledgements]The authors would like to acknowledge David Youngs for
providing useful advice and insight on the formulation of the initial
conditions, as well as Jeffrey Jacobs for helpful discussions on the
computational setup and interpreting experimental results.

The authors would also like to acknowledge the computational resources at the
National Computational Infrastructure provided through the National
Computational Merit Allocation Scheme, as well as the Sydney Informatics Hub
and the University of Sydney’s high performance computing cluster Artemis,
which were employed for all cases presented here.

[Declaration of interests]The authors report no conflict of interest.

[Author ORCID]M. Groom, https://orcid.org/0000-0003-2473-7229; B. Thornber,
https://orcid.org/0000-0002-7665-089X

##  Appendix A Grid Convergence of Direct Numerical Simulations

Following the methodology presented in Olson & Greenough (2014) for
demonstrating grid convergence, two quantities are used that depend on
gradients of the velocity and concentration fields and which are therefore
sensitive to the smallest scales in the flow. These are the domain-integrated
enstrophy $\Omega$ and scalar dissipation rate $\chi$, given by

$\Omega(t)=\iiint\rho\omega_{i}\omega_{i}\>\mathrm{d}x\>\mathrm{d}y\>\mathrm{d}z$
(29)

and

$\chi(t)=\iiint D_{12}\frac{\partial Y_{1}}{\partial x_{i}}\frac{\partial
Y_{1}}{\partial x_{i}}\>\mathrm{d}x\>\mathrm{d}y\>\mathrm{d}z$ (30)

where $\omega_{i}$ is the vorticity in direction $i$ (summation over $i$ is
implied). Figures 21 and 22 demonstrate grid convergence in the domain-
integrated enstrophy and scalar dissipation rate for both DNS cases. Each case
is shown to be suitably converged for both of these integral quantities at the
finest grid resolution considered, even during the early-time period prior to
the shock exiting the domain.

(a) $Re_{0}=261$.

(b) $Re_{0}=526$.

Figure 21: Temporal evolution of domain integrated enstrophy for each grid
resolution employed in the DNS cases.

(a) $Re_{0}=261$.

(b) $Re_{0}=526$.

Figure 22: Temporal evolution of domain integrated scalar dissipation rate for
each grid resolution employed in the DNS cases.

##  Appendix B Integral Definitions of Bubble and Spike Heights

In Youngs & Thornber (2020a) novel definitions were given for the bubble and
spike heights $h_{b}$ and $h_{s}$ as weighted average distances from the
mixing layer centre, h_s^(p) & = [(p+1)(p+2)2∫-∞xc—x—p(1-⟨f1⟩) dx∫-∞xc(1-⟨f1⟩)
dx ]^1/p
h_b^(p) = [(p+1)(p+2)2∫xc∞—x—p⟨f1⟩ dx∫xc∞⟨f1⟩ dx ]^1/p. Figures 23 and 24 plot
the bubble and spike heights (with $p=3$), while figure 25 plots their ratio
$h_{s}/h_{b}$. The results are quite similar to the VF-based bubble and spike
heights shown in figures 12 to 14, albeit smoother and therefore more suitable
for estimating $\theta_{b}$ and $\theta_{s}$. While the main purpose of this
paper is to compare the quantities typically measured in experiments based on
thresholds of the TKE or VF profiles, it is recommended that future studies
focus on using integral definitions such as the ones given here.

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$ and $m=0$.

Figure 23: Temporal evolution of the bubble height $h_{b}$ based on the
integral definitions of Youngs & Thornber (2020a).

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$ and $m=0$.

Figure 24: Temporal evolution of the spike height $h_{s}$ based on the
integral definitions of Youngs & Thornber (2020a).

(a) $m=-1$.

(b) $m=-2$.

(c) $m=-3$ and $m=0$.

Figure 25: Temporal evolution of the ratio of spike to bubble heights based on
the integral definitions of Youngs & Thornber (2020a).

## References

  * Abgrall (1996) Abgrall, R. 1996 How to prevent oscillations in multicomponent flow calculations: A quasi conservative approach. J. Comput. Phys. 125 (125), 150–160.
  * Allaire et al. (2002) Allaire, G., Clerc, S. & Kokh, S. 2002 A Five-Equation Model for the Simulation of Interfaces between Compressible Fluids. J. Comput. Phys. 181, 577–616.
  * Andrews & Spalding (1990) Andrews, M.J. & Spalding, D.B. 1990 A simple experiment to investigate two-dimensional mixing by Rayleigh-Taylor instability. Phys. Fluids A 2, 922.
  * Arnett (2000) Arnett, D. 2000 The Role of Mixing in Astrophysics. Astrophys. J., Suppl. Ser. 127, 213.
  * Aslangil et al. (2020) Aslangil, Denis, Farley, Zachary, Lawrie, Andrew G. W. & Banerjee, Arindam 2020 Rayleigh–Taylor Instability With Varying Periods of Zero Acceleration. J. Fluids Eng. 142 (12), 121103.
  * Balakumar et al. (2008) Balakumar, B. J., Orlicz, G. C., Tomkins, C. D. & Prestridge, K. P. 2008 Simultaneous particle-image velocimetry-planar laser-induced fluorescence measurements of Richtmyer-Meshkov instability growth in a gas curtain with and without reshock. Phys. Fluids 20, 124103.
  * Balasubramanian et al. (2012) Balasubramanian, S., Orlicz, G. C., Prestridge, K. P. & Balakumar, B. J. 2012 Experimental study of initial condition dependence on Richtmyer-Meshkov instability in the presence of reshock. Phys. Fluids 24 (3).
  * Banerjee & Andrews (2009) Banerjee, Arindam & Andrews, Malcolm J. 2009 3D Simulations to investigate initial condition effects on the growth of Rayleigh–Taylor mixing. Int. J. Heat Mass Transf. 52 (17), 3906–3917.
  * Barenblatt et al. (1983) Barenblatt, G I, Looss, G & Joseph, D D 1983 Nonlinear Dynamics and Turbulence. Pitman Publishing.
  * Clark et al. (2016) Clark, D. S., Weber, C. R., Eder, D. C., Haan, S. W., Hammel, B. A., Hinkel, D. E., Jones, O. S., Kritcher, A. L., Marinak, M. M., Milovich, J. L., Patel, P. K., Robey, H. F., Salmonson, J. D., Sepke, S. M. & Edwards, M. J. 2016 Three-dimensional simulations of low foot and high foot implosion experiments on the National Ignition Facility. Phys. Plasmas 23, 56302.
  * Cook & Dimotakis (2001) Cook, A. W. & Dimotakis, P. E. 2001 Transition stages of Rayleigh-Taylor instability between Miscible Fluids. J. Fluid Mech. 443, 69–99.
  * Cook & Zhou (2002) Cook, A. W. & Zhou, Y. 2002 Energy Transfer in Rayleigh-Taylor Instability. Phys. Rev. E 66, 26312.
  * Dimotakis (2000) Dimotakis, P. E. 2000 The mixing transition in turbulent flows. J. Fluid Mech. 409, 69–98.
  * Duff et al. (1962) Duff, R. E., Harlow, F. H. & Hirt, C. W. 1962 Effects of diffusion on interface instability between gases. Physics of Fluids 5 (4), 417–425.
  * Elbaz & Shvarts (2018) Elbaz, Y. & Shvarts, D. 2018 Modal model mean field self-similar solutions to the asymptotic evolution of Rayleigh-Taylor and Richtmyer-Meshkov instabilities and its dependence on the initial conditions. Phys. Plasmas 25, 062126.
  * Groom & Thornber (2019) Groom, M. & Thornber, B. 2019 Direct numerical simulation of the multimode narrowband Richtmyer–Meshkov instability. Comput. Fluids 194, 104309.
  * Groom & Thornber (2020) Groom, M. & Thornber, B. 2020 The influence of initial perturbation power spectra on the growth of a turbulent mixing layer induced by Richtmyer–Meshkov instability. Physica D 407, 132463\.
  * Groom & Thornber (2021) Groom, M. & Thornber, B. 2021 Reynolds number dependence of turbulence induced by the richtmyer–meshkov instability using direct numerical simulations. J. Fluid Mech. 908, A31.
  * Hill et al. (2006) Hill, D. J., Pantano, C. & Pullin, D. I. 2006 Large-eddy simulation and multi-scale modelling of a Richtmyer-Meshkov instability with reshock. J. Fluid. Mech 557, 29–61.
  * Jacobs et al. (2013) Jacobs, J. W., Krivets, V. V., Tsiklashvili, V. & Likhachev, O. A. 2013 Experiments on the Richtmyer-Meshkov instability with an imposed, random initial perturbation. Shock Waves 23 (4), 407–413.
  * Kim & Kim (2005) Kim, K. H. & Kim, C. 2005 Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows Part II: Multi-dimensional limiting process. J. Comput. Phys. 208, 570–615.
  * Krivets et al. (2017) Krivets, V. V., Ferguson, K.J. & Jacobs, J.W. 2017 Turbulent mixing, induced by the Richtmyer-Meshkov instability. In AIP Conference Proceedings, , vol. 1793, p. 150003.
  * Lindl et al. (2014) Lindl, J., Landen, O., Edwards, J. & Moses, E. 2014 Review of the national ignition campaign 2009-2012. Phys. Plasmas 21, 020501.
  * Livescu (2020) Livescu, Daniel 2020 Turbulence with Large Thermal and Compositional Density Variations. Annu. Rev. Fluid Mech. 52 (1), 309–341.
  * Lombardini et al. (2012) Lombardini, M., Pullin, D. I. & Meiron, D. I. 2012 Transition to turbulence in shock-driven mixing: A Mach number study. J. Fluid Mech. 690, 203–226.
  * Massoni et al. (2002) Massoni, J., Saurel, R., Nkonga, B. & Abgrall, R. 2002 Proposition de méthodes et modèles eulériens pour les problèmes à interfaces entre fluides compressibles en présence de transfert de chaleur. Int. J. Heat Mass Transf. 45 (6), 1287–1307.
  * Meshkov (1969) Meshkov, E. E. 1969 Instability of the interface of two gases accelerated by a shock wave. Fluid Dyn. 43, 101–104.
  * Meyer & Blewett (1972) Meyer, K.A. & Blewett, P.J. 1972 Numerical investigation of the stability of a shock-accelerated interface between two fluids. Phys. Fluids 15, 753.
  * Mohaghar et al. (2017) Mohaghar, M., Carter, J., Musci, B., Reilly, D., McFarland, J. & Ranjan, D. 2017 Evaluation of turbulent mixing transition in a shock-driven variable-density flow. J. Fluid Mech. 831, 779–825.
  * Mohaghar et al. (2019) Mohaghar, M., Carter, J., Pathikonda, G. & Ranjan, D. 2019 The transition to turbulence in shock-driven mixing: Effects of Mach number and initial conditions. J. Fluid Mech. 14, 595–635.
  * Olson & Greenough (2014) Olson, B. J. & Greenough, J. 2014 Large eddy simulation requirements for the Richtmyer-Meshkov instability. Phys. Fluids 26, 044103.
  * Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.
  * Ramaprabhu et al. (2005) Ramaprabhu, P., Dimonte, Guy & Andrews, M. J. 2005 A numerical study of the influence of initial perturbations on the turbulent Rayleigh–Taylor instability. J. Fluid Mech. 536, 285–319.
  * Reese et al. (2018) Reese, D. T., Ames, A. M., Noble, C. D., Oakley, J. G., Rothamer, D. A. & Bonazza, R. 2018 Simultaneous direct measurements of concentration and velocity in the Richtmyer-Meshkov instability. J. Fluid Mech. 849, 541–575.
  * Richtmyer (1960) Richtmyer, R. D. 1960 Taylor instability in shock acceleration of compressible fluids. Comm. Pure Appl. Math. 13, 297–319.
  * Sewell et al. (2021) Sewell, E.G., Ferguson, K.J., Krivets, V.V & Jacobs, J.W. 2021 Time-resolved particle image velocimetry measurements of the turbulent Richtmyer–Meshkov instability. J. Fluid Mech. 917, A41.
  * Soulard & Griffond (2022) Soulard, O. & Griffond, J. 2022 Permanence of large eddies in Richtmyer–Meshkov turbulence for weak shocks and high Atwood numbers. Phys. Rev. Fluids 7, 014605.
  * Soulard et al. (2018) Soulard, O., Guillois, F., Griffond, J., Sabelnikov, V. & Simoëns, S. 2018 Permanence of large eddies in Richtmyer-Meshkov turbulence with a small Atwood number. Phys. Rev. Fluids 3, 104603.
  * Spiteri & Ruuth (2002) Spiteri, R. J. & Ruuth, S. J. 2002 A class of optimal high-order strong-stability preserving time discretization methods. SIAM J. Num. Anal. 40, 469–491.
  * Thornber (2016) Thornber, B. 2016 Impact of domain size and statistical errors in simulations of homogeneous decaying turbulence and the Richtmyer-Meshkov instability. Phys. Fluids 28, 45106.
  * Thornber et al. (2010) Thornber, B., Drikakis, D., Youngs, D. L. & Williams, R. J. R. 2010 The influence of initial conditions on turbulent mixing due to Richtmyer-Meshkov instability. J. Fluid Mech. 654, 99–139.
  * Thornber et al. (2011) Thornber, B, Drikakis, D, Youngs, D L & Williams, R J R 2011 Growth of a Richtmyer-Meshkov turbulent layer after reshock. Physics of Fluids 23 (9), 95107.
  * Thornber et al. (2017) Thornber, B., Griffond, J., Poujade, O., Attal, N., Varshochi, H., Bigdelou, P., Ramaprabhu, P., Olson, B., Greenough, J., Zhou, Y., Schilling, O., Garside, K. A., Williams, R. J. R., Batha, C. A., Kuchugov, P. A., Ladonkina, M. E., Tishkin, V. F., Zmitrenko, N. V., Rozanov, V. B. & Youngs, D. L. 2017 Late-time growth rate, mixing, and anisotropy in the multimode narrowband Richtmyer-Meshkov instability: The $\theta$-group collaboration. Phys. Fluids 29, 105107.
  * Thornber et al. (2018) Thornber, B., Groom, M. & Youngs, D. 2018 A five-equation model for the simulation of miscible and viscous compressible fluids. J. Comput. Phys. 372, 256–280.
  * Thornber et al. (2008) Thornber, B., Mosedale, A., Drikakis, D., Youngs, D. & Williams, R. 2008 An Improved Reconstruction Method for Compressible Flows with Low Mach Number Features. J. Comput. Phys. 227, 4873–4894.
  * Thornber & Zhou (2012) Thornber, B. & Zhou, Y. 2012 Energy transfer in the Richtmyer-Meshkov instability. Phys. Rev. E 86 (5), 1–11.
  * Toro et al. (1994) Toro, E. F., Spruce, M. & Speares, W. 1994 Restoration of the contact surface in the HLL-Riemann solver. Shock Waves 4, 25–34.
  * Tritschler et al. (2014a) Tritschler, V. K., Olson, Britton J., Lele, Sanjiva K., Hickel, S., Hu, X. Y. & Adams, Nikolaus A. 2014a On the Richtmyer-Meshkov instability evolving from a deterministic multimode planar interface. J. Fluid Mech. 755, 429–462.
  * Tritschler et al. (2014b) Tritschler, V. K., Zubel, M., Hickel, S. & Adams, N. A. 2014b Evolution of length scales and statistics of Richtmyer-Meshkov instability from direct numerical simulations. Phys. Rev. E 90, 1–11.
  * Vandenboomgaerde et al. (1998) Vandenboomgaerde, M., Mügler, C. & Gauthier, S. 1998 Impulsive model for the Richtmyer-Meshkov instability. Phys. Rev. E 58 (2), 1874–1882.
  * Vetter & Sturtevant (1995) Vetter, B. & Sturtevant, B. 1995 Experiments on the Richtmyer-Meshkov instability on a air/sf$_6$ interface. Shock Waves 4 (5), 247–252.
  * Walchli & Thornber (2017) Walchli, B. & Thornber, B. 2017 Reynolds number effects on the single-mode Richtmyer-Meshkov instability. Phys. Rev. E 95, 013104.
  * Weber et al. (2012) Weber, C., Haehn, N., Oakley, J., Rothamer, D. & Bonazza, R. 2012 Turbulent mixing measurements in the Richtmyer-Meshkov instability. Phys. Fluids 24, 074105.
  * Weber et al. (2013) Weber, Christopher R., Cook, Andrew W. & Bonazza, Riccardo 2013 Growth rate of a shocked mixing layer with known initial perturbations. J. Fluid Mech. 725, 372–401.
  * Weber et al. (2014) Weber, C. R., Haehn, N. S., Oakley, J. G., Rothamer, D. A. & Bonazza, R. 2014 An experimental investigation of the turbulent mixing transition in the Richtmyer–Meshkov instability. J. Fluid Mech. 748, 457–487.
  * Wong et al. (2022) Wong, M.L., Baltzer, J.R., Livescu, D. & Lele, S.K. 2022 Analysis of second moments and their budgets for Richtmyer-Meshkov instability and variable-density turbulence induced by reshock. Phys. Rev. Fluids 7, 044602.
  * Wong et al. (2019) Wong, Man Long, Livescu, Daniel & Lele, Sanjiva K. 2019 High-resolution Navier-Stokes simulations of Richtmyer-Meshkov instability with reshock. Phys. Rev. Fluids 4 (10), 104609.
  * Yang et al. (1993) Yang, J., Kubota, T. & Zukoski, E. 1993 Applications of shock-induced mixing to supersonic combustion. AIAA Journal 31 (5), 854–862.
  * Yang et al. (2014) Yang, Q., Chang, J. & Bao, W. 2014 Richtmyer-Meshkov Instability Induced Mixing Enhancement in the Scramjet Combustor with a Central Strut. Adv. Mech. Eng. 2014, 614189\.
  * Youngs (1994) Youngs, D. L. 1994 Numerical simulation of mixing by Rayleigh-Taylor and Richtmyer-Meshkov instabilities. Laser Part. Beams 12, 725–750.
  * Youngs (2004) Youngs, D. L. 2004 Effect of initial conditions on self-similar turbulent mixing. In Proceedings of the 9th International Workshop on the Physics of Compressible Turbulent Mixing.
  * Youngs & Thornber (2020a) Youngs, D. L. & Thornber, B. 2020a Buoyancy-Drag modelling of bubble and spike distances for single-shock Richtmyer-Meshkov mixing. Physica D: Nonlinear Phenomena 410, 132517\.
  * Youngs & Thornber (2020b) Youngs, D. L. & Thornber, B. 2020b Early time modifications to the buoyancy-drag model for Richtmyer-Meshkov mixing. J. Fluids Eng. 142, 121107–1.
  * Zhou (2001) Zhou, Y. 2001 A scaling analysis of turbulent flows driven by Rayleigh-Taylor and Richtmyer-Meshkov instabilities. Phys. Fluids 13, 538–543.
  * Zhou (2017a) Zhou, Y. 2017a Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. I. Phys. Reports 720-722, 1–136.
  * Zhou (2017b) Zhou, Y. 2017b Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. II. Phys. Reports 723-725, 1–160.
  * Zhou & Cabot (2019) Zhou, Y. & Cabot, W. H. 2019 Time-dependent study of anisotropy in Rayleigh-Taylor instability induced turbulent flows with a variety of density ratios. Phys. Fluids 31, 084106.
  * Zhou et al. (2016) Zhou, Y., Cabot, W. H & Thornber, B. 2016 Asymptotic behavior of the mixed mass in Rayleigh-Taylor and Richtmyer-Meshkov instability induced flows. Phys. Plasmas 23, 52712.
  * Zhou et al. (2003) Zhou, Y., Robey, H. F. & Buckingham, A. C. 2003 Onset of turbulence in accelerated high-Reynolds-number flow. Phys. Rev. E 67, 56305.
  * Zhou et al. (2021) Zhou, Y., Williams, R.J.R., Ramaprabhu, P., Groom, M., Thornber, B., Hillier, A., Mostert, W., Rollin, B., Balachandar, S., Powell, P.D., Mahalov, A. & Attal, N. 2021 Rayleigh–Taylor and Richtmyer–Meshkov instabilities: A journey through scales. Physica D 423, 132838\.