4.1. Phases of mixing-layer growth

Figure 5 plots time histories of the mixing-layer width $\mathcal {W}$ as defined by (3.20) using a $1\,\%\text {--}99\,\%$ mass-fraction criterion. At ${\sim }11\ \textrm {ns}$, the main shock arrives at the $\mathbb {M}_H$–$\mathbb {M}_L$ interface, initiating perturbation growth via the RM instability mechanism. At ${\sim }31\ \textrm {ns}$, the reshock arrives at and compresses the developing mixing layer. At ${\sim }33\ \textrm {ns}$, rapid post-reshock growth begins. The growth rate decreases at later times, with the maximum width reached at ${\sim }45\ \textrm {ns}$. Finally, from ${\sim }45$ to $50\ \textrm {ns}$, a weak compression event occurs. It corresponds to the impact of the reflection (moving in the $-x$ direction) of the main shock off the reshock ablator $\mathbb {M}_R$, which was initially accelerated towards the $\mathbb {M}_H$–$\mathbb {M}_L$ interface. This late-time event is analysed further in § 4.3.

Figure 5.

Evolution of the mixing-layer width $\mathcal {W}$ from (3.20) in the baseline simulations (Sim), along with NIF experimental measurements (Exp). For the simulations, the legend states $\mathcal {N}_{yz,2}$ as defined in the caption of table 1. The experimental data are tabulated and further discussed in the supplementary material.

Figure 5 also includes six experimental measurements of the mixing-layer width from NIF shots. As explained in § 3.6.2, the coarsest-resolution baseline simulation was tuned to these data $\{\mathcal {W}_{{exp},i}\}$ by adjusting the perturbation-component standard deviations $S^{{\dagger} }$ and $S^{*}$. Once fixed, the same parameters were used in all simulations. Accordingly, $\mathcal {W}$ agrees well with $\{\mathcal {W}_{{exp},i}\}$ in the 180-zone case. The higher-resolution simulations yield somewhat higher values for $\mathcal {W}$, including some values outside the bounds of estimated experimental error, although the same qualitative behaviour is observed in all cases. No experimental data are available prior to 30 ns, so it is not known how accurately the simulations capture the early post-first-shock behaviour of $\mathcal {W}$.

The mixing-layer development can be visualized via contours of the instantaneous $\mathbb {M}_H$ mass fraction $Y_H$. Figures 6 and 7 provide views of these contours at selected pre-reshock and post-reshock times, respectively, in the finest-resolution baseline simulation. For corresponding animations of the $Y_H$ contours and shock positions over time, see movies 1 and 2 in the supplementary material. The first pre-reshock figure 6(a) shows the $\mathbb {M}_H$–$\mathbb {M}_L$ interface at 10 ns, before the main shock arrives. The 2-D character of the initial perturbation is evident: in the centre-plane cross-section view of $Y_H$, horizontal bands correspond to the principal component $\delta _x^{{\dagger} }$, while smaller-scale irregularities correspond to the lower-amplitude noise component $\delta _x^{*}$. At the second pre-reshock time, 20 ns (figure 6b), the flow exhibits the mushroom-cap bubble and spike structures typically associated with the RM and RT instabilities. Principal-perturbation growth appears to be dominant. By the third pre-reshock time, 30 ns (figure 6c), vortices of various sizes and orientations have formed. The first post-reshock figure 7(a) shows the mixing layer just after reshock compression, with significant interpenetration of the $\mathbb {M}_H$ and $\mathbb {M}_L$ fluids. The second post-reshock time, 45 ns (figure 7b), corresponds roughly to the point of maximum mixing-layer width. The flow exhibits a spectrum of length scales and many zones that are well mixed (with $Y_H \approx Y_L \approx 0.5$). At the third post-reshock time at the end of the simulation, 50 ns (figure 7c), the mixing layer is slightly compressed, relative to its state at 45 ns. Again, many zones are well mixed, 3-D structures abound and the flow appears chaotic and turbulent.

Figure 6.

Contours of the instantaneous $\mathbb {M}_H$ mass fraction $Y_H$ in the $\mathbb {M}_H$–$\mathbb {M}_L$ mixing layer of the finest-resolution baseline simulation at three times before reshock: (a) 10 ns, (b) 20 ns and (c) 30 ns. The left frame of each figure depicts the 3-D field $Y_H( x, y, z)$ near the mixing-layer centre-plane $x=x_c$. Zones with $Y_H < 0.05$ are not shown. The orientation is rotated from the orientation of figure 1. The right frame of each figure depicts the 2-D cross-section $Y_H( x_c, y, z)$. Equation (4.3) defines the coordinate $x_c$. Movies 1 and 2 in the supplementary material accompany this figure.

Figure 7.

Contours of the instantaneous $\mathbb {M}_H$ mass fraction $Y_H$ in the $\mathbb {M}_H$–$\mathbb {M}_L$ mixing layer of the finest-resolution baseline simulation at three times after reshock: (a) 35 ns, (b) 45 ns and (c) 50 ns. The same conventions used in figure 6 apply here. Movies 1 and 2 in the supplementary material accompany this figure.

The quantity $\mathcal {W}$ is an important metric of the mixing-layer size, because it corresponds most naturally to the experimental data $\{\mathcal {W}_{{exp},i}\}$: calculated bubble and spike fronts based on mass-fraction thresholds correspond to observable edges in X-ray radiographs. The edges appear where X-ray transmission through the plasma changes, e.g. due to density and species variations (Nagel et al. Reference Nagel2017; Huntington et al. Reference Huntington, Raman, Nagel, MacLaren, Baumann, Bender, Prisbrey, Simmons, Wang and Zhou2020). However, $\mathcal {W}$ can be sensitive to the structure of individual jets and other flow features at the edges of the mixing layer (Zhou & Cabot Reference Zhou and Cabot2019). Alternate metrics based on volume integrals may suffer less from these sensitivities. One such metric, advocated by Zhou, Cabot & Thornber (Reference Zhou, Cabot and Thornber2016) in the context of ICF research, is the mixed mass

(4.1)\begin{equation} \mathcal{M}( t)= \int_{x_b}^{x_s} \int_{0}^{L} \int_{0}^{L} 4 \rho Y_H Y_L\, {\textrm{d}y}\, \textrm{d}z\, {\textrm{d}\kern0.7pt x}. \end{equation}

For a perfectly mixed simulated layer consisting only of zones with $Y_H = Y_L = 0.5$, observe that $\mathcal {M}$ is equal to the total mass of the layer. Figure 8(a) plots $\mathcal {M}$ versus time. The mixed mass is monotonically increasing and is notably not reduced by the reshock compression, illustrating its potential usefulness as a progress variable for mixing-layer development. Growth of $\mathcal {M}$ is much more rapid after reshock than after first shock, and a discernible increase in the magnitude of $d \mathcal {M}/\textrm {d}t$ occurs as soon as the reshock enters the mixing layer, i.e. before the time of maximum compression when the mixing-layer width $\mathcal {W}$ is locally minimal. As expected, $\mathcal {M}$ demonstrates slightly better grid-convergence behaviour than $\mathcal {W}$.

Figure 8.

Evolution of (a) the mixed mass $\mathcal {M}$ from (4.1) and (b) the normalized mixed mass $\varPsi$ from (4.2). All results are from the baseline simulations. The legends state $\mathcal {N}_{yz,2}$ as defined in the caption of table 1.

Figure 8(b) plots a related metric, the normalized mixed mass (Zhou et al. Reference Zhou, Cabot and Thornber2016)

(4.2)\begin{equation} \varPsi( t)= \frac{\displaystyle\int_{x_b}^{x_s} \overline{ \rho Y_H Y_L }\,{\textrm{d}\kern0.7pt x} } { \displaystyle\int_{x_b}^{x_s} \bar{\rho} (\overline{Y_H}) (\overline{Y_L}) \,{\textrm{d}\kern0.7pt x}}. \end{equation}

This quantity is related to the molecular mixing fraction introduced by Youngs (Reference Youngs1991, Reference Youngs1994) and analysed in many computational studies of shock-induced mixing (e.g. Latini, Schilling & Don Reference Latini, Schilling and Don2007; Thornber et al. Reference Thornber, Drikakis, Youngs and Williams2011; Lombardini et al. Reference Lombardini, Pullin and Meiron2012; Tritschler et al. Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014). It expresses a ratio of sub-zonal mixing to larger-scale entrainment; it is equal to unity for a perfectly mixed simulated layer, and it is equal to zero for a simulated layer consisting only of zones with either $Y_H= 1$ or $Y_L= 1$ in any distribution. (Note that quantities like $\varPsi$ have been described as ratios of atomic mixing to chunk mixing. However, those terms are misleading if the fluid-dynamical action of atomic-scale processes like viscous dissipation is not resolved on the computational mesh, as is the case in many investigations of complex mixing.) Before arrival of the first shock, $\varPsi$ increases slightly due to mass diffusion, which is not constrained by the WISI approximation described in § 3.4. The first shock initiates perturbation growth, stretching coherent packets of the $\mathbb {M}_H$ and $\mathbb {M}_L$ plasmas across the interfacial region and sharply decreasing $\varPsi$. Then, as smaller-scale structures develop and physical and numerical dissipation ensues, $\varPsi$ rises towards values between approximately 0.80 and 0.85, with an interruption due to reshock. Those values are roughly in agreement with the asymptotic behaviour of simulations of both singly shocked and reshocked non-HED mixing layers (Thornber et al. Reference Thornber, Drikakis, Youngs and Williams2011; Lombardini et al. Reference Lombardini, Pullin and Meiron2012; Tritschler et al. Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014). The figure indicates that the relative amount of sub-zonal mixed fluid in the earliest phases of post-first-shock perturbation growth is particularly sensitive to the mesh.

4.2. Mass-fraction profiles and definition of the mixing-layer centre-plane

Figure 9 shows instantaneous profiles of the Reynolds-averaged $\mathbb {M}_L$ mass fraction $\overline {Y}_L(x)$ at a late pre-reshock time in figure 9(a) and at approximately the time of maximum mixing-layer width in figure 9(b). As $\overline {Y}_L$ increases, the mass fraction $\overline {Y}_H$ of the heavy material $\mathbb {M}_H$ decreases and the density generally decreases. At each of the two times, root-mean-square (r.m.s.) simulation-to-simulation deviations between the profiles are ${<}0.04$. At 30 ns, $\overline {Y}_L$ is not monotonically increasing with $x$, a consequence of mostly single-species fluid entrained in large bubbles and spikes, e.g. as seen in figure 6(c). At the late post-reshock time, the mixing layer exhibits an inner core with approximately linear variation in the averaged $\mathbb {M}_L$ mass fraction, for all three meshes.

Figure 9.

Instantaneous profiles of the Reynolds-averaged mass fraction $\overline {Y}_{L}$ of the light material $\mathbb {M}_L$ versus axial distance at two times in the baseline simulations. Symbols indicate data extracted directly from the simulations (Sim). The data are under-sampled for clarity in the figure, i.e. for each case, there are more data available than there are symbols drawn. Lines indicate fits (Fit) to the simulation data using the analytic form (4.3). The abscissa limits are chosen such that the centre of each plot is $x=x_c$, as defined in the text, and the abscissa range is $250\ \mathrm {\mu }\textrm {m}$. The legends state $\mathcal {N}_{yz,2}$ as defined in the caption of table 1.

Recall from § 3.5 that $\overline {Y}_L$ underwrites the definitions of the bubble front $x_b$ and spike front $x_s$. Many of the analyses described here also require definition of a mixing-layer centre-plane $x_c$. Due to the non-monotonicity in some $\overline {Y}_L$ profiles, a definition based simply on a mass-fraction threshold (e.g. $\overline {Y}_L = 0.5$) is not robust and may be undesirably sensitive to local bubbles and spikes. Instead, we define $x_c$ by fitting the following error-function-based form to the $\overline {Y}_L$ data (compare (3.25)),

(4.3)\begin{equation} \overline{Y}_{L,{fit}}( x)= \frac{1}{2} \left[ 1+ \textrm{erf} \left( \frac{x - x_c}{ x_w} \right) \right], \end{equation}

where $x_w$ is the mixing-layer fitting width. The parameters $x_c$ and $x_w$ are determined via nonlinear least-squares fitting methods (Jones Reference Jones2001), and this definition for $x_c$ is used throughout the present study. For some analyses, we may replace $x_c$ with the $x$ coordinate of the nearest plane of zone centres or nodes on the finest AMR level. Figure 9 includes plots of the fitted curves at 30 and 45 ns. At each of the two times, r.m.s. simulation-to-simulation deviations between the curves are ${<}0.008$; the $\overline {Y}_{L,{fit}}$ curves are better converged than the $\overline {Y}_{L}$ profiles. Over all post-first-shock times, r.m.s. simulation-to-simulation deviations in $x_c$ are ${<}1\ \mathrm {\mu }\textrm {m}$. The post-reshock $\overline {Y}_L$ profiles in figure 9(b) are better represented by the analytic function $\overline {Y}_{L,{fit}}$ than are the pre-reshock profiles in figure 9(a). However, deviations from the fit near the edges of the mixing layer persist to late times. At the lower-density spike front, the simulations predict a sharper rise to $\overline {Y}_L = 1$ (as $x$ increases) than the analytic form. Conversely, at the higher-density bubble front, the simulations predict a shallower fall to $\overline {Y}_L = 0$ (as $x$ decreases) than the analytic form. Thus, relative to the simulation data, the error function fit mostly underestimates the light-material mass fraction at both mixing-layer edges.

Although not investigated in detail here, other definitions for the mixing-layer centre-plane are plausible. Consider three alternate definitions for $x_c$ at a given time: the $x$ coordinate of the midpoint of a line segment fit to the $\overline {Y}_L$ data from $x_b$ to $x_s$; the $x$ value for which the accumulated mixed volume $\int _{x_b}^{x} \iint 4 Y_H Y_L\, {\textrm {d}y}\,\textrm {d}z\,\textrm {d}s$ is half of its maximum; and the $x$ value for which the accumulated mixed mass $\int _{x_b}^{x} \iint 4 \rho Y_H Y_L\, {\textrm {d}y}\,\textrm {d}z\,\textrm {d}s$ is half of its maximum. Each of these three methods returns values of $x_c$ usually smaller than those derived from (4.3), with method-to-method r.m.s. deviations – calculated over all post-first-shock times and as percentages of the instantaneous mixing-layer width – less than $1.5\,\%$, $9.6\,\%$ and $14\,\%$, respectively. Note that the choice of definition for $x_c$ has no effect on integrated quantities computed from (3.21) or (3.22).

4.3. Thermodynamic properties and bulk velocity at the centre-plane

Figure 10 presents time histories of the Reynolds averages of various thermodynamic variables and of the axial component of velocity at the mixing-layer centre-plane in the finest-resolution baseline simulation. The first shock drives the centre of the mixing layer into the HED regime. From $\sim$11 to 31 ns, the averaged total pressure is $\sim$2–4 Mbar and the averaged electron temperature holds relatively steady at $\sim$15 eV. The first shock also imparts a large positive axial velocity to the mixing layer. From $\sim$11 to 31 ns, the averaged axial velocity decreases, i.e. the interface decelerates as it moves downstream towards unmixed light material $\mathbb {M}_L$. Hence, we expect $\mathbb {M}_H$–$\mathbb {M}_L$ perturbation growth due to both the RM instability, activated by the first-shock impact, and the RT instability, activated by the heavy–light interface deceleration. The presence of both mechanisms may explain the slightly super-linear growth of the mixing-layer width in figure 5 from $\sim$20 to 31 ns. The results imply that the flows realized in the NIF experiments cannot be accurately characterized as either ‘pure RM’ or ‘pure RT’.

Figure 10.

Evolution of various Reynolds-averaged flow variables at the mixing-layer centre-plane $x_c$, as defined by (4.3), in the finest-resolution baseline simulation: $p$ is the total pressure, $u_1$ is the axial component of velocity, $\rho$ is the density and $T_e$ is the electron temperature. Overbars are omitted to simplify the notation.

Reshock causes multiple-factor increases in the averaged total pressure, density and electron temperature at the centre-plane. The averaged axial velocity falls to approximately zero and remains small for the duration of the simulation. Hence, the reshock effectively halts the bulk forward motion of the mixing layer. The averaged total pressure, density and electron temperature all decay from their peak values as the doubly shocked plasma decompresses. It is notable that $p$ decreases by a factor of $\sim$2.8 from 33 to 45 ns; the post-reshock decompression observed here is generally much more intense than in the non-HED studies discussed in § 1. Near the end of the simulation, all three thermodynamic variables abruptly rise again, although by smaller amounts than immediately after reshock. This late-time event, previously mentioned in § 4.1, corresponds to the impact of the reflection of the main shock off the reshock ablator $\mathbb {M}_R$. Indeed, a high-density region of $\mathbb {M}_R$ plasma is created by the reshock at early times, moves in the $-x$ direction towards the $\mathbb {M}_H$–$\mathbb {M}_L$ interface, and serves as a reflection surface for the main shock (after the main shock has transmitted into the $\mathbb {M}_L$ region).

The simulations predict a small precursor in the radiation temperature, extending $\lesssim 20\ \mathrm {\mu }\textrm {m}$ ahead of the reshock in $\mathbb {M}_L$ before its impact on the $\mathbb {M}_H$–$\mathbb {M}_L$ interface. However, the precursor has a negligible effect on the ion temperature, electron temperature, and total pressure. Ahead of the main shock in $\mathbb {M}_H$ before its impact on the $\mathbb {M}_H$–$\mathbb {M}_L$ interface, temperature and pressure precursors are negligible. Indeed, the NIF experiments were intentionally designed to keep energies low enough such that the shocks could be reasonably interpreted as simple discontinuities in the thermodynamic quantities. (Our computational framework does predict strong radiative precursors at higher energies than those in the present study, although accurate assessment of such phenomena would likely require more-sophisticated treatments of radiation transport and opacities.)

Additional insight into the thermodynamics of the flow comes from Reynolds-averaged densities and electron temperatures at the mixing-layer edges. Figures 11(a) and 11(b) plot those quantities, respectively, at both the bubble and spike fronts. The time variations of the bubble-front density $\rho _b$ and the spike-front density $\rho _s$ are qualitatively similar to the time variation of the centre-plane density. Since the reshock and the reflected main shock strike the spike front before the bubble front, the corresponding rises in $\rho _s$ occur before those in $\rho _b$. Figure 11(b) indicates that there is a significant temperature gradient across the mixing layer. After first shock, the lower-density spike front remains at a much higher electron temperature than the higher-density bubble front. Throughout the shocked and reshocked mixing layers, we observe that packets of foam $\mathbb {M}_L$ tend to be much hotter than nearby packets of main ablator $\mathbb {M}_H$, with important fluid-dynamical consequences that will be analysed in § 5.

Figure 11.

Evolution of (a) the Reynolds-averaged density $\rho$ and (b) the Reynolds-averaged electron temperature $T_e$ at the mixing-layer bubble front $x_b$ and spike front $x_s$, as defined in § 3.5, in the finest-resolution baseline simulation. The subscripts $b$ and $s$ denote values at the bubble and spike fronts, respectively. Overbars are omitted to simplify the notation.

Corresponding time evolutions of ion temperature and radiation temperature are nearly identical to those of electron temperature in both figures 10 and 11(b); considering the plotted quantities after first shock, r.m.s. $T_n$–$T_e$ deviations are $<$0.02 eV and r.m.s. $T_r$–$T_e$ deviations are $<$0.07 eV. Indeed, throughout the mixing layer, $T_n$, $T_e$ and $T_r$ remain approximately equal, suggesting that both ion–electron energy coupling via (3.14) and electron–radiation energy coupling via the second term of the right-hand side of (3.6) are very fast compared to other dynamics in the mixing plasmas. Differences between the three temperatures are more significant near the $-x$ and $+x$ domain boundaries in the pure $\mathbb {M}_H$ and $\mathbb {M}_R$ materials, both during initial formation of the main shock and reshock at early times and after the materials have expanded to very low densities at late times.

While this section only presented results for the finest-resolution baseline simulation, the qualitative trends observed here apply to all three baseline simulations. Post-first-shock r.m.s. simulation-to-simulation deviations in the Reynolds-averaged centre-plane, bubble-front and spike-front quantities plotted in figures 10 and 11 are less than 0.2 Mbar for $p$, $0.2\ \mathrm {\mu }\textrm {m}\,\textrm {ns}^{-1}$ for $u_1$, $0.07\ \textrm {g}\,\textrm {cm}^{-3}$ for $\rho$ and 1 eV for $T_e$.

4.4. Evolution of turbulent kinetic energy

This section considers the kinetic energy in the fluctuating fluid motion in the developing mixing layers. Using the formalisms of § 3.5, define the local turbulent kinetic energy

(4.4)\begin{equation} \mathcal{I}= \tfrac{1}{2} u_i'' u_i'' \end{equation}

and the mean density-weighted turbulent kinetic energy

(4.5)\begin{equation} \mathcal{K}= \tfrac{1}{2} \overline{ \rho u_i'' u_i''} = \tfrac{1}{2} \bar{\rho} \widetilde{ u_i'' u_i''} = \bar{\rho} \widetilde{ \mathcal{I}}, \end{equation}

which are here abbreviated LTKE and MDTKE, respectively. Further discussion of $\mathcal {I}$, $\mathcal {K}$ and related quantities can be found in Sagaut & Cambon (Reference Sagaut and Cambon2008), Chassaing et al. (Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2010), Gatski & Bonnet (Reference Gatski and Bonnet2013) and the supplementary material. In the present study, $\mathcal {I}$ is a function of $x$, $y$, $z$ and $t$, while $\mathcal {K}$ is a function of $x$ and $t$ only.

Figure 12(a) plots $\lfloor {\mathcal {K}} \rfloor$, the mixing-layer integral of MDTKE, versus time. It rises slowly after first shock, stabilizes, increases sharply by over one order of magnitude after reshock and finally decays as the mixing layer grows and decompresses. To quantify the reshock-induced increase and the post-reshock decay, table 3 reports mixing-layer integral values at the representative times 30 ns (shortly before the reshock arrival), 35 ns (shortly after the reshock has moved through the mixing-layer and after the moment of peak compression) and 45 ns (when the mixing-layer width is approximately largest). The simulations predict that $\lfloor {\mathcal {K}} \rfloor$ increases by a factor of $\sim$30–40 from 30 to 35 ns. The magnitude of this increase, along with much of the history of $\lfloor {\mathcal {K}} \rfloor$, is only weakly sensitive to the mesh. Grid sensitivity is most noticeable just after first shock and at very late times. The large and rapid reshock-induced increase in MDTKE seen here is consistent with findings from analogous non-HED computational work (Hill et al. Reference Hill, Pantano and Pullin2006; Schilling & Latini Reference Schilling and Latini2010; Tritschler et al. Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014).

Figure 12.

In (a), evolution of the mixing-layer integral of MDTKE $\mathcal {K}$ from (4.5). In (b), turbulent energy spectra $\mathcal {R}$ from (E11) at the mixing-layer centre-plane at 30, 35 and 45 ns. The spectra are compensated via division by $f^{-5/3}$, where $f$ is the spectroscopic wavevector magnitude. Note the differences in the ordinate limits of the plots of (b). All results are from the baseline simulations. Equation (4.3) defines the mixing-layer centre-plane, and (3.22) defines the mixing-layer integral. The legends state $\mathcal {N}_{yz,2}$ as defined in the caption of table 1.

Table 3.

Selected values of the mixing-layer integral of MDTKE ($\textrm {g}\ \textrm {cm}^{-3}\, \mathrm {\mu }\textrm {m}^{2}\,\textrm {ns}^{-2}$) at three different times (ns). The ratio of the early post-reshock value to the late pre-reshock value and the ratio of the early post-reshock value to the late post-reshock value are given. Compare with figure 12(a).

Compensated turbulent energy spectra – derived from the Favre fluctuation of velocity per appendix E – are shown in figure 12(b) at the mixing-layer centre-plane at the same three times discussed above. On each plot, horizontal lines correspond to the canonical $f^{-5/3}$ behaviour of fully developed turbulence (Pope Reference Pope2000; Davidson Reference Davidson2015). At the two post-reshock times, the turbulent energy spectra $\mathcal {R}$ are larger at high wavevector magnitudes – relative to corresponding values at low wavevector magnitudes – than at the pre-reshock time, indicating that energy was transferred to smaller scales as the mixing evolved. However, even at the time when the integrated MDTKE is near its maximum, there are not broad intervals exhibiting $f^{-5/3}$ behaviour. At high wavevector magnitudes, the spectra are clearly not grid-converged, suggesting that the simulations are not close to resolving the scales associated with the physical viscosity $\mu$. It is not known whether more-distinct inertial ranges would appear in finer-resolution versions of the present simulations.

Appendix F expands on the discussion in this section. Various terms in an evolution equation for $\mathcal{K}$ are analysed, providing further insight into how energy exchanges between the mean and fluctuating flows.

4.5. Anisotropy, length scales and Reynolds numbers

The previous section considered only the total MDTKE, i.e. the sum of three terms corresponding to fluctuating motion in each coordinate direction. The flow in the present study is neither homogeneous nor isotropic. Anisotropy arises not only from the distinction between the axial ($x$) and spanwise ($y$, $z$) directions (which are, respectively, parallel and normal to the shock velocity vectors), but also from the experimentally based asymmetries in the interface initial perturbation described in § 3.6.2. This section analyses how the directional components of the MDTKE and related quantities evolve. Throughout this section, we refer to the coordinate directions interchangeably by the indices $i= 1, 2 \text{ or }3$ or by the letters $x,y \text { or }z$.

Define the $i$th component of MDTKE by

(4.6)\begin{equation} \mathcal{K}_{i} = \tfrac{1}{2} \overline{ \rho u_{(i)}'' u_{(i)}''}, \end{equation}

with no summation on $i$ implied. Note that $\mathcal {K}= \mathcal {K}_{1}+ \mathcal {K}_{2}+ \mathcal {K}_{3}$. A metric of the relative contribution of axial-direction fluctuations to $\mathcal {K}$, integrated across the mixing layer, is

(4.7)\begin{equation} \mathcal{Y}_{1}= \frac{ \lfloor { \mathcal{K}_{1}} \rfloor}{ \lfloor { \mathcal{K}} \rfloor} - \frac{ 1}{3}, \end{equation}

which equals 2/3 when the only contributions to $\mathcal {K}$ are the $u_1''$ fluctuations and equals zero in the limit of perfect isotropy. Figure 13(a) plots $\mathcal {Y}_{1}$ versus time, showing that the mixing layer is strongly anisotropic before reshock, with axial fluctuations dominating over spanwise fluctuations. This imbalance lessens after reshock, with $\mathcal {Y}$ rapidly decreasing towards its isotropic value at later times. Figure 13(b) complements figure 13(a), plotting the ratio of the mixing-layer integrals of the two spanwise components of $\mathcal {K}$. When $\lfloor {\mathcal {K}_{2}} \rfloor /\lfloor {\mathcal {K}_{3}} \rfloor \gg 1$, the $u_2''$ fluctuations are larger in magnitude than the $u_3''$ fluctuations, on average; conversely, in a perfectly isotropic flow, $\lfloor {\mathcal {K}_{2}} \rfloor /\lfloor {\mathcal {K}_{3}} \rfloor = 1$. The curves in figure 13(b) indicate that, just after first shock, $y$-direction contributions to the MDTKE drastically outweigh $z$-direction contributions – a trend expected given that the standard deviation of the principal perturbation is much larger than that of the noise perturbation per § 3.6.2. The energy distribution among the two spanwise components changes rapidly, with the ratio $\lfloor {\mathcal {K}_{2}} \rfloor /\lfloor {\mathcal {K}_{3}} \rfloor$ dropping almost monotonically. The ratio appears to reach an asymptotic value greater than unity at late times, i.e. the $u_2''$ fluctuations persistently account for a larger fraction of $\lfloor {\mathcal {K}} \rfloor$ than the $u_3''$ fluctuations, even after the mixing layer appears well developed. However, caution is necessary when interpreting figure 13(b): the spanwise-component ratio is clearly not grid converged. The finer computational meshes are better capable of resolving the small scales of the noise perturbation (e.g. scales comparable to $\lambda _{{min}}^{*}= 2\ \mathrm {\mu }\textrm {m}$) than the coarsest mesh. Since only the noise perturbation – not the principal perturbation – seeds initial growth of $\mathcal {K}_{3}$, the finer-resolution simulations generally predict larger $u_3''$-fluctuation magnitudes than the coarsest-resolution simulation.

Figure 13.

Anisotropy metrics, Taylor microscales and Taylor Reynolds numbers in the baseline simulations: (a) plots $\mathcal {Y}_{1}$ from (4.7); (b) plots the ratio $\lfloor {\mathcal {K}_2} \rfloor /\lfloor {\mathcal {K}_3} \rfloor$ of the mixing-layer integrals of the spanwise MDTKE components from (4.6) and (3.22); (c) and (d) plot, at 30 and 45 ns, respectively, the mixing-layer centre-plane autocorrelation functions (4.8a,b) calculated directly from the finest-resolution simulation (Sim), along with parabolic fits (Fit); (e) plots the $y$- and $z$-direction Taylor microscales $\lambda _{t,i}$ from (4.10); and (f) plots the effective spanwise-direction Taylor Reynolds number $Re_{t,23}$ from (4.12). In (c,d), only values of $R_i$ at $r<2\ \mathrm {\mu }\textrm {m}$ were used to fit the analytic form (4.10). In all legends and headers, numbers state $\mathcal {N}_{yz,2}$ as defined in the caption of table 1, and the letters $y$ and $z$ correspond to the coordinate indices 2 and 3, respectively.

To further explore the structure of the mixing layer in the two spanwise directions, we consider length scales and dimensionless numbers based on the $u_2$ and $u_3$ velocity components, at the mixing-layer centre-plane $x=x_c$ defined by (4.3). We draw on the analytical approaches of Weber et al. (Reference Weber, Haehn, Oakley, Rothamer and Bonazza2014b) in their experimental study of the RM instability and Cook & Dimotakis (Reference Cook and Dimotakis2001) in their computational study of the RT instability. Define the centre-plane $y$- and $z$-direction autocorrelation functions

(4.8a,b)\begin{equation} R_{2}( r)= \frac{ \overline{ u_2( x_c, y, z) u_2( x_c, y + r, z)} }{ \overline{ u_2( x_c, y, z)^2}}, \quad R_{3}( r)= \frac{ \overline{ u_3( x_c, y, z) u_3( x_c, y , z+ r)} }{ \overline{ u_3( x_c, y, z)^2}}. \end{equation}

The corresponding $y$- and $z$-direction Taylor microscales (Tennekes & Lumley Reference Tennekes and Lumley1972; Pope Reference Pope2000; Davidson Reference Davidson2015) can be defined from the curvature of $R_2$ and $R_3$, respectively

(4.9)\begin{equation} \lambda_{t,i} = \left( -\frac{1}{2} \left.\frac{ \textrm{d}^2 R_i}{\textrm{d} r^2} \right|_{r=0} \right)^{-1/2}. \end{equation}

In the present study, we calculate $\lambda _{t,i}$ by fitting a parabola

(4.10)\begin{equation} R_{i,{fit}} = 1 - \frac{r^2}{ \lambda_{t,i}^2} \end{equation}

to $R_i(r)$ over a small interval near $r=0$. Figures 13(c) and 13(d) show examples of $R_i$ and $R_{i,{fit}}$ calculated from the finest-resolution baseline simulation at a late pre-reshock time and a late post-reshock time. The Taylor microscales evoke definitions of centre-plane Taylor Reynolds numbers based on MDTKE-component velocities

(4.11)\begin{equation} Re_{t,i} = \frac{ \bar{\rho} \lambda_{t,(i)}}{ \overline{\mu}} \left( \frac{ 2 \mathcal{K}_{(i)}}{\bar{\rho}} \right)^{1/2}, \end{equation}

with no summation on $i$ implied. An effective spanwise-direction Taylor Reynolds number can also be defined by

(4.12)\begin{equation} Re_{t,23}= \frac{Re_{t,2} + Re_{t,3}}{2}. \end{equation}

Compare (39a) and (39b) in Cook & Dimotakis (Reference Cook and Dimotakis2001).

Figures 13(e) and 13(f) plot time histories of the $y$- and $z$-direction Taylor microscales and the effective spanwise-direction Taylor Reynolds number in the three baseline simulations. For each simulation, the $y$-direction microscale $\lambda _{t,2}$ is always larger than the corresponding $z$-direction microscale $\lambda _{t,3}$. Hence, as also observed in figure 13(b), there are persistent differences in the mixing-layer structure along the two spanwise directions, even at late times. The post-reshock Taylor Reynolds number $Re_{t,23}$ is significantly larger than 100–140, suggesting that the flows meet the criterion proposed by Dimotakis (Reference Dimotakis2000) for transition to a well-mixed turbulent state. (Zhou (Reference Zhou2007) argues that additional, more stringent conditions must also be met for non-stationary flows.) However, both the microscales and Reynolds number are moderately sensitive to the mesh (particularly after reshock) and have not reached their DNS limits in the present study. Any definitive conclusions based on absolute values of $\lambda _{t,2}$, $\lambda _{t,3}$ and $Re_{t,23}$ would require greater resolution and/or more advanced computational schemes, e.g. less-dissipative, higher-order numerics or explicit SGS models. We also reiterate that the turbulent energy spectra of figure 12(b) do not exhibit broad inertial ranges; such ranges may or may not appear with greater resolution and/or more advanced computational schemes.

4.6. Evolution of enstrophy

The vorticity $\omega _i$ is the curl of the velocity $u_i$, i.e. $\omega _i = \epsilon _{ijk} \partial u_k/\partial x_j$, where $\epsilon _{ijk}$ is the antisymmetric Levi-Civita symbol (also called the permutation symbol or the alternating unit tensor). The enstrophy is a scalar measure of vorticity magnitude

(4.13)\begin{equation} \varOmega= \tfrac{1}{2} \omega_i \omega_i . \end{equation}

An evolution equation for $\varOmega$ can be derived from the Navier–Stokes equations (3.3)

(4.14)\begin{align} \frac{ \partial}{\partial t} \left( \rho \varOmega \right) &+ \frac{ \partial}{\partial x_j} \left( \rho \varOmega u_j \right) =\underbrace{ \left( \rho \omega_j S_{ij} \omega_i \right) }_{ \mathbb{E}_{I}} + \underbrace{\left( -2\rho \varOmega \frac{ \partial u_j}{\partial x_j} \right)}_{ \mathbb{E}_{II}} \nonumber\\ &+ \underbrace{ \left( \frac{\omega_i}{\rho} \epsilon_{ijk} \frac{ \partial \rho}{\partial x_j} \frac{\partial p}{\partial x_k} \right) }_{ \mathbb{E}_{III}} + \underbrace{\left( \rho \omega_i \epsilon_{ijk} \frac{\partial}{\partial x_j} \left[ \frac{1}{\rho} \frac{\partial \sigma_{kl}}{\partial x_l} \right] \right)}_{ \mathbb{E}_{IV}} . \end{align}

Vorticity dynamics is important in instability growth, variable-density turbulent mixing and shock-wave–turbulence interaction (Andreopoulos, Agui & Briassulis Reference Andreopoulos, Agui and Briassulis2000; Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2010). Analysis of vorticity and enstrophy in simulated non-HED mixing layers with reshock can be found in Schilling & Latini (Reference Schilling and Latini2010), Grinstein et al. (Reference Grinstein, Gowardhan and Wachtor2011), Hahn et al. (Reference Hahn, Drikakis, Youngs and Williams2011) and Tritschler et al. (Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014), which consider 3-D layers arising from multimode perturbations; Latini et al. (Reference Latini, Schilling and Don2007) and Schilling, Latini & Don (Reference Schilling, Latini and Don2007), which consider 2-D layers arising from single-mode perturbations; and Latini & Schilling (Reference Latini and Schilling2020), which considers both 2-D and 3-D layers arising from single-mode perturbations. Vorticity analysis also features prominently in the works on shock-wave–turbulence interaction by Larsson & Lele (Reference Larsson and Lele2009), Larsson, Bermejo-Moreno & Lele (Reference Larsson, Bermejo-Moreno and Lele2013), Ryu & Livescu (Reference Ryu and Livescu2014) and Livescu & Ryu (Reference Livescu and Ryu2016). The supplementary material provides a complete derivation of (4.14).

The terms on the right-hand side of (4.14) are labelled. The first term $\mathbb {E}_{I}$ is the vortex-stretching term, accounting for stretching, shrinking and tilting of vortices. It vanishes for a 2-D flow. The second term $\mathbb {E}_{II}$ is the enstrophy–dilatation term, which vanishes for a constant-density flow (in which $\partial u_j/\partial x_j \equiv 0$). The third term $\mathbb {E}_{III}$ is the baroclinic source term, which generates enstrophy via misalignment of density and pressure gradients, e.g. in a flow with RM and RT instabilities. The fourth term $\mathbb {E}_{IV}$ is the dissipation term, which captures the dissipative action of viscosity via the viscous stress tensor $\sigma _{kl}$.

Figure 14(a) plots $\lfloor {\rho \varOmega } \rfloor$, the mixing-layer integral of density-weighted enstrophy, versus time. Like $\lfloor {\mathcal {K}} \rfloor$ in figure 12(a), $\lfloor {\mathcal {\rho \varOmega }} \rfloor$ rises slowly after first shock, stabilizes, increases sharply by over one order of magnitude after reshock and finally decays as the mixing layer grows and decompresses. Table 4 quantifies the reshock-induced increase and the post-reshock decay at selected times. From 30 to 35 ns, the integrated density-weighted enstrophy increases by a factor of $\sim$30–70, much like the integrated MDTKE as reported in table 3. The sharp reshock-induced increase in enstrophy seen here is consistent with findings from the non-HED computational studies by Latini et al. (Reference Latini, Schilling and Don2007), Schilling et al. (Reference Schilling, Latini and Don2007) and Tritschler et al. (Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014). As also observed by those authors, enstrophy is much more sensitive to the mesh than turbulent kinetic energy. Indeed, at any given time, $\lfloor {\rho \,\varOmega } \rfloor$ increases with the density of computational zones. The ratio of post-reshock to pre-reshock $\lfloor {\rho \,\varOmega } \rfloor$ is less sensitive to the mesh than instantaneous values of $\lfloor {\rho \,\varOmega } \rfloor$.

Figure 14.

In (a), evolution of the mixing-layer integral of density-weighted enstrophy $\rho \varOmega$ from (4.13) in the three baseline simulations. The legend states $\mathcal {N}_{yz,2}$ as defined in the caption of table 1. In (b,c), evolution of the mixing-layer integrals of each term $\mathbb {E}_{I}, \ldots , \mathbb {E}_{IV}$ in (4.14), normalized by $\lfloor {\rho \varOmega } \rfloor$, in the finest-resolution baseline simulation: (b) displays early-time results before reshock, and (c) displays late-time results after reshock. Note the difference in the ordinate limits of the two figures. Equation (3.22) defines the mixing-layer integral.

Table 4.

Selected values of the mixing-layer integral of density-weighted enstrophy ($\textrm {g}\,\textrm {cm}^{-3}\,\textrm {ns}^{-2}$) at three different times (ns). The ratio of the early post-reshock value to the late pre-reshock value and the ratio of the early post-reshock value to the late post-reshock value are given. Compare with figure 14(a).

To obtain deeper insight into growth and decay of the enstrophy, we examine mixing-layer integrals of each term on the right-hand side of (4.14). When normalized by $\lfloor {\mathcal {\rho \,\varOmega }} \rfloor$, such integrals have units of inverse time and quantify the relative magnitude of the various mechanisms of enstrophy creation or destruction. Figure 14(b) plots the four quantities $\lfloor {\mathbb {E}_{I}} \rfloor /\lfloor {\rho \,\varOmega } \rfloor ,\ldots ,\lfloor {\mathbb {E}_{IV}} \rfloor /\lfloor {\rho \varOmega } \rfloor$ in the finest-resolution baseline simulation before reshock. The baroclinic source term dominates the enstrophy production just after first shock. This early-time production is offset by destruction due to the enstrophy–dilatation term, which is negative in areas of local fluid expansion (where $\partial u_j/\partial x_j > 0$). The contribution from the vortex-stretching term is initially small but grows steadily, eventually outweighing the contribution from the baroclinic source term. Enstrophy destruction due to the dissipation term occurs throughout pre-reshock growth, but its magnitude is small.

Figure 14(c) plots the same four quantities after reshock. As the reshock traverses the mixing layer, there is a sharp increase in enstrophy production due to both the baroclinic source and enstrophy–dilatation terms. The latter term is positive in areas of local fluid compression (where $\partial u_j/\partial x_j < 0$). In fact, $\lfloor {\mathbb {E}_{II}} \rfloor$ exceeds $\lfloor {\mathbb {E}_{III}} \rfloor$ during the reshock traversal, suggesting that fluid compressibility has a significant role in vorticity generation during reshock. After the reshock exits the mixing layer, enstrophy production is principally due to the vortex-stretching term. At very late time, enstrophy–dilatation production increases again as the reflected main shock compresses the mixing layer. The sum of the four mixing-layer integrals $\lfloor {\mathbb {E}_{I}} \rfloor , \ldots , \lfloor {\mathbb {E}_{IV}} \rfloor$ is positive at all times after reshock, suggesting that the post-reshock decay of $\lfloor {\rho \varOmega } \rfloor$ is attributable to numerical dissipation and/or redistributive mechanisms involving the left-hand side of (4.14), e.g. entrainment of zero-enstrophy unmixed fluid into the developing mixing layer. The fact that the mixing layer undergoes intense post-reshock decompression, as discussed in § 4.3 and in contrast to typical non-HED analogues, also points to the importance of advective redistribution.

Figures 14(b) and 14(c) only plot results from the finest-resolution baseline simulation. The supplementary material provides corresponding results from the two lower-resolution baseline simulations. Qualitative trends in the normalized mixing-layer integrals $\lfloor {\mathbb {E}_{I}} \rfloor /\lfloor {\rho \varOmega } \rfloor ,\ldots ,\lfloor {\mathbb {E}_{IV}} \rfloor /\lfloor {\rho \varOmega } \rfloor$ in the two lower-resolution cases are mostly the same as those observed in the finest-resolution case, with one important exception: $\lfloor {\mathbb {E}_{I}} \rfloor /\lfloor {\rho \varOmega } \rfloor$ is significantly smaller with coarser meshes, particularly at late pre-reshock and early post-reshock times. We suspect that the finest-resolution simulation better resolves the 3-D physics of vortex stretching, just as it better resolves the growth of initially small $z$-direction velocity fluctuations as discussed in § 4.5.

Moreover, it is important to emphasize that (4.14) is not exact for an under-resolved 3-D simulation. Numerical errors arise both from the discretization of the governing equations within Ares and from the use of finite differences during post-processing, e.g. to compute $\varOmega$ from $u_i$. For a detailed examination of numerical errors, $\varOmega$ could be replaced in (4.14) with $\breve { \varOmega }+ \mathfrak {T}_{\varOmega }$, where $\varOmega$ is the exact solution to the system of partial differential equations (3.1)–(3.6) with (4.13), $\breve {\varOmega }$ is the numerical solution obtained from the post-processed simulation and $\mathfrak {T}_{\varOmega }$ is the truncation error. Similar substitutions could be made for all flow variables and their derivatives. Then, a new equation could be derived that would resemble (4.14), except for the addition of truncation error terms that vanish in the DNS limit. The procedure outlined here draws on the theory of equivalent or modified differential equations, foundational to computational fluid dynamics (Hirsch Reference Hirsch2007, § 7.1). Truncation errors associated with enstrophy evolution were recently analysed by Zhou, Groom & Thornber (Reference Zhou, Groom and Thornber2020) in simulations of shocked non-HED mixing layers. In the present study, the grid independence of several trends mentioned above – e.g. the larger production during reshock via the enstrophy–dilatation term than via the baroclinic source term – suggests that they are physically meaningful and not merely numerical artefacts. Nevertheless, further analysis of numerical errors is warranted, particularly of their role in the post-reshock enstrophy decay.