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. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013) and Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021). 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$ and $3.2$ mm in the experiments. As the effects of different initial impulses are of primary interest, it does not make 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\,{\rm m}\,{\rm s}^{-1}$, while in the experiments $A^+=0.71$ and $\Delta u=74\,{\rm m}\,{\rm s}^{-1}$. Therefore the non-dimensional time is given by

(3.1)\begin{equation} \tau=\frac{(t-t_0)A^+\Delta u}{\lambda_{min}}, \end{equation}

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\le \tau \le 118$ at the latest time prior to reshock in the experiments of Jacobs et al. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013) ($t-t_0=6.5$ ms) and $73.9\le \tau \le 81.5$ at the latest time prior to reshock in the experiments of Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021) ($t-t_0=4.5$ ms), assuming the same range of values for $\lambda _{min}$ of $2.9$–$3.2$ mm. Figure 7 shows a subset of the image sequence taken from a typical vertical shock tube experiment in Jacobs et al. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013) 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$ 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. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013).

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. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013) 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 Reference Cook and Dimotakis2001; Cook & Zhou Reference Cook and Zhou2002; Zhou & Cabot Reference Zhou and Cabot2019), 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)\begin{equation} h=x\left(\langle\, f_1\rangle=0.01\right)-x\left(\langle\, f_1\rangle=0.99\right). \end{equation}

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 choice of cutoff location as well as statistical fluctuations (Zhou & Cabot Reference Zhou and Cabot2019). For that purpose, an integral definition is typically used such as the integral width (Andrews & Spalding Reference Andrews and Spalding1990)

(3.3)\begin{equation} W=\int\langle\, f_1\rangle\langle\, f_2\rangle\,\mathrm{d} x. \end{equation}

If $f_1$ varies linearly with $x$ then $h=6W$ (Youngs Reference Youngs1994). See also the recent paper by Youngs & Thornber (Reference Youngs and Thornber2020a) 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. (Reference Sewell, Ferguson, Krivets and Jacobs2021), 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 (TKE) 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 TKE for each of the four initial conditions, defined as

(3.4)\begin{equation} \mathrm{TKE} = \tfrac{1}{2}\overline{u_i^{\prime}u_i^{\prime}}, \end{equation}

where $\psi ^{\prime }=\psi -\bar {\psi }$ indicates a fluctuating quantity and the ensemble average $\bar {\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.

Figure 8.

Plane-averaged profiles of 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); (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

In Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021) a definition for the mixing layer centre is given as the centroid of the mean TKE profile, i.e.

(3.5)\begin{equation} x_c = \frac{\displaystyle\int xf(x)\,\mathrm{d} x}{\displaystyle\int f(x)\,\mathrm{d} x}, \end{equation}

where $f(x)$ is the mean TKE 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

(3.6)\begin{equation} \int_{-\infty}^{x_c}\langle\, f_2\rangle\,\mathrm{d} x=\int_{x_c}^{\infty}\langle\, f_1\rangle\,\mathrm{d} x, \end{equation}

which has been used previously in both computational (Walchli & Thornber Reference Walchli and Thornber2017; Groom & Thornber Reference Groom and Thornber2021) and experimental (Krivets et al. Reference Krivets, Ferguson and Jacobs2017) 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.

Figure 9.

Temporal evolution of the mixing layer centre $x_c$, comparing the definition based on the centroid of the mean TKE profile with the definition based on the $x$-location of equal mixed volumes; (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

Figure 10 shows the temporal evolution of the mixing layer width, using both the visual width definition based on the mean volume fraction (VF) profile (referred to as the VF-based width) as well as the definition from Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021) based on the distance between the 5 % cutoff locations in the mean TKE profile (referred to as the TKE-based width). The mean VF $f_1$ at these 5 % cutoff locations is $\ge 0.997$ on the spike side ($x< x_c$) and $\le 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 VF 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. (Reference Sewell, Ferguson, Krivets and Jacobs2021), 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$.

Figure 10.

Temporal evolution of mixing layer width $h$ based on the distance between cutoff locations using either the mean TKE or mean VF 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$).

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

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 (Reference Elbaz and Shvarts2018), whereas the TKE-based width is growing at a much faster 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. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013) and Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021), 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.

An alternative method for estimating $\theta$ is also given in Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021), 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. (Reference Thornber, Drikakis, Youngs and Williams2010). 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 are 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. (Reference Thornber, Drikakis, Youngs and Williams2010, Reference Thornber2017). However, Groom & Thornber (Reference Groom and Thornber2020) 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.

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; (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

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 with the perturbations analysed in Groom & Thornber (Reference Groom and Thornber2020) 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. (Reference Sewell, Ferguson, Krivets and Jacobs2021) 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 were homogeneous (Barenblatt, Looss & Joseph Reference Barenblatt, Looss and Joseph1983)). 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 Reference Youngs2004; Thornber et al. Reference Thornber, Drikakis, Youngs and Williams2010; Soulard et al. Reference Soulard, Guillois, Griffond, Sabelnikov and Simoëns2018; Soulard & Griffond Reference Soulard and Griffond2022) 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 % VF cutoff (referred to as VF-based heights).

Figure 12.

Temporal evolution of the bubble height $h_b$ based on the distance between cutoff locations using either the mean TKE or mean VF 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 TKE or mean VF 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$).

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 (Reference Youngs and Thornber2020a) 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 (Reference Youngs and Thornber2020a), 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.

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; (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

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 Reference Groom and Thornber2021). This can also be observed for the integral definitions of $h_b$ and $h_s$ given in Appendix B.

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.

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. (Reference Sewell, Ferguson, Krivets and Jacobs2021). The first is a global measure of anisotropy, defined as

(3.7)\begin{equation} \mathrm{TKR}=\frac{2\times\mathrm{TKX}}{\mathrm{TKY}+\mathrm{TKZ}}, \end{equation}

where $\mathrm {TKX}=\tfrac {1}{2}\overline {u^{\prime }u^{\prime }}$, $\mathrm {TKY}=\tfrac {1}{2}\overline {v^{\prime }v^{\prime }}$ and $\mathrm {TKZ}=\tfrac {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

(3.8)\begin{equation} {\mathsf{b}}_{ij}=\frac{\overline{u_i^{\prime}u_j^{\prime}}}{\dfrac{1}{2}\overline{u_i^{\prime}u_i^{\prime}}}-\frac{1}{3}\delta_{ij}. \end{equation}

This tensor, specifically the $x$-direction principal component ${\mathsf{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 ${\mathsf{b}}_{11}$ as

(3.9) \begin{equation} \frac{2\overline{u^{\prime}u^{\prime}}}{\overline{v^{\prime}v^{\prime}}+\overline{w^{\prime}w^{\prime}}}=\frac{2{\mathsf{b}}_{11}+2/3}{2/3-{\mathsf{b}}_{11}}, \end{equation}

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 with the equivalent figure 13 in Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021) 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 (Reference Groom and Thornber2019) 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. (Reference Sewell, Ferguson, Krivets and Jacobs2021). 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. (Reference Sewell, Ferguson, Krivets and Jacobs2021) 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. Reference Thornber2017); 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.

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; (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

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. (Reference Sewell, Ferguson, Krivets and Jacobs2021) for ${\mathsf{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 ${\mathsf{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.

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); (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

Figure 17 shows the temporal evolution of ${\mathsf{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 ${\mathsf{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. (Reference Wong, Livescu and Lele2019) prior to reshock.

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; (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

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

(3.10)\begin{equation} E_i(\kappa)=\widehat{u^\prime_i}^{\dagger}\widehat{u^\prime_i}, \end{equation}

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 2-D 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 (2.9), 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 Reference Groom and Thornber2020). 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{\rm \pi} )=32$. Thornber et al. (Reference Thornber, Drikakis, Youngs and Williams2010) modified the analysis of Zhou (Reference Zhou2001) 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.

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; (a) $m=-1$, (b) $m=-2$, (c) $m=-3$, (d) $m=0$ ($R=2$).

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 Reference Thornber2016; Groom & Thornber Reference Groom and Thornber2019). 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. (Reference Thornber2017). 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 a $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. (Reference Sewell, Ferguson, Krivets and Jacobs2021) 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 set-up.

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. (Reference Jacobs, Krivets, Tsiklashvili and Likhachev2013) and Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021), 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

(3.11)\begin{equation} {Re}_h=\frac{\overline{\rho^+} h \dot{h}}{\bar{\mu}}. \end{equation}

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 ${Re}_h\gtrsim 1$–$2\times 10^4$ for fully developed turbulence (Dimotakis Reference Dimotakis2000). 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.

Figure 19.

Outer-scale Reynolds numbers vs time; (a) $m=-1$, (b) $m=-2$.

Dimotakis (Reference Dimotakis2000) showed that, for stationary flows, fully developed turbulence is obtained when $\lambda _L/\lambda _V\ge 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

(3.12)\begin{gather} \lambda_L = 5 {Re}_h^{{-}1/2}h, \end{gather}

(3.13)\begin{gather}\lambda_V = 50{Re}_h^{{-}3/4}h, \end{gather}

from which it can be shown that ${Re}_h\geq 10^4$ for fully developed turbulence. For a time-dependent flow, Zhou et al. (Reference Zhou, Robey and Buckingham2003) 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

(3.14)\begin{equation} \min(\lambda_L,\lambda_D)>\lambda_V. \end{equation}

Figure 20 shows the temporal variation of each length scale in (3.14), 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 ${Re}_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.

Figure 20.

The Liepmann–Taylor (circles), inner-viscous (squares) and diffusion length scales vs time for both definitions of the outer-scale Reynolds number; (a) $m=-1$, (b) $m=-2$.

Sewell et al. (Reference Sewell, Ferguson, Krivets and Jacobs2021) 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 ${Re}_h$, for $\theta \le 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 \le 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 \le 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.