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

The effects of initial conditions on the evolution of the Richtmyer-Meshkov instability (RMI) at early to intermediate times are analysed, using numerical simulations of an idealised version of recent shock tube experiments performed at the University of Arizona. The experimental results are bracketed by performing both implicit large-eddy simulations of the high-Reynolds-number limit as well as direct numerical simulations (DNS) at Reynolds numbers lower than those observed in the experiments. Various measures of the mixing layer width, based on both the plane-averaged turbulent kinetic energy and volume fraction profiles are used to explore the effects of initial conditions on $\theta$ and are compared with the experimental results. The decay rate 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. 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.


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(Yang et al. , 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 ℎ. At late time ℎ scales as ∼ 2 for RTI and ∼ for RMI where the exponent ⩽ 1 has been shown to depend on initial conditions (Youngs 2004;Thornber et al. 2010). Various approaches have been taken to define ℎ, 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 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 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 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 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 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(Weber et al. , 2014Reese et al. 2018;Mohaghar et al. 2017Mohaghar et al. , 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 that are higher than the 0.25-0.33 expected for narrowband initial conditions. For example Weber et al. (2012Weber et al. ( , 2014 measured in the range 0.43-0.58, while later experiments on the same facility by Reese et al. (2018) obtained = 0.34 ± 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 ( = 0.51 and = 0.54), while at later times the growth rates were smaller but also more different ( = 0.38 and = 0.29 respectively). Krivets et al. (2017) also found a wide range of for the integral width prior Figure 1: A schematic of the problem setup. The major ticks correspond to a grid spacing of Δ = 1.0 m. The interface is initially located at = 3.0 m and the shock is initially located at = 2.5 m in the light fluid and travels from light to heavy.
to reshock, ranging from = 0.18 to = 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 = 0.45 ± 0.08 and = 0.51 ± 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 . 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.

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 Focus on Fluids articles must not exceed this page length 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 SF 6 , 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 In (2.1), is the mass density, = [ , , ] is the mass-weighted velocity vector, is the pressure, is the volume fraction of species and = + is the total energy per unit mass, where = 1 2 · is the kinetic energy and the internal energy is given by the equation of state. Note that only (2.1 ) 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 where is the ratio of specific heats of the mixture. For the five-equation model this is given by which is an isobaric closure (individual species temperatures are retained in the mixture). The viscous stress tensor for a Newtonian fluid is where is the dynamic viscosity of the mixture. Note that in (2.4) the bulk viscosity is assumed to be zero according to Stokes' hypothesis. The heat flux = + , with the conductive heat flux given by Fourier's law where is the thermal conductivity of the mixture, and is the temperature. The thermal conductivity of species is calculated using kinetic theory as = 5 4 R + , , while the thermal conductivity of the mixture (as well as the mixture viscosity) is calculated using Wilke's rule. The enthalpy flux , arising from changes in internal energy due to mass diffusion, is given by where ℎ = , is the enthalpy of species and , the specific heat at constant pressure. The diffusion flux on the RHS of (2.1 ) invokes Fick's law of binary diffusion, written in terms of volume fraction. is the molecular weight of species , is the molecular weight of the mixture and the binary diffusion coefficient 12 is calculated by assuming both species have the same Lewis number (Le 1 = Le 2 = Le), such that 12 = Le¯( 2.7) with¯the specific heat at constant pressure for the mixture. Finally in (2.1 ), M = 1 − 2 1 1 + 2 2 and = / is the number density.

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 finitevolume 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 fifthorder 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 The additional terms in 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 fiveequation 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 shockinduced turbulent mixing problems (see Thornber et al. 2010Thornber et al. , 2011Groom & Thornber 2019).

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 × × = × × where = 2 m is used for all simulations. The extent of the domain in the -direction is either = 1.5 for the ILES cases or = 0.75 for the DNS cases. Periodic boundary conditions are used in the -and -directions, while in the -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 = 2.5 m and 0 = 3.0 m respectively and the initial pressure and temperature of both (unshocked) fluids is = 0.915 atm and = 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 for each species. A schematic of the initial condition is shown in Figure 1.
The shock Mach number is = 1.5, which is higher than the = 1.2 shock used in Jacobs et al. (2013); Krivets et al. (2017) and the = 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 postshock 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 SF 6 are 1 = 1.083 kg/m 3 and 2 = 5.465 kg/m 3 and the post-shock densities are + 1 = 2.469 kg/m 3 and + 2 = 15.66 kg/m 3 respectively. This gives a post-shock Atwood number of + = 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 and 1 across the interface are computed based on the surface perturbation described in (2.8) below. The evolution of the interface is solved in the post-shock frame of reference by applying a shift of Δ = −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. 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. 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 = 298 K, while the diffusivities are calculated under the assumption of Lewis number equal to unity (hence Pr = Sc ). 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 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 (2.9) Three broadband initial conditions are simulated, containing length scales in the range = /2 to = /32 and with a spectral exponent = −1, −2 and −3 respectively. The choice of bandwidth = / = 16 is based on estimates of the minimum initial wavelength performed in Jacobs et al. (2013) of = 2.9 to 3.2 mm, relative to a test section width of = 8.9 × 10 −2 m. When scaled to the dimensions of the experiment, the perturbations in this study all have a minimum wavelength of = 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 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 /2 and , specifically = 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). In Figure  2 it can be seen that the = −3 initial condition is the closest match to the experiments (with an estimated slope of = −2.99 over the same range of modes), with the other perturbations included to study the effects of varying . A fourth perturbation (not shown) is also considered; a narrowband perturbation with a constant power spectrum (i.e. = 0) and length scales in the range = /16 to = /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 = 1, which is at the limit of the linear regime. Note that in the experiments of Jacobs et al. (2013), 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 δ = /4. The volume fractions 1 and 2 = 1 − 1 are computed as being the amplitude perturbation satisfying the specified power spectrum and 0 the mean position of the interface. The amplitude perturbation ( , ) is given by 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 ; only the amplitudes are varied. For full details on the derivation of the surface perturbation see Thornber et al. (2010Thornber et al. ( , 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 = −1 and = −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.
For each perturbation, the weighted-average wavelength can be defined as¯= 2 /¯, wherē ( 2.12) Similarly, the initial growth rate of the perturbation variance is given by  where + 0 = (1 − Δ / ) 0 is the post-shock standard deviation, 0 is the initial standard deviation and is a correction factor to account for the diffuse interface (Duff et al. 1962;Youngs & Thornber 2020b). Here = ( − + + )/(2 + ) is an additional correction factor that is applied to the Richtmyer compression factor = (1 − Δ / ) to give the impulsive model of Vandenboomgaerde et al. (1998). For the present gas combination and configuration, = 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 = ∫ ⟨ 1 ⟩⟨ 2 ⟩ d is equal to 0.564 and therefore 0 = 0.564 0 . For the DNS cases, the initial Reynolds number is calculated in line with previous studies as (2.14) + = 9.065 kg/m 3 is the mean post-shock density. Table 2 gives the initial growth rate and weighted-average wavelength for each perturbation.

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× the bandwidth of initial modes compared to the narrowband perturbation analysed in Groom & Thornber (2019, 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 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 Re ℎ = ℎℎ/ where the width ℎ ∝ then the Reynolds number, and hence the grid resolution requirements, can either increase or decrease in time depending on the value of since (2.15) Therefore for < 1/2 the Reynolds number is decreasing and vice versa for > 1/2. Youngs (2004); Thornber et al. (2010) showed that the value of depends on both the bandwidth and spectral slope of the initial condition, which was recently demonstrated in Groom & Thornber (2020) using ILES for perturbations of the form given by (2.9) with = −1, −2 and −3. For the largest bandwidths simulated, these perturbations gave values of = 0.5, 0.63 and 0.75 respectively, which for the = −1 and −2 cases are quite close to the theoretical values of = 1/2 and = 2/3. What these results imply is that the Reynolds number of a broadband perturbation with ⩽ −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 ⩽ 1/3 (Elbaz & Shvarts 2018;Soulard et al. 2018).
For DNS of narrowband RMI the number of cells per 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 Re 9/4 and the total number of floating point operations required to perform a simulation to a given time scales as 3 (Pope 2000). For transitional RMI, empirically the scaling appears to be less severe (closer to Re 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 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.

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 Case Re 0 Simulation time (s) Domain size (m 3 ) Grid resolution initial conditions as well as two DNS; one for the = −1 initial condition and one for the = −2 initial condition. The viscosity used in these DNS is = 0.3228 Pa·s, which corresponds to initial Reynolds numbers of Re 0 = 261 and Re 0 = 526 for the = −1 and = −2 cases respectively. While this viscosity is much higher than would occur experimentally, it is equivalent to using a much smaller value of 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)   In the narrowband case the mixing layer has remained relatively uniform over the span of the domain, whereas in the broadband cases, particularly the = −2 and = −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 1 = 0.001 and 1 = 0.999 are plotted for both the = 0 narrowband case and the = −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. 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.

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 , 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

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.
( 3.2) 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 If 1 varies linearly with then ℎ = 6 (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 -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 where ′ = − indicates a fluctuating quantity and the ensemble average = ⟨ ⟩ is calculated as a plane average taken over the statistically homogeneous directions (in this case and ). The volume fraction profile ⟨ 1 ⟩⟨ 2 ⟩ is also shown on the right axis of each plot, as well as the (outermost) -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 ( < ) is further from the mixing layer centre than in the = −1 and = −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 ⟨ 1 ⟩⟨ 2 ⟩, 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 , is discussed in further detail §3.3.
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.
where ( ) 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 -location of equal mixed volumes, 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 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 = −3 and = 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. 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 1 at these 5% cutoff locations is ⩾ 0.997 on the spike side ( < ) and ⩽ 0.003 on the bubble side ( > ) 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 ℎ = ( − 0 ) , the growth rate exponent 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¯ 0 > 1 so that the flow is sufficiently developed. The estimated value of 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 .
Analysing the results in table 4, there is good agreement between the values of 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 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 = −3 ILES case. For the narrowband ILES case however, the VF-based (and integral) width is growing at close to the theoretical value of = 1/3 for self-similar decay proposed by Elbaz & Shvarts (2018), whereas the TKE-based width is growing at a much fast rate of = 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 TKEbased 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 = −1 case, the value of obtained from the integral width is Reynolds number independent, while for = −2 the value of 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 = −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 (a) = −1.
Figure 10: Temporal evolution of mixing layer width ℎ 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.
For the narrowband case the estimate of 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 = −1 and = −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 = 3 − 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 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 with those in Sewell et al. (2021) using both the TKE-based width and TKE decay rate, the = −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 = 0.45 and the TKE decay rate measurements gave = 0.68 (which would correspond to no decay of TKE if the layer was homogeneous (Barenblatt et al. 1983)). The equivalent results in the = −3 simulation were = 0.493 and = 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 = 0.51, indicating that the turbulence in the mixing layer is more developed and closer to self-similar prior to reshock. The = −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.

Bubble and Spike Heights
In order to help better explain the estimates for given in table 4, it is useful to decompose the TKE-based and VF-based widths into separate bubble and spike heights, ℎ and ℎ , defined as the distance from the mixing layer centre 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 -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 ℎ and ℎ 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).
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 ℎ and ℎ are greater than the corresponding VFbased 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 ℎ /ℎ vs. time and shows that ℎ /ℎ > 1 for all cases. The same trend was observed in Youngs & Thornber (2020a) for both = 0.5 and = 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 (a) = −1.
(d) = 0 ( = 2). Figure 12: Temporal evolution of the bubble height ℎ 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.
This analysis provides evidence that, prior to reshock, ℎ and ℎ 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 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 , 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 TKR = 2 × TKX TKY + TKZ (3.7) where TKX = 1 2 ′ ′ , TKY = 1 2 ′ ′ and TKZ = 1 2 ′ ′ , 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 This tensor, specifically the -direction principal component 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 -direction) can be written in terms of 11 as 2 ′ ′ ′ ′ + ′ ′ = 2 11 + 2/3 2/3 − 11 (3.9) 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 = −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 increases. For the = 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 -group quarter-scale case ; 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.
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 = −3 case. This mirrors the results shown in Sewell et al. (2021) for 11 , with quite good agreement observed between the = −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 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. Figure 17 shows the temporal evolution of 11 at the mixing layer centre, both for the definition of 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 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.

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 is the complex conjugate of this transform. As isotropy is expected in the transverse directions, the ( ) and ( ) spectra are averaged to give a single transverse spectrum ( ). 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 (2.9), a scaling of ( ) ∼ ( +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 = ( /2 ) = 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 ( ) ∼ ( −6)/4 . This scaling is observed for the transverse spectra at wavenumbers greater than , while for the normal spectra a scaling of ( ) ∼ ( −5)/4 is observed, the reason for which is currently unclear.
For wavenumbers less than the normal spectra scale as −3/2 in the = −2 and = −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 = −1 case also has a scaling that is less steep than −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 ( ) ∼ ( +2)/2 and a −3/2 scaling. For the transverse spectra in each of the broadband cases at wavenumbers less than a similar trend is observed, with each spectrum having a scaling that is shallower than −3/2 . The same argument of transition between an ( ) ∼ ( +2)/2 and a −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 = 47. This makes it difficult to compare and verify the current findings with their existing experimental setup.

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 ℎ from which to compute the outer-scale Reynolds number as (3.11) 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 ℎ 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 Re ℎ ≳ 1-2 × 10 4 for fully developed turbulence (Dimotakis 2000). For both the = −1 and = −2 perturbations the VF-based Reynolds number is approximately constant in time, consistent with the measured values of given in table 4. Dimotakis (2000) showed that for stationary flows, fully developed turbulence is obtained when / ⩾ 1 where = 5 is the Liepmann-Taylor length scale and = 50 is the inner-viscous length scale, with and the Taylor and Kolmogorov length scales respectively. These length scales may be related to the outer-scale Reynolds number by  from which it can be shown that Re ℎ ⩾ 10 4 for fully developed turbulence. For a timedependent flow, Zhou et al. (2003) showed that an additional length scale = 5( ) 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 min( , ) > .
(3.14)  Figure 20 shows the temporal variation of each length scale in (3.14), with and calculated from the outer-scale Reynolds number using both definitions for ℎ. In both cases there is good agreement between the length scales calculated from either definition of Re ℎ . 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. Sewell et al. (2021) also observed < at all times prior to reshock in their lowamplitude experiments. The authors note that, because of the different dependence of each length scale on Re ℎ , for ⩽ 0.5 the flow can never transition to turbulence as will grow faster than . Furthermore, the definition for implies that it will be 0 at time = 0, which would seem to imply that an RMI-induced flow with ⩽ 0.5 can never become turbulent. However, the virtual time origin is neglected in the original definition for ; if it is included then this allows for the possibility that < at early time. In that situation, transition to turbulence will occur provided the initial velocity jump is strong enough to produce > for some period of time. The turbulence will still be decaying over time if ⩽ 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.

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. 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 were analysed by curve-fitting the expected power law behaviour for the mixing layer width ℎ to two different definitions of ℎ; 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 was also considered, based on the relationship between the total fluctuating kinetic energy decay rate and that is derived under the assumption that the mixing layer growth is self-similar.
In general, estimates of using either definition for ℎ were found to be in good agreement with one another, particularly for the = −3 broadband perturbation that is the most representative of the initial conditions used in the experiments of Sewell et al. (2021). The estimates of based on ℎ 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 = 0.45 ± 0.08 for their low-amplitude cases and = 0.51 ± 0.04 for their high-amplitude cases prior to reshock. When the TKE decay rate was used to estimate the results were generally close to the estimates based on ℎ, 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 ℎ and ℎ , it was found that ℎ ∼ and ℎ ∼ with ≠ at early time. However, it was shown that ≈ by the end of each simulation by examining the ratio of ℎ /ℎ 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 for the two definitions of ℎ in the narrowband case also leads to some interesting observations. For the TKE-based mixing layer width the value of 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 ℎ than the bubble height ℎ . 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 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 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 . 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 ≈ 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 . 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 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.

Appendix B. Integral Definitions of Bubble and Spike Heights
In Youngs & Thornber (2020a) novel definitions were given for the bubble and spike heights ℎ and ℎ as weighted average distances from the mixing layer centre, Figures 23 and 24 plot the bubble and spike heights (with = 3), while figure 25 plots their ratio ℎ /ℎ . 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 and . 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.