1. Introduction
The Richtmyer–Meshkov instability (RMI) (Richtmyer Reference Richtmyer1960; Meshkov Reference Meshkov1969) is generated when a shock passes through a perturbed interface separating two gases with different densities. The instability is formed by misalignment of the pressure gradient across the shock and the density gradient across the interface. This misalignment results in the formation of baroclinic vorticity, which then leads to the growth of perturbations at the interface, and mixing of the two gases (Brouillette Reference Brouillette2002). As these perturbations grow, secondary, shear-driven instabilities develop, which cause smaller-scale mixing to occur (Vorobieff et al. Reference Vorobieff, Tomkins, Kumar, Goodenough, Mohamed and Benjamin2004). See the review by Brouillette (Reference Brouillette2002) for a comprehensive description of the basic physical processes involved in the RMI. The recent reviews by Zhou (Reference Zhou2017a,Reference Zhoub) provide an extensive discussion on the applications, properties and analytical methods for flows exhibiting the RMI as well as the Rayleigh–Taylor and Kelvin–Helmholtz instabilities.
The process of shock-driven mixing by RMI occurs in flows at a range of scales. At the largest scales, shocks from supernova explosions propagate through gas and cosmic dust (Mendis & Rosenberg Reference Mendis and Rosenberg1994; Woitke Reference Woitke2006; Bocchio, Jones & Slavin Reference Bocchio, Jones and Slavin2014), leading to formation of structures on the scale of light-years (Chevalier, Blondin & Emmering Reference Chevalier, Blondin and Emmering1992; Kane, Drake & Remington Reference Kane, Drake and Remington1999) and allowing for the creation of heavy elements. In contrast, inertial confinement fusion experiments produce flows on the scale of micrometres, where shock-driven mixing has the undesirable effect of dissipating energy, thereby reducing fuel compression and lowering fusion yield (Lindl, McCrory & Campbell Reference Lindl, McCrory and Campbell1992). The RMI can also be used to enhance the mixing of gaseous or liquid fuels in scramjet engines (Yang, Kubota & Zukoski Reference Yang, Kubota and Zukoski1993) and is essential for the operation of pulse detonation and rotating detonation engines (Huang et al. Reference Huang, Tang, Li and Zhang2012).
Earlier experimental studies of the RMI have focused on the interaction of a planar shock with a single, sinusoidally perturbed interface (Meshkov Reference Meshkov1969). Variations from this configuration include shock interactions with perturbed, planar interfaces such as those in the experiments of Rasmus et al. (Reference Rasmus2019) and V-shaped interfaces that were studied experimentally by Zhai et al. (Reference Zhai, Dong, Si and Luoa2016). Experimental studies have also been conducted at the Atomic Weapons Establishment on the interaction of a shock with a finite-width interface featuring a ‘chevron’ shaped oblique perturbation by Smith et al. (Reference Smith, Holder, Barton, Morris and Youngs2001), and of half-height interfaces by Holder & Barton (Reference Holder and Barton2004). The related numerical work by Hahn et al. (Reference Hahn, Drikakis, Youngs and Williams2011) studies the effect of additional perturbations on the oblique interfaces. These configurations are essentially two-dimensional (2-D), and the presence of the oblique interfaces leads to shock-driven deposition of vorticity on the interfaces, which leads to shear-driven instabilities.
Additional variations of the RMI include experiments conducted at Los Alamos National Labs, involving a shock interacting with one or more circular cylinders with their axes oriented parallel to the shock (Tomkins et al. Reference Tomkins, Prestridge, Rightley, Marr-Lyon, Vorobieff and Benjamin2003; Vorobieff et al. Reference Vorobieff, Mohamed, Tomkins, Goodenough, Marr-Lyon and Benjamin2003, Reference Vorobieff, Tomkins, Kumar, Goodenough, Mohamed and Benjamin2004; Kumar et al. Reference Kumar, Orlicz, Tomkins, Goodenough, Prestridge, Vorobieff and Benjamin2005). These flow configurations do not feature shear-driven instabilities as their primary mechanism, but they do result in the formation of strong counter-rotating vortex pairs.
The variation of interest in this work has been studied in shock tube experiments at the University of New Mexico conducted by Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b), which involved the interaction of a planar shock with an inclined circular cylinder of heavy gas. This study found that the RMI caused a counter-rotating vortex pair in the cross-sectional plane normal to the cylinder axis. However, on the column surface, it was discovered that the shock caused a Kelvin–Helmholtz instability (KHI) to form. This type of KHI, driven directly by the passage of the shock, has been named the shock-driven KHI (SDKHI). Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b) described an apparent correlation between the Mach number and the KHI wavelength, and proposed a scaling mechanism that could be used to compare experiments at various initial column angles. An additional study by Olmstead et al. (Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a) used the same experimental data to analyse structure functions of scalar intensity maps. The power-law scalings exhibited by these structure functions deviated from those proposed by the Obukhov–Corrsin and Kolmogorov theories (Sreenivasan Reference Sreenivasan1996; Villermaux, Innocenti & Duplat Reference Villermaux, Innocenti and Duplat2001; Celani et al. Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005). However, in-depth characterization of the effects of the SDKHI on the transition to turbulence and mixing was not performed in the previous studies due to limitations of existing experimental techniques caused by the large range of scales involved and the impulsive nature of the RMI and SDKHI.
Popular experimental techniques used to study flows exhibiting RMI and SDKHI include planar laser-induced florescence (PLIF) and particle image velocimetry (PIV). These techniques present challenges to detailed studies of transition and mixing because they are limited to collecting data from a single plane at a given time. Certain important quantities, e.g. pressure fields, also cannot be measured. PLIF requires diluting the test gas with a suitable tracer and has limitations in quantifying planar velocity fields with velocimetry techniques (Palmer & Hanson Reference Palmer and Hanson1994). PIV can provide velocity measurements but can suffer from particles not following the flow at high speed and interfering with the flow physics (Martins et al. Reference Martins, Kirchmann, Kronenburg and Beyrau2021): adding tracer-like particles to shock-accelerated flow can produce shock-driven multiphase instability (Vorobieff et al. Reference Vorobieff, Anderson, Conroy, White, Truman and Kumar2011), and may not capture fine-scale fluctuating velocities. The experimental studies of Olmstead et al. (Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a) in particular could not use PIV at higher Mach numbers due to both particle lag and seeding effects, and could not extract velocity fields using appropriate techniques specifically developed for PLIF, such as image correlation velocimetry (Tokumaru & Dimotakis Reference Tokumaru and Dimotakis1995; Asay-Davis et al. Reference Asay-Davis, Marcus, Wong and De Pater2009) because of insufficient temporal resolution.
Numerical simulations are often used to complement experimental studies by providing additional data inaccessible to experiments. Direct numerical simulations can produce detailed flow quantities at all scales of interest, but such simulations are still unfeasible for these types of flows due to their high computational costs. A more reasonable, commonly used alternative is large-eddy simulations (LES).
Applicability of LES to the flows with RMI was tested by Thornber, Groom & Youngs (Reference Thornber, Groom and Youngs2018), including comparison of the results obtained with different codes. The study found excellent agreement amongst all the codes at high grid resolutions for various considered flow characteristics. Tritschler et al. (Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014) determined that gas composition in the mixing layer was predicted accurately with different LES approaches when compared to experiments. Some of the highest resolution LES of RMI to date have been conducted by Wong, Livescu & Lele (Reference Wong, Livescu and Lele2019). In particular, they found that at late times, after re-shock, the flow resembled that of the Batchelor-type decaying turbulence. The early development of RMI was also considered in high-resolution simulations by Groom & Thornber (Reference Groom and Thornber2019).
However, only a few simulations of RMI appearing in the shock–gas column interaction have been conducted previously (Palekar, Vorobieff & Truman Reference Palekar, Vorobieff and Truman2007; Yang et al. Reference Yang, Kubota and Zukoski1993). These studies have mostly analysed large-scale flow features such as the interface width and the growth rate for 2-D circular gas clouds.
In our group, implicit LES (ILES) have been conducted for 2-D curtains (Romero et al. Reference Romero, Poroseva, Vorobieff and Reisner2021a,Reference Romero, Poroseva, Vorobieff and Reisnerc) and three-dimensional (3-D) heavy-gas columns (Romero et al. Reference Romero, Vorobieff, Poroseva and Reisner2021b, Reference Romero, Poroseva, Vorobieff and Reisner2022) to verify the mechanism leading to SDKHI and to analyse the effect of the flow dimensionality on the shock–column interaction. Simulations were run at Mach numbers ranging from 1.13 to 2.0. The gas curtains and columns were inclined at angles varying from 
$10^\circ$ to 
$45^\circ$ with respect to the shock. The simulation results showed that the vorticity field produced by the shock after it passed through the gas column is consistent with a shear-driven flow, proving that the waves formed on the column surface were, in fact, Kelvin–Helmholtz waves. These wavelengths were found to be sensitive to the post-shock column diameter and were thus influenced implicitly by parameters such as the Mach number and the initial column/curtain angle. Furthermore, while SDKHI was observed to occur in both 2-D and 3-D simulations, 3-D effects were found to be important for the correct reproduction of the flow morphology and the interface growth rates seen in experiments.
The grid resolutions used in our previous simulations were not sufficient to conduct a detailed statistical characterization of the effects of SDKHI on transition to turbulence and mixing. The present study aims to overcome this limitation by conducting high-resolution ILES for the shock passing through a 3-D inclined cylindrical column of heavy gas. The simulation results are used to investigate the effect of SDKHI on the turbulent transition and mixing that occur in such a flow. In the simulations, a column of heavy gas (sulfur hexafluoride, 
${\rm SF}_6$) surrounded by air is inclined at a range of angles 
$\alpha_0 = 1^\circ,5^\circ,10^\circ,30^\circ$ with respect to an oncoming shock with Mach number 2.0. The considered range of angles and the Mach number are chosen to facilitate comparison with the experiments of Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b). The ensuing analysis compares mixing statistics with an emphasis on the mixing width and other integral mixing metrics. Energy spectra and anisotropy are also analysed for each case. In addition, scalar statistics in the vertical plane are compared to published experimental data to investigate the anomalous scaling reported by Olmstead et al. (Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a).
Simulations were performed with FIESTA, a graphics processing unit (GPU) accelerated flow solver developed at the University of New Mexico (Romero Reference Romero2021) for multi-species compressible flows.
2. Numerical method and set-up
2.1. Numerical method
In this study, results were obtained using the 3-D, two-species, compressible, viscous transport equations in their conservative form. These equations consist of the continuity equations for each gas species, equations for each momentum component, and the conservation equation for specific total energy. The resulting system is represented by the following three equations:
$$\begin{gather} \frac{\partial\rho Y_{i}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho \boldsymbol{u}Y_{i}) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \rho \boldsymbol{u}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho \boldsymbol{uu^{\rm T}} + p\boldsymbol\delta) = \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\tau}, \end{gather}$$
$$\begin{gather}\frac{\partial\rho e_{t}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \left[(\rho e_t + p)\boldsymbol{u}\right] = \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{u} - \boldsymbol{q}). \end{gather}$$
Here, 
$\rho$ is the density of the mixture; 
$Y_i$ is the mass fraction of species 
$i$, where 
$Y_i$ is 
$Y_a$ for air and 
$Y_s$ for SF
$_6$; 
$\boldsymbol {u}$ is the 3-D velocity vector 
$\boldsymbol {u} = [u_1,u_2,u_3]$, where 
$1$, 
$2$ and 
$3$ are the streamwise (
$x$), vertical (
$y$) and spanwise (
$z$) directions, respectively; 
$p$ is the static pressure; 
$e_t$ is the specific total energy; 
$\boldsymbol {\tau }$ is the viscous stress tensor; and 
$\boldsymbol {q}$ is the conductive heat flux vector.
The pressure of the gas mixture is computed with Dalton's mixing law (Dalton Reference Dalton1802) assuming ideal gas behaviour:
\begin{equation} p = \rho \left( Y_s R_s T + (1-Y_s) R_a T \right). \end{equation}The mixture temperature is obtained from the internal energy as follows:
\begin{equation} T = \frac{e}{C_{v_{mix}}}, \end{equation}
where 
$e$ is the internal energy, and the mixture specific heat 
$C_{v_{mix}}$ is computed from the species specific heats and mass fractions as
\begin{equation} C_{v_{mix}} = Y_sC_{v,s} + (1-Y_s)C_{v,a}. \end{equation}Expressions for the specific heats of individual species are
\begin{equation} C_{v,a} = \frac{R_{a}}{\gamma_{a}-1},\quad C_{v,s} = \frac{R_{s}}{\gamma_{s}-1}. \end{equation}The gas constants for each individual species are
\begin{equation} R_{a} = \frac{R}{M_{a}},\quad R_{s} = \frac{R}{M_{s}}, \end{equation}
where 
$R$ is the universal gas constant, 
$M_{a}$ is the molecular weight of air, and 
$M_{s}$ is the molecular weight of SF
$_6$.
The viscous stress tensor 
$\boldsymbol {\tau }$ is
\begin{equation} \boldsymbol{\tau} = 2\mu\boldsymbol{S} + \tfrac{1}{3}\, \mu\boldsymbol{\delta}(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}), \end{equation}
where 
$\mu$ is the mixture viscosity computed as 
$\mu =Y_{a} \mu _{a} + Y_{s} \mu _{s}$, and the strain rate tensor is given by
\begin{equation} \boldsymbol{\mathsf{S}} = \tfrac{1}{2}\left[\boldsymbol{\nabla}\boldsymbol{u} + (\boldsymbol{\nabla}\boldsymbol{u})^{{\rm T}}\right]. \end{equation}Subsection 3.1 presents a detailed discussion of the relevant viscous and scalar dissipation length scales. Analysis of the results from the simulations indicates that the viscous length scales may be resolved partially by the grid, while length scales relevant to mass diffusion are not resolved. Therefore, the omission of mass diffusion terms in the governing equations presented above is appropriate for the flow regime under consideration.
The code solves fully dimensional equations. The specific values of parameters utilized in the code are 
$R = 8.314462\ {\mathrm {J}}\ \mathrm {kg}^{-1}\ \mathrm {mol}^{-1}$, 
$M_{a} = 0.028966\ {\mathrm {kg}}\ {\mathrm {mol}^{-1}}$, 
$M_{s} = 0.14606\ {\mathrm {kg}}\ {\mathrm {mol}^{-1}}$, 
$\mu _{a} = 2.928\times 10^{-5}\ {\mathrm {kg}}\ {\mathrm {m}}^{-1}\ {\mathrm {s}}^{-1}$, 
$\mu _{s} = 1.610\times 10^{-5}\ {\mathrm {kg}}\ {\mathrm {m}}^{-1}\ {\mathrm {s}}^{-1}$, 
$\gamma _{a} = 1.402$ and 
$\gamma _{s} = 1.092$.
Advective terms are approximated with the fifth-order weighted essentially non-oscillatory (WENO5) finite difference scheme (Jiang & Shu Reference Jiang and Shu1996; Ramani, Reisner & Shkoller Reference Ramani, Reisner and Shkoller2019a,Reference Ramani, Reisner and Shkollerb). The pressure gradient term and the strain rate tensor are approximated using a fourth-order central difference scheme. The time scheme used is a low-storage, second-order, explicit Runge–Kutta integrator (Williamson Reference Williamson1980).
2.2. Computational domain and discretization
Figure 1(a) depicts the computational domain used for the inclined gas column simulations. In the figure, 
$L_x=0.4$ m is the length of the domain, 
$L_y=0.08$ m is its height, and 
$L_z=0.05$ m is its width.
Figure 1. (a) Domain dimensions and boundary conditions for the inclined gas column simulations. (b) Illustration of the shock, the gas column, and the key cross-sectional planes.
The simulations used a uniform 
$10\,000\times 2000\times 1250$ Cartesian grid with 
${{\rm d}x}={{\rm d}y}={\rm d}z=4.0\times 10^{-5}$ m. This resolution exceeds that of the simulations of Palekar et al. (Reference Palekar, Vorobieff and Truman2007), and is similar to the values used in Wong et al. (Reference Wong, Livescu and Lele2019) and Groom & Thornber (Reference Groom and Thornber2019).
The time step 
${\rm d}t=1\times 10^{-8}$ s was chosen so that the CFL number, based on the maximum wave speed, did not exceed 
$0.1$. Simulations were conducted for 80 000 time steps, which is equivalent to 
$0.8$ ms, similar to the experiments of Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b).
2.3. Initial conditions
Solution of (2.1)–(2.3) requires defining pre-shock and post-shock initial states. In this study, the pre-shock state was defined as atmospheric pressure and room temperature similar to the conditions observed in the experiments (Olmstead et al. Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b). The post-shock conditions in air were determined from the inviscid normal shock relations for Mach number 
$M_{0}=2.0$. Table 1 lists the pre- and post-shock conditions used to initialize these simulations. Throughout this paper, post-shock quantities are indicated with the subscript ‘
${post}$’ where necessary to remove ambiguity.
Table 1. Pre- and post-shock values for the initial conditions on the primary variables in (2.1)–(2.3).
The gas column is positioned ahead of the shock at an angle 
$\alpha _0$ to the shock, with the distribution of SF
$_6$ mass fraction, 
$Y_s$, across the diameter of the gas column described by a Gaussian-based distribution with maximum concentration 
$1.0$ at 
$x=x_0$: 
$Y_s=\exp [-(x-x_0)^2/2\sigma ^2]$, where 
$x_0$ is the location of the centre of the cylindrical column. In experiments, steady-state initial conditions were established and well characterized (Wayne et al. Reference Wayne, Olmstead, Vorobieff, Truman and Kumar2015; Olmstead et al. Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a,Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieffb), allowing value 
$\sigma =1.71$ to be chosen so that the initial gas distribution used in the current simulations matched closely that seen in experiments. In the absence of well-characterized experimental initial conditions, or when steady-state initial conditions cannot be established experimentally, an alternative technique would be to pre-compute an appropriate interface, such as was done in the gas curtain experiments of Gowardhan & Grinstein (Reference Gowardhan and Grinstein2011). The heavy-gas column diameter 
$\delta _0$ is 
$10.4$ mm as measured at the location where the mass fraction of SF
$_6$ has value 0.01. Figure 2 compares the gas distributions used in the simulations and in the experiments.
Figure 2. Comparison of the gas species distribution used in the simulations with the experimental data from Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b). The location of the column diameter 
$\delta$ is indicated by vertical dashed lines.
2.4. Boundary conditions
Boundary conditions for (2.1)–(2.3) include inflow/outflow boundaries (
$BC_{IO}$) at 
$x=0$, 
$x=L_x$, 
$z=0$, 
$z=L_z$, and the reflective wall conditions (
$BC_{W}$) at 
$y=0$ and 
$y=L_y$ as seen in figure 1(a).
$BC_{IO}$ conditions are defined as
\begin{equation} \partial_x\rho(x,t)=0, \quad \partial_x\rho\,u(x,t)=0, \quad \partial_x\rho\,v(x,t)=0, \quad \partial_x \rho\,e_t(x,t)=0, \end{equation}
where 
$\partial _x={\partial }/{\partial x}$ and 
$\partial _y={\partial }/{\partial y}$.
$BC_{W}$ conditions are defined as
\begin{equation} \partial_y\rho(x,t)=0,\quad\partial_y\rho\,u(x,t)=0,\quad\rho\,v(x,t)=0,\quad\partial_y \rho\, e_t(x,t)=0. \end{equation}3. Results
The initial density distribution in the domain is shown in figure 3 for the centreline and vertical planes (as defined in figure 1a). The green region to the left corresponds to the shock wave and post-shock conditions, while the ambient, pre-shock conditions are seen in blue. Pressure is distributed uniformly within each of the two regions. The inclined gas column is seen to the right of the shock in the ambient region. The gas column is inclined at angle 
$\alpha _0=30^\circ$ relative to the shock, similar to the experimental set-up used by Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b).
Figure 3. Initial conditions of the density in (a) the centreline plane and (b) the vertical plane used for simulations of a shock interacting with the gas column initially inclined at 
$\alpha _0=30^\circ$.
As the simulation progresses, the shock travels to the right through the heavy-gas cylinder. Figures 4 and 5 depict the interaction of the shock with the cylinder at various times. The dimensionless time is defined as 
$\tau =\kappa Au_{post}(t-t_0)$ (Olmstead et al. Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b). Here, 
$\kappa =2{\rm \pi} /\delta _0$ is the wavenumber based on the cylinder diameter in the experiments 
$\delta _0$, 
$A=(\rho _{\textrm {SF}_6}-\rho _{air})/(\rho _{\textrm {SF}_6}+\rho _{air})$ is the Atwood number, 
$u_{post}$ is the post-shock velocity, and 
$t_0$ is defined as the time when the shock interaction is halfway down the length of the initial column and corresponds to 
$\tau =0$.
Figure 4. Density plots of the vertical plane at various dimensionless times 
$\tau$ for a Mach 2.0 shock interacting with an inclined circular column of heavy gas at initial angle 
$\alpha _0=30^\circ$ to the plane of the shock. Arrows indicate the transmitted (
$T$) and reflected (
$R$) shocks. Times are (a) 
$\tau =5.7$, (b) 
$\tau =20.7$, (c) 
$\tau =35.7$, (d) 
$\tau =50.7$ and (e) 
$\tau =110.8$.
Figure 5. Pressure plots of the vertical planes at various dimensionless times 
$\tau$ for a Mach 2.0 shock interacting with an inclined circular column of heavy gas at initial angle 
$\alpha _0=30^\circ$ to the plane of the shock. Arrows indicate the transmitted (
$T$) and reflected (
$R$) shocks. Times are (a) 
$\tau =5.7$, (b) 
$\tau =20.7$, (c) 
$\tau =35.7$, (d) 
$\tau =50.7$ and (e) 
$\tau =110.8$.
At 
$\tau =5.7$, the shock has changed the inclination of most of the column to be 
$\alpha _{post}=16.2^\circ$. Due to differences in the density of the heavy gas and the air, the shock is refracted and bent inwards where it intersects the column (figures 4a and 5a). The shock becomes planar again after passing through the column. The column remains at this angle, 
$\alpha _{post}=16.2^\circ$, defined as the post-shock angle, for the remainder of its development.
High-pressure regions are seen behind the main shock and behind a reflected shock as seen in figure 5(a). These high-pressure regions intersect at the foot of the column and are likely to be responsible for an instability observed to form at the foot of the column. The instability is more evident by 
$\tau =20.7$ (figure 4b). Similar phenomena were observed in 2-D simulations of SDKHI (Romero et al. Reference Romero, Poroseva, Vorobieff and Reisner2021a). In the 3-D flow, this phenomenon still occurs, but is less pronounced.
After passage of the shock through the column, variations in the pressure field diminish until later times where pressure variations are coincident with spots of high vorticity.
The surface of the gas column begins to move upwards along its axis at 
$\tau =20.7$. At a later time, 
$\tau =35.7$, perturbations appear on the leading edge of the column. These perturbations first appear near the top of the column and progress quickly downwards towards the column foot. These perturbations are similar to the Kelvin–Helmholtz (KH) waves observed in the experiments of Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b). In the experiments, the KH perturbations formed first on the column leading edge as they do in the 3-D simulations. The KH waves in 3-D simulations and in the experiments also have similar wavelengths and develop at the same times. The KH wavelengths 
$\lambda _{KH}$ were identified by measuring the average distance between wave peaks. Wavelengths are well defined and have regular spacing for each column angle between 
$\tau =20$ and 
$\tau =25$ before additional secondary instabilities begin to appear. See Romero et al. (Reference Romero, Poroseva, Vorobieff and Reisner2022) for a detailed description of wavelength measurements and comparisons with experiments.
When results from the current 3-D simulations (figure 6a) are compared to those from the 2-D simulations of Romero et al. (Reference Romero, Poroseva, Vorobieff and Reisner2021a) (figure 6b), there are notable differences. The lengths of KH waves are much smaller in the 3-D case (table 2), and they form at first near the top of the column. In the 2-D case, the wavelengths are larger and appear first near the column midpoint and at later times (figure 6b). The 3-D case also contains finger-like features downstream of the column. These features are caused by the counter-rotating vortex pair moving material into the vertical plane. Thus while 2-D simulations are important to emphasize or isolate certain flow features, 3-D simulations provide a more accurate description of the flow.
Figure 6. Comparison of the KH waves in the vertical plane of (a) the 3-D heavy-gas column at 
$\tau =28.7$, and (b) a quasi-2-D heavy-gas curtain at 
$\tau =73$ from Romero et al. (Reference Romero, Poroseva, Vorobieff and Reisner2021a).
Table 2. 2-D and 3-D average KH wavelengths normalized by pre-shock diameter 
$\lambda _{KH}/\delta _0$ for different angles.
The turbulent state of the flow continues to evolve to 
$\tau =50.7$. At later time 
$\tau =110.8$, larger structures begin to dominate the flow.
Figure 7 depicts experimental results from Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b). This image shows pixel intensity measurements collected with PLIF diagnostics of the gas column seeded with an acetone tracer; it can be compared to those in figure 4. In particular, the qualitative similarities near 
$\tau \approx 20$ and 
$\tau \approx 35$ are noted and show the similar timing of simulations with experiments.
Figure 7. Experimental data for a Mach 2.0 shock interacting with a cylindrical column of heavy gas with initial tilt angle 
$30^{\circ }$, reproduced with permission from Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b).
Quantitative comparisons of KH wavelengths with experiments can be made using the wavelength selection mechanism introduced in Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b):
\begin{equation} \lambda_{KH} = D_{C}\tan(\alpha_{2})\,M^{1/2}A^{1/2}, \end{equation}
where 
$D_{C}\propto \delta _0$ is the compressed column diameter defined in Olmstead et al. (Reference Olmstead, Wayne, Yoo, Kumar, Truman and Vorobieff2017b), 
$\alpha _{2}$ is the post-shock angle of the column, 
$M$ is the Mach number, and 
$A$ is the Atwood number. Figure 8 compares mean wavelengths normalized by (3.1) for each column angle to mean wavelengths from experimental data. There is exceptional agreement at 
$\alpha _{0}=30^{\circ }$ where the data points overlap. However, this particular scaling does not seem to hold for smaller angles. For the cases 
$\alpha _{0}=1^{\circ }$ and 
$\alpha _{0}=5^{\circ }$, the wavelengths measured in simulations were identical and were resolved only by about five grid cells. This may indicate that the simulation resolution is insufficient to resolve the short wavelengths that occur at smaller initial tilt angles and is not necessarily a failing of the proposed wavelength selection mechanism. For detailed comparisons of simulations to experiment, see Romero et al. (Reference Romero, Vorobieff, Poroseva and Reisner2021b, Reference Romero, Poroseva, Vorobieff and Reisner2022).
Figure 8. Normalized KH wavelengths at different angles for experiments (
$\times$) and simulations (
$\bullet$).
3.1. Statistics
Time averaging is not applicable in this flow to collect statistics for turbulence characterization. However, the central portion of the column appears to be homogeneous in the direction of the column tilt. This portion of the flow is used in the current study for collecting statistical data. Data are averaged at the vertical plane through the column axis and in the direction of the mean post-shock tilt angle of the column. For the analysis below, the flow field was cropped vertically to retain the centre third of the column. Horizontally, the domain was cropped to include the area where the SF
$_6$ mass fraction first exceeds 0.001. Figure 9 shows the cropped flow field in the vertical plane for the 
$30^\circ$ column at 
$\tau =20.7$, along with an axis indicating the direction in which averages were taken.
Figure 9. Vertical cross-section of the cropped density field for a 
$30^\circ$ column used for collecting statistics at 
$\tau =35.7$. The dashed line indicates the mean post-shock tilt angle along which averages were taken.
Favre averaging (Favre et al. Reference Favre, Kovasznay, Dumas, Gaviglio and Coantic1976) is used when extracting statistics. In this flow representation, the velocity field is decomposed into the mean and fluctuating components: 
$u_i = \widetilde {u_i} + u''_i$, where 
$\widetilde {u_i}={\left \langle \rho u_i\right \rangle }/{\left \langle \rho \right \rangle }$ is the Favre-averaged mean velocity, and 
$u''_{i}$ is the Favre-averaged velocity fluctuation. The spatial average of a variable 
$\phi$ is defined by
\begin{equation} \left\langle\phi\right\rangle=\frac{1}{N}\sum_i \phi_i, \end{equation}
where 
$N$ is the number of computational cells in the direction of the column tilt. For these simulations, 
$N=650$ for the statistical region described above.
Below, data are presented with respect to the 
$x$-coordinate normalized as
\begin{equation} x^*(y) = \frac{x-x_0(y)}{w}, \end{equation}
where 
$x_0(y)$ is the location of the leading edge of the column and depends on the 
$y$-coordinate due to the post-shock angle, and 
$w$ is the width of the column defined by the locations at the leading and trailing edges where the mass fraction is 0.001.
Profiles of the mean mass fraction 
$\left \langle Y_s\right \rangle$ and of the mean velocity components 
$\widetilde {u_{i}}$ were computed at different times in the streamwise direction and are shown in figure 10. In the figure, the mean velocity components are normalized by the post-shock velocity 
$u_{post}$.
Figure 10. Mean mass fraction and mean velocity profiles for a column with initial tilt 
$\alpha _0=30^\circ$ at different times 
$\tau$: 5.7 (solid line), 20.7 (dashed line), and 50.7 (dotted line).
Figure 10(a) shows that at early time, the concentration of SF
$_6$ is highest on the leading edge. As time progresses, the SF
$_6$ concentration becomes more uniform inside the column.
Figures 10(b)–10(d) show that the flow motion in the vertical plane of the gas column occurs mainly in the 
$x$- and 
$y$-directions. In the streamwise direction, the shock accelerates the gas near the edges of the column to a higher velocity than the gas in the interior of the column. With time, the minimum value of 
$u_1$ moves within the column, finally localizing between the column centre and the leading edge, as shown in figure 10(b) at 
$\tau = 50.7$. After that time, the profile of this velocity component tends to the post-shock velocity (the value 1.0 in figure 10b) at all locations (not shown here).
In the 
$y$-direction, the gas column inclination causes the dense gas to accelerate upwards almost everywhere through the column width (figure 10c). The less dense gas at the column edges moves downwards. As time progresses, the dense gas motion in the area between the column centre and its trailing edge reverses its direction, with all gas moving downwards. At 
$\tau = 20.7$, only the gas near the leading edge is still moving upwards. By 
$\tau = 50.7$, there is almost no gas motion in the vertical direction.
The velocity in the 
$z$-direction, 
$u_3$, is close to zero at all locations within the column and at all times (figure 10d). This is due to the column cross-sectional symmetry being preserved with time for the chosen sample plane.
The column inclination angle was found to affect the gas 
$x$- and 
$y$-velocity components (figure 11), but not the velocity component in the 
$z$-direction. Specifically, in the 
$x$-direction, columns with initial angles up to 
$\alpha _0=10^\circ$ have similar profiles (figure 11a). The gas is accelerated mostly at the column leading edge and in an area behind the trailing edge. At angle 
$30^\circ$, the areas of high gas velocity are more narrow, with the second peak shifting closer to the column trailing edge. At any column angle, the gas at the trailing edge accelerates less than the rest of the column when experiencing a shock.
Figure 11. Mean (a) 
$x$-velocity and (b) 
$y$-velocity at 
$\tau \approx 5.7$ for columns with different initial angles: 
$1^\circ$ (solid line), 
$5^\circ$ (dashed line), 
$10^\circ$ (dash-dotted line), and 
$30^\circ$ (dotted line).
In the 
$y$-direction, the gas was found to accelerate upwards to a higher velocity for larger initial angles 
$\alpha _{1}$, as shown in figure 11(b). At any angle, the gas velocity in the 
$y$-direction is less than in the 
$x$-direction.
Overall, we can infer that the gas motion in such a flow geometry is mainly 2-D in the vertical plane, and with time becomes predominantly one-dimensional (in the 
$x$-direction).
Normalized Reynolds stresses are presented in figure 12 as 
${\sqrt {{\mathsf{R}}_{ii}}}/{u_{post}}$, where the Favre-averaged Reynolds stresses are defined as
\begin{equation} {\mathsf{R}}_{ij} = \frac{\left\langle\rho u_i'' u_j''\right\rangle}{\left\langle\rho\right\rangle}. \end{equation}
Figure 12. Evolution of the Reynolds stresses with time for a column with initial tilt 
$\alpha _0=30^\circ$ at different times 
$\tau$: 5.7 (solid line), 20.7 (dashed line), and 50.7 (dotted line).
These figures show that initially, the intensities of all fluctuating velocities and of the shear stresses have a maximum near the trailing edge of the column. Over time, the peak values of all Reynolds stresses move towards the leading edge of the column (
$x^*=0$), with stresses generally growing in time across the entire column.
The Reynolds number describing this flow must be defined using the most physically relevant velocity and length scales. The flow features observed and characterized in this work evolve after the shock interaction with the density interface produces baroclinic vorticity. Therefore, neither the piston velocity of the shocked flow nor the shock front speed are relevant – the patterns observed are produced by the fluctuating components of the velocity field due to the deposition of baroclinic vorticity. Here, this characteristic velocity is represented by the average magnitude of fluctuating velocity, 
$u_{turb}=\sqrt {R_{11}+R_{22}+R_{33}}$. The appropriate length scale can be determined from the same considerations. According to Richtmyer (Reference Richtmyer1960), the initial (linearized) perturbation amplitude growth 
$a(t)$ is proportional to the instability wavenumber 
$\kappa$ and to the initial instability amplitude. For a fixed amplitude, the instability growth is faster for higher wavenumbers (and smaller wavelengths). The feature that is associated with the highest initial amplitude and the lowest wavelength in the centreline plane is the heavy-gas column diameter 
$\delta _0$, and in the vertical plane it is the KH wavelength 
$\lambda _{KH}$ (which is related to 
$\delta _0$ via (3.1)).
For the cases considered, the characteristic velocity is 
$u_{turb}=65.1\ {\rm m}\ {\rm s}^{-1}$, the mean post-shock density is 
$\rho _{post}=3.6\ {\rm kg}\ {\rm m}^{-3}$, and the mean mixture viscosity is 
$\mu _{mix}=2.65\times 10^{-5}\ {\rm kg}\ {\rm m}^{-1}\ {\rm s}^{-1}$. For the characteristic length 
$\delta _0$, the Reynolds number 
$Re_\delta = {\rho _{mix}\, u_{turb} \, \delta _0}/{\mu _{mix}}$ can be estimated as 
$\sim$90 000. This is perhaps the most realistic estimate, because it is associated with the scale that also produces the strongest vorticity deposition in the form of counter-rotating vortex pairs in the centreline plane. In direct experimental measurements of the flow field in a similar geometry (periodically perturbed gas curtain with wavelength 
$\lambda =6\ {\rm mm}$) but at an appreciably lower Mach number 1.2, the Reynolds number varied in the range 
$10\,000 < Re < 22\,500$ (Prestridge et al. Reference Prestridge, Rightley, Vorobieff, Benjamin and Kurnit2000).
However, the present study is concerned with the flow evolution in the vertical plane, so the scale 
$\lambda _{KH}$, varying in our simulations from 
$0.66\ {\rm mm}$ to 
$1.68\ {\rm mm}$, becomes relevant. An estimate for the Reynolds number based on such a length scale, 
$Re_\lambda = {\rho _{mix}\, u_{turb}\,\lambda _{KH}}/{\mu _{mix}}$, is in the range between 6000 and 17 000.
Information about the Reynolds number range can be used further to estimate characteristic length scales for viscous and mass diffusion processes. For viscous processes, the Taylor and Kolmogorov length scales can be approximated by 
$\lambda _{T} = \lambda _{KH} \, Re_{\lambda _{KH}}^{-1/2}$ and 
$\eta = \lambda _{KH} \, {Re}_{\lambda _{KH}}^{-3/4}$, where 
$\lambda _{T}$ is the Taylor microscale, and 
$\eta$ is the Kolmogorov microscale. For the range of Reynolds numbers 
${Re}_{\lambda _{KH}}$ identified above (going with the more conservative estimate), this results in 
$\lambda _{T}=10\ \mathrm {\mu } {\rm m}$ to 
$40\ \mathrm {\mu } {\rm m}$, and 
$\eta =1\ \mathrm {\mu } {\rm m}$ to 
$10\ \mathrm {\mu } {\rm m}$. The grid spacing is 
${{\rm d}x}=40\ \mathrm {\mu } {\rm m}$, which indicates that the viscous dissipation range may be partially resolved, so viscous effects must be accounted for.
The scalar dissipation length scale is relevant to the mass diffusion process and can be computed as 
$\eta _{D} = \eta \, {Sc}^{-3/4}$, where 
${Sc}$ is the Schmidt number (Sreenivasan Reference Sreenivasan2019). The Schmidt number is defined as 
${Sc}={\mu _{mix}}/({\rho _{mix} \, D})$, where 
$D=2.23\times 10^{-6} \ {\rm m}^2\ {\rm s}^{-1}$ is the mean diffusion coefficient between air and SF
$_6$ at the post-shock temperature and pressure, resulting in 
$Sc=2.8$. These estimates give a range for the mass diffusion length scale 
$\eta _{D}$ between 
$0.5\ \mathrm {\mu } {\rm m}$ and 
$5\ \mathrm {\mu } {\rm m}$. This is smaller than the grid resolution 
${{\rm d}x}=40\ \mathrm {\mu } {\rm m}$, so the length scales at which mass diffusion processes occur are not resolved. Notably, the Péclet number 
${Pe}={Sc}\,{Re}_{\lambda _{KH}}$ ranges from 17 000 to 48 000, which is much larger than unity, indicating that the effects of mass diffusion occur on a significantly larger time scale than those of the advective transport. For these reasons, the effects of mass diffusion are not included in this study and can be disregarded safely on the time scales of interest.
3.2. Anisotropy
In this subsection, the flow anisotropy is studied. For this purpose, the anisotropy tensor is defined as
\begin{equation} b_{ij}=\frac{{\mathsf{R}}_{ij}}{{\mathsf{R}}_{kk}}-\frac{1}{3}\,\delta_{ij}, \end{equation}
where 
$\delta _{ij}$ is the Kronecker delta, and 
${\mathsf{R}}_{ij}$ is the Favre-averaged Reynolds stress tensor defined above.
The diagonal components of the anisotropy tensor give an estimate of the contribution of each fluctuating velocity component to turbulent kinetic energy. These values range from 
$-1/3$, corresponding to zero energy from that component, to 
$2/3$, corresponding to all energy coming from that component. If all three diagonal components of the anisotropy tensor 
$b_{ij}$ are zero, then the flow is isotropic.
Figure 13 shows diagonal components of the anisotropy tensor for different initial column angles at approximately the same time 
$\tau = 85.1$. Turbulent energy is distributed non-uniformly through the column at each angle. For the smaller angles (
${\leq }10^\circ$), the flow is anisotropic throughout, dominated by the streamwise velocity fluctuations near the leading and trailing edges of the column, and by the spanwise fluctuations in the column interior. For the 
$30^\circ$ column, the flow becomes more isotropic in the interior of the column, with the transverse and vertical fluctuations having equal contributions to turbulent kinetic energy.
Figure 13. Anisotropic tensor components for columns with initial angles (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$. The components are 
$b_{11}$ (solid line), 
$b_{22}$ (dashed line) and 
$b_{33}$ (dotted line).
In addition to examining the components, invariants of the anisotropy tensor can be plotted to examine the turbulence state over time. The anisotropy tensor invariants 
$\eta$ and 
$\xi$ are defined as
$$\begin{gather} 6\eta^{2}=b_{ij}b_{ji}, \end{gather}$$
$$\begin{gather}6\xi^{3}=b_{ij}b_{jk}b_{ki}, \end{gather}$$
where 
$b$ is the anisotropy tensor defined in (3.5), and repeated indices indicate summation. For realizable turbulence, 
$\eta$ and 
$\xi$ fall within the Lumley triangle, with the sides and vertices of the triangle representing special states of the Reynolds stress tensor (Pope Reference Pope2000).
Anisotropy-invariant maps are plotted for each column initial angle in figure 14. Figure 14(a) provides an overview of the Lumley triangle with an invariant map at three locations (
$x^{*}=0.1,0.5,0.9$) within the 
$30^{\circ }$ column. For all cases, the trajectories are contained within only a small portion of the Lumley triangle, so figures 14(b–e) are restricted to the region indicated by the dotted box in figure 14(a). Arrows indicate the direction of the invariant trajectories with time.
Figure 14. Reynolds-stress-invariant maps with (a) an overview of the Lumley triangle and (b–e) columns with different initial angles for three locations within the column: (b) 
$1^\circ$, (c) 
$5^\circ$, (d) 
$10^\circ$ and (e) 
$30^\circ$. Circles indicate 
$x^{*}=0.1$, triangles indicate 
$x^{*}=0.5$, and squares indicate 
$x^{*}=0.9$.
For each case, the initial values of 
$\eta$ and 
$\xi$ are near the limit of single-component turbulence as fluctuations are initially present in only the 
$x$-direction (figure 13). As time progresses, the flow behaviour changes in the direction of isotropy while maintaining a state close to that of axisymmetric turbulence.
3.3. Mixedness
The molecular mixedness ratio over time can also be compared for columns with different initial angles. In the current study, molecular mixedness is defined as
\begin{equation} \theta =\frac{\left\langle\chi_s(1-\chi_s)\right\rangle_V}{\left\langle\chi_s\right\rangle_V\left\langle(1-\chi_s)\right\rangle_V}, \end{equation}
where 
$\left \langle {\cdot }\right \rangle _V$ is the volumetric mean defined below, and 
$\chi _s$ is the mole fraction of SF
$_6$. The mole fraction is given by 
$\chi _s={n_{s}}/({n_{s}+n_{a}})$, where 
$n_{i}={\rho Y_{i}}/{M_{i}}$ is the molar count of species 
$i$ in a given grid cell with mass fraction 
$Y_{i}$ and molecular mass 
$M_{i}$. The volumetric mean is defined as
\begin{equation} \left\langle\phi\right\rangle_V = \frac{1}{N_xN_yN_z}\sum_i\sum_j\sum_k \phi_{ijk}, \end{equation}
where 
$N_x$, 
$N_y$ and 
$N_z$ are the numbers of grid cells in the 
$x$-, 
$y$- and 
$z$-directions, respectively, and 
$\phi _{ijk}$ is the value at the grid cell with coordinates 
$(i,j,k)$.
A value 
$\theta =1$ means that the flow is mixed as much as possible within the volume, which means that the SF
$_6$ is distributed uniformly.
Figure 15 shows the molecular mixedness ratio over time for different initial column angles. For each angle, the column has initial molecular mixedness 
$\theta =0.62$ due to the diffuse interface across the column defined by the initial conditions (see figure 2). Over time, the mixedness tends to value 1 for each initial angle, while the initial angle has very little effect on the mixing rate. At times beyond 
$\tau =60$, large-scale structures begin to impact the statistical region, and this may affect the mixedness behaviour. The apparent deviation for the 
$10^{\circ }$ column near 
$\tau \approx 0$ is a numerical artefact due to the presence of the shock wave within the statistical region under consideration that happens in this particular case. The shock is not parallel to the column tilt angle and disrupts the averaging procedure in (3.8). In other cases, solution data files do not capture the shock wave within the statistical region due to the I/O frequency of the simulation.
Figure 15. Volume mean molecular mixedness ratio of SF
$_6$ over time for different initial column angles.
3.4. Spectra
Figure 16 shows energy spectra at different times for heavy-gas columns with various initial tilt angles. Here, spectral energy is defined as 
$E=\left \langle \frac {1}{2}\widehat {\nu _{i}}^*\widehat {\nu _{i}}\right \rangle$, where 
$\nu_i = \sqrt {\rho }\,{u''_{i}}$, 
$\hat {({\cdot })}$ is the one-dimensional Fourier transform in the homogeneous direction, and 
$({\cdot })^*$ is the complex conjugate. The non-dimensional time 
$\tau$ depends on 
$t_0$, which itself depends on initial tilt angle 
$\alpha _0$. Simulation data files may not capture 
$t_0$ exactly for each angle, resulting in slightly different values of 
$\tau$ for each solution file set, so spectra for different initial column angles are compared at close, but non-equal times.
Figure 16. Mean turbulent kinetic energy spectra in the vertical plane at different times for columns with initial angles (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$.
At every angle, an apparent power-law scaling may be observed between wavenumbers 
$\kappa =30$ and 
$\kappa =100$ (indicated by the arrows in figure 16). These wavenumbers correspond to length scales 
$0.4$ mm and 
$0.12$ mm, respectively. The values for the power-law exponent 
$p$ were computed by fitting a line through the apparent power-law region of the energy spectra. These values are presented in table 3. These power-law exponents are all close to 
$-1$ at early times, but decrease slightly in magnitude over time for all considered angles. These values are in agreement with the power-law scalings observed by Olmstead et al. (Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a) for scalar structure functions. In ideal Kolmogorov turbulence, exponents of the energy spectrum and the spectrum of a diffusive scalar (passive or even reactive) should be equal, provided that the same conditions are met as those necessary for the equivalence between spectral and structure function representation of the flow (Monin & Yaglom Reference Monin and Yaglom2013).
Table 3. Energy spectra power-law exponents at different times for different initial column angles.
3.5. Structure functions
To examine further the power-law behaviour of this flow, second-order, scalar structure functions were assembled for the SF
$_6$ mass fraction in a manner similar to those presented in Olmstead et al. (Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a):
\begin{equation} {Y_s}_2(r)=\left\langle[Y_s(\boldsymbol{x}) - Y_s(\boldsymbol{x}+\boldsymbol{r})]^2\right\rangle, \end{equation}
where 
$r=|\boldsymbol {r}|$ is the magnitude of radius vector 
$\boldsymbol {r}$, and 
$\langle {\cdot } \rangle$ denotes averaging over all pairs of points in the same region as the statistics defined previously. The structure function computations have quadratic efficiency, 
$O(n^2)$, and can be prohibitively costly for high-resolution datasets. An algorithm was developed for the current study to produce structure function data that utilized the Kokkos framework (Edwards, Trott & Sunderland Reference Edwards, Trott and Sunderland2014) to take advantage of the performance benefits of GPU systems.
Structure functions evolving in time are shown in figure 17 for different initial gas column angles. For reference purposes, two lines are added to each plot: one corresponds to the power-law exponent 
$2/3$ that is expected theoretically in homogeneous isotropic turbulence (Celani et al. Reference Celani, Cencini, Vergassola, Villermaux and Vincenzi2005; Monin & Yaglom Reference Monin and Yaglom2013), and the other corresponds to the power law with the exponent close to unity observed in some experimental runs of Olmstead et al. (Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a). As the figure demonstrates, numerical results resemble a power law with the exponent close to unity only at the inclination angle 
$\alpha _0=30^\circ$ at early times, and for large scales only. There is no single power-law exponent fitting all scales at any given inclination angle at any time. Overall, the results are remarkably similar to the experiment (figure 3 of Olmstead et al. Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a): the plots can be interpreted as showing the transfer of energy from large scales (here, the integral scale would be related to the geometry of the initial conditions – representative gas column size after compression) to small scales, but an inertial range spanning multiple orders of magnitude does not develop.
Figure 17. Structure functions of mass fraction at various initial inclination angles: (a) 
$\alpha _0=1^\circ$, (b) 
$\alpha _0=5^\circ$, (c) 
$\alpha _0=10^\circ$ and (d) 
$\alpha _0=30^\circ$.
This is not surprising as, strictly speaking, this flow is transitioning to turbulence (at least in the phenomenological sense) in a non-canonical case that is distinct from the homogeneous, isotropic, fully developed classical Kolmogorov turbulence.
First, the entire energy input driving the turbulent transition is provided by the shock interaction with density interfaces over a finite time. Accordingly, the spectral properties of kinetic energy fluctuations in RMI-driven mixing differ from those of Kolmogorov turbulence (Abarzhi & Sreenivasan Reference Abarzhi and Sreenivasan2022).
Second, it is important how the flow is driven on the injection scale. For the flows considered here, a reasonable assumption is that the driver is shear on several scales (dictated by the initial conditions geometry) due to RMI. A recent theory of shear-driven turbulence suggests that the second-order velocity structure function is affected by both the isotropic (
$r^{2/3}$) and non-isotropic (
$r^{4/3}$) contributions (Kumar, Meneveau & Eyink Reference Kumar, Meneveau and Eyink2022), whose respective roles may also be scale- and time-dependent. This could explain some of the behaviours observed in figure 17. A different way to look at decaying turbulence involves considering how the effective length scale varies with time, producing an unsteady cascade (Goto & Vassilicos Reference Goto and Vassilicos2016). A somewhat similar approach in an earlier theory (George Reference George1992) considers the possibility that a universal scaling for spectra may exist not for the wavenumber 
$\kappa$, but for a different object (
$\kappa l$, where 
$l$ is the time-dependent length scale, in the George theory). The best theoretical approach to the flow considered in this work could be one taking into consideration the influence of the initial conditions (irrelevant for the fully developed turbulence case) on scalings. Whether the turbulence is affected by how it is driven in time, in space, or both, the appearance of the spectra (in terms of 
$\kappa$) or structure functions (
$r$) may change.
Third, there is the possibility of a mismatch between numerics and experiment due to diagnostics. In experimental PLIF results, imaging shows the laser-induced fluorescence intensity from the tracer gas that usually closely follows mass concentration; however, possibilities for several diagnostic-related discrepancies remain (Melton & Lipp Reference Melton and Lipp2003). An example of such a discrepancy is discussed in § 3.6. Some subtler differences between experiment and numerics could also be attributable to differences in initial conditions: in the experiments, the diffuse gas column always has small perturbations that the numerical initial conditions do not attempt to reproduce. In a sense, this would mean that in the experiments, there is an additional initial contribution to disorder in the flow.
For possible future studies, it might be worthwhile to consider the approach described by Vorobieff et al. (Reference Vorobieff, Mohamed, Tomkins, Goodenough, Marr-Lyon and Benjamin2003), where ensembles of velocity fields from multiple experimental realizations were collected for the same experiment with nominally identical initial conditions. This made it possible to produce ensemble-averaged velocities and then extract velocity fluctuation fields for each individual experiment. Notably, while the structure functions for the full velocity fields manifested steeper slopes, the fluctuation field slopes were closer to 
$2/3$. In numerics, small fluctuations could be added to the initial conditions to produce similar ensembles (Palekar et al. Reference Palekar, Vorobieff and Truman2007).
3.6. Probability density functions
In this subsection, we present probability density functions (p.d.f.s) for mass fraction, vorticity components and velocity components in order to further characterize the flow development.
To construct discrete p.d.f.s, a bin size is defined to be 
$\Delta \phi =(\phi _{max}-\phi _{min})/N_b$ for some quantity 
$\phi$, where 
$N_b$ is the number of bins. A bin count 
$N_b=96$ was selected to produce smooth p.d.f. profiles with an adequate number of samples within each bin. Values of 
$\phi$ are then distributed into the 
$N_b$ bins, yielding the sample count of 
$N_k$ per bin. The p.d.f. is then defined by 
$P(\phi )={N_k}/{\Delta \phi N}$, where 
$N$ is the total number of samples of 
$\phi$. This definition results in 
$\sum _{k=1}^{N_k}P(\phi )\,\Delta \phi =1$.
The p.d.f.s of the 
${\rm SF}_6$ mass fraction 
$Y_s$ (figure 18) were constructed using the mass fraction data from the same region as the statistics discussed in previous subsections. That is, the vertical cross-section is cropped vertically to retain the central one-third of the column using 
$Y_s=0.001$ as a criterion to determine the upstream and downstream edges of the column. Within this region, 
$Y_s$ falls between 0.0 and 0.27 at all times for each initial tilt angle.
Figure 18. P.d.f.s of SF
$_6$ mass fraction 
$Y_s$ for columns with different initial tilt angles 
$\alpha _0$, at different times 
$\tau$: 5.7 (solid line), 35.7 (dashed line), and 110.8 (dotted line). Angles are (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$.
At early times, immediately after passage of the shock, lower values of the mass fraction (
$\approx 0.0\unicode{x2013}0.15$) are approximately equally likely at all tilt angles of the gas column, while the probability of higher mass fractions decreases rapidly (figures 18a–d). At later times, the p.d.f.s take on a bell-like shape. Tables 4(a)–4(d) present the mean, standard deviation (SD), skewness and kurtosis over time for each initial column angle. These values were computed using the mean(), std(), skewness() and kurtosis() subroutines from the pandas Python library (pandas development team 2020).
Table 4. Mass fraction p.d.f. statistics over time for different initial column angles.
With time, the distribution becomes more narrow at all considered column angles. The mean, corresponding to the peak of the distribution, decreases with time for all initial column angles and approaches nearly the same value (
${\approx }0.063$) for each column angle. That is, at later times, most gas exists in a lower concentration independent of the initial gas column angle. Similar tendencies are observed for the standard deviation.
More variation is observed with time and at different column angles in the skewness and kurtosis values (table 4). As time progresses, the skewness tends towards negative values, and more so as the column angle increases. Somewhat similar dynamics are observed for the kurtosis.
Figure 19 shows p.d.f.s of the 
$z$-component of vorticity that were constructed in the same way as the p.d.f.s of mass fraction. As the figure demonstrates, the column angle has little effect on the p.d.f. of this parameter and how its shape changes with time. At any time, the mean vorticity value is located around zero, with the p.d.f. shape being symmetric. With time, the distribution widens. However, at the largest column angle, 
$30^\circ$, the range of vorticity values does not change much between 
$\tau \approx 35.7$ and 
$\tau \approx 110.8$, which explains the larger mean value at 
$\tau \approx 110.8$ to compare with those at the other angles at the same time.
Figure 19. P.d.f.s of vorticity 
$\omega _z$ for columns with different initial tilt angles 
$\alpha _0$, at different times 
$\tau$: 5.7 (solid line), 35.7 (dashed line), and 110.8 (dotted line). Angles are (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$.
Plots of p.d.f.s of the 
$x$-component of the gas velocity are presented in figure 20. At all considered gas column angles, a bimodal distribution emerges over time with two peaks: an upper peak corresponding to 
$u/u_{post}=1$, and a lower peak between 
$0.90$ and 
$0.95$. The positions of the peaks are not influenced significantly by the column angle, and have similar behaviour over time for each initial column angle. However, their amplitudes vary with time and depend on the column angle. In particular, smaller column angles have a higher initial concentration of velocity around the lower peak, but over time, the higher column angles have a larger concentration of velocity around the lower peak. The opposite is observed for the upper peak.
Figure 20. P.d.f.s of the 
$x$-component of velocity normalized by the post-shock velocity 
$u_{post}$ at different times 
$\tau$: 5.7 (solid line), 35.7 (dashed line), and 110.8 (dotted line). Angles are (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$.
P.d.f.s of the 
$y$-component of the gas velocity are shown in figure 21. At the smallest column angle, there is a single peak value for 
$y$-velocity around 
$v/u_{post}=0$. With time, the distribution widens slightly, with gas starting to move in the vertical direction.
Figure 21. P.d.f.s of the 
$y$-component of velocity normalized by the post shock velocity 
$u_{post}$ at different times 
$\tau$: 5.7 (solid line), 35.7 (dashed line), and 110.8 (dotted line). Angles are (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$.
As the column angle increases, the initial distribution widens. Over time, larger initial column angles develop a bimodal distribution of velocity in the 
$y$-direction, with a positive and negative mean. At 
$\alpha \geqslant 10^\circ$, the second peak becomes obvious with time. Both peaks have similar values at 
$\tau =21.4$, and their locations are shifted from 
$0\ {\rm m}\ {\rm s}^{-1}$ to negative and positive values. At the latest time, 
$\tau =110.8$, only one peak remains in the p.d.f., corresponding to 
$\alpha = 10^\circ$, with the peak value being positive (
$v/u_{post}\approx 0.003$). This again indicates gas moving upwards. At the largest column angle, 
$\alpha = 30^\circ$, both peaks are seen at 
$\tau =110.8$ (figure 21d), but the peak with the positive value is larger (
$v/u_{post}\approx 0.013$), showing that more gas is moving upwards than downwards.
The component of gas velocity in the 
$z$-direction (shown in figure 22) maintains a narrow distribution around 
$w/u_{post}=0$ for all initial column angles. Over time, the distribution widens slightly for smaller initial column angles. The zero mean is due to the symmetry of the sample plane, as discussed previously.
Figure 22. P.d.f.s of the 
$z$-component of velocity normalized by the post shock velocity 
$u_{post}$ at different times 
$\tau$: 5.7 (solid line), 35.7 (dashed line), and 110.8 (dotted line). Angles are (a) 
$1^\circ$, (b) 
$5^\circ$, (c) 
$10^\circ$ and (d) 
$30^\circ$.
It is informative to compare the p.d.f.s of the SF
$_6$ mass fraction (figure 18) with the experimental p.d.f.s of PLIF intensity in a very similar flow (figure 16 of Olmstead et al. Reference Olmstead, Wayne, Simons, Monje, Yoo, Kumar, Truman and Vorobieff2017a). While the relationships of the overall trends for the mixing-driven p.d.f. time evolution in the flow are clear, the PLIF intensity p.d.f.s are dominated by a peak corresponding to the tracer-free (no fluorescence) flow, which is understandably absent from figure 18.
4. Conclusion
High-resolution viscous simulations were conducted to model the 3-D interaction between a planar shock and an inclined cylindrical column of heavy gas with an initially diffuse density interface. The overall flow morphology was similar to the inviscid case presented in our previous study (Romero et al. Reference Romero, Poroseva, Vorobieff and Reisner2022). However, the increased resolution allowed for more detailed statistics to be collected. The effects of initial column angle on statistical properties of the flow were analysed.
In particular, the analysis of mean velocity profiles showed that the streamwise velocity component is at the maximum near the column edges at all times. The magnitude of the vertical component initially has a maximum in the interior of the column, but approaches zero everywhere over time. The spanwise component is close to zero at all times everywhere due to the symmetry of the cross-sectional plane.
The intensities of all fluctuating velocities and of the shear stresses initially have a maximum near the centre (interior) of the column. Over time, the peak values of all Reynolds stresses move towards the leading edge of the column. The trailing edge of the column does not experience significant velocity fluctuations at any time.
Analysis of the anisotropy tensor shows that turbulent kinetic energy is non-uniformly distributed through the column. At larger column angles, turbulence tends towards isotropy in the column interior in the 
$y$- and 
$z$-components mainly. Anisotropy-invariant maps show that the turbulent state is near to that of single-component turbulence at early times, and that the flow moves in the direction of isotropy over time.
Analysis of the mixedness over time showed that the initial angle has very little effect on the mixing rate.
A power-law scaling with the exponent close to 
$-1$ was observed in the turbulent kinetic energy spectra at large wavenumbers. This scaling is similar to the power-law scaling for scalar structure functions obtained in some experiments. However, the scalar structure functions for the mass fraction based on the simulation data do not follow the power law except for large initial column angles at early times. The scalar field (mass fraction) structure function behaviour in our simulation is overall quite similar to that observed in experiment (tracer fluorescence intensity fields).
The mean and the standard deviation of the mass fraction p.d.f.s were found to be independent of the column angle, while both the skewness and the kurtosis generally decrease as the initial column angle increases. Vorticity is distributed symmetrically around 0 and does not have a significant dependence on initial column angle. Both 
$x$- and 
$y$-components of velocity tend towards a bimodal distribution as the initial column tilt angle is increased. The 
$z$-component of velocity maintains a narrow distribution around zero, as expected due to the symmetry of the sample plane.
Funding
This work is supported in part by subcontract 594693 ‘Development of a hydrodynamic code for large-scale atmospheric event simulations on modern and future computational platform’, from Los Alamos National Laboratory to the University of New Mexico, and DTRA grant HDTRA-18-1-0022. A part of the work conducted at the University of New Mexico was also supported from DOE grant DE-NA0004108, ‘The Rio Grande Consortium for Advanced Research on Exascale Simulation (Rio Grande CARES)’. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the US Department of Energy National Nuclear Security Administration under contract no. 89233218CNA000001.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Grid convergence
Here, grid sensitivity analysis results are presented for simulations of the gas column initially inclined at 
$30^{\circ }$ and with Mach number 
$2.0$. The grid resolutions considered are 
${{\rm d}x}=40\ \mathrm {\mu } \mathrm {m}$, 
$80\ \mathrm {\mu } \mathrm {m}$ and 
$160\ \mathrm {\mu } \mathrm {m}$ – that is, two additional grid resolutions, which are coarser than the grid used for the results presented in the main text. Simulations on finer grids could not be conducted due to insufficient computational resources.
The behaviour of the KH wavelengths with respect to grid cell size is shown in figure 23. As described previously, the KH wavelength 
$\lambda _{KH}$ was determined by visually identifying wave peaks when they first become visible after shock passage, and measuring the distance between them. For each case, at least ten consecutive wave peaks were identified. The figure shows the average wavelength, along with error bars indicating one standard deviation. As the plot shows, the two finest grids considered, 
$40\ \mathrm {\mu } \mathrm {m}$ and 
$80\ \mathrm {\mu } \mathrm {m}$, agree with regard to 
$\lambda _{KH}$, while the 
$160\ \mathrm {\mu }{\rm m}$ grid overpredicts wavelengths by a factor of about two.
Figure 23. KH wavelengths 
$\lambda _{KH}$ at different grid sizes for a column with initial tilt angle 
$30^{\circ }$. Error bars indicate one standard deviation.
Effects of the grid cell size on the molecular mixedness were also considered and are presented in figure 24. As the figure demonstrates, the coarser grids predict a higher mixing rate than the finest grid at early times. This may be due to the averaging effect of the larger grid cells, which are unable to resolve fine-scale concentration fluctuations. At later times, however, the grid size has no distinguishable effect on the mixing efficiency.
Figure 24. Molecular mixedness over time for different grid resolutions.
P.d.f.s of 
${\rm SF}_6$ mass fraction are shown in figure 25 for the grids at two different times. At early times (
$\tau = 5.7$), there is a higher probability of SF
$_6$ mass fraction in the range from 0.05 to 0.2 for the 
${{\rm d}x}=40\ \mathrm {\mu } {\rm m}$ grid when compared to the coarser grids. This peak at higher mass fractions results in the lower mixedness of the fine grid at early times, as seen in figure 24. At 
$\tau =35.7$ (figure 25b), SF
$_6$ distributions for 
${{\rm d}x}=40\ \mathrm {\mu } {\rm m}$ and 
$80\ \mathrm {\mu } {\rm m}$ are very similar, while the mass fraction distribution for the coarsest grid, though it has the same mean, is much wider. Nonetheless, the coarse grid shows identical mixedness to the finer grids despite the larger variance. The SF
$_6$ mass fraction distribution converges faster with respect to grid size at late times than it does at early times, where it is more sensitive to grid size.
Figure 25. P.d.f.s of 
${\rm SF}_6$ mass fraction 
$Y_s$ for different grid resolutions at (a) 
$\tau =5.7$ and (b) 
$\tau =110.8$.
A similar phenomenon is observed for energy spectra. At early times, in this case 
$\tau =35.7$ (figure 26a) after the flow has become turbulent, the coarse grids underestimate the spectral energy when compared to the finest grid. At later times (
$\tau =65.7$, figure 26b), the spectra for the two finest grids are in better agreement, while the coarse grid severely underpredicts spectral energy by an order of magnitude.
Figure 26. Mean turbulent kinetic energy spectra in the vertical plane for different grid resolutions at (a) 
$\tau =35.7$ and (b) 
$\tau =65.7$.
In general, the flow geometry, in terms of 
$\lambda _{KH}$, is well represented by the 
${{\rm d}x}=40\ \mathrm {\mu } {\rm m}$ grid, and late-time mixedness and mass fraction p.d.f.s are converged. Early-time statistics, however, are more sensitive to grid resolution.
Open access
