Density effects on post-shock turbulence structure and dynamics

Abstract

Turbulence structure resulting from multi-fluid or multi-species, variable-density isotropic turbulence interaction with a Mach 2 shock is studied using turbulence-resolving shock-capturing simulations and Eulerian (grid) and Lagrangian (particle) methods. The complex roles that density plays in the modification of turbulence by the shock wave are identified. Statistical analyses of the velocity gradient tensor (VGT) show that density variations significantly change the turbulence structure and flow topology. Specifically, a stronger symmetrization of the joint probability density function (PDF) of second and third invariants of the anisotropic VGT, PDF$(Q^{\ast },R^{\ast })$, as well as the PDF of the vortex stretching contribution to the enstrophy equation, are observed in the multi-species case. Furthermore, subsequent to the interaction with the shock, turbulent statistics also acquire a differential distribution in regions having different densities. This results in a nearly symmetric PDF$(Q^{\ast },R^{\ast })$ in heavy-fluid regions, while the light-fluid regions retain the characteristic tear-drop shape. To understand this behaviour and the return to ‘standard’ turbulence structure as the flow evolves away from the shock, Lagrangian dynamics of the VGT and its invariants is studied by considering particle residence times and conditional particle variables in different flow regions. The pressure Hessian contributions to the VGT invariants transport equations are shown to be not only affected by the shock wave, but also by the density in the multi-fluid case, making them critically important to the flow dynamics and turbulence structure.

1 Introduction

The interaction of a normal shock wave with multi-fluid or multi-species isotropic turbulence is an extension of the canonical shock–turbulence interaction (STI) problem which includes strong variable-density effects. This extended configuration can enhance our understanding of more complex flow problems such as fuel–air mixing in supersonic combustion, the interaction of supernova remnants with interstellar clouds, shock propagation through foams and bubbly liquids, inertial confinement fusion and re-shock problem in Richtmyer–Meshkov instability. Most of the previous theoretical, numerical and experimental studies of STI have been dedicated to the original canonical problem.

The early theoretical study by Ribner (1954) restricted the STI to the linear interaction regime with a large-scale separation between the shock and turbulence, so that the nonlinear and viscous effects are assumed to be negligible during the interaction. By decomposing the pre-shock turbulence into independent modes (acoustic, vortical and entropy) using Kovasznay decomposition (Kovasznay 1953), the post-shock turbulence statistics can be theoretically derived from the linearized Rankine–Hugoniot jump conditions. This approach is referred to as the linear interaction approximation (LIA) and represents an important limiting case, since it provides analytical predictions for the jumps of fluctuating quantities across the shock.

Due to the challenges of accurate experimental measurements of the smallest time and length scales around the shock wave, numerical simulations have been widely employed to investigate this interaction. Researchers have used both shock-capturing and shock-resolving simulations to understand the post-shock amplification of Reynolds stress, vorticity variance and turbulent length scales (Lee, Lele & Moin 1993; Hannappel & Friedrich 1995; Mahesh et al. 1995; Lee, Lele & Moin 1997; Mahesh, Lele & Moin 1997; Jamme et al. 2002; Larsson & Lele 2009; Larsson, Bermejo-Moreno & Lele 2013). Earlier numerical studies showed limited agreement with the LIA predictions because the parameter range was outside the linear regime. More recently, Ryu & Livescu (2014) have considered a wide range of parameters in their shock-resolving direct numerical simulations (DNS) to show that the DNS results converge to the LIA solutions when the ratio of the shock thickness ( $\unicode[STIX]{x1D6FF}$ ) to the pre-shock Kolmogorov length scale ( $\unicode[STIX]{x1D702}$ ) becomes small. Replacing the actual shock interaction with the LIA relations can extend the reach of DNS to arbitrarily high shock Mach numbers and much larger Taylor Reynolds number ( $Re_{\unicode[STIX]{x1D706}}$ ) than otherwise computationally feasible, provided that the interaction parameters correspond to the linear regime. This method (named Shock-LIA by the authors) was used for detailed studies of post-shock turbulent energy flux and vorticity dynamics (Livescu & Ryu 2016; Quadros, Sinha & Larsson 2016). Sethuraman, Sinha & Larsson (2018) used shock-capturing simulation and LIA to study the thermodynamic field generated by STI. In a recent study (Tian et al. 2017a ), we showed, using shock-capturing turbulence-resolving simulations, that the LIA predictions for the Reynolds stresses can be approached provided that the scale separation between numerical shock thickness ( $\unicode[STIX]{x1D6FF}_{n}$ ) and Kolmogorov length scale is sufficient. Thus, when the ratio of turbulent to shock scales is large enough, so that the numerical artifacts near the shock do not influence the flow, the shock-capturing method can correctly simulate the STI.

As mentioned above, in many practical applications, STI may occur in a mixture of fluids of very different densities. This motivated our extension of the canonical STI problem to include variable-density effects (Tian et al. 2017a , c ) by considering the pre-shock turbulence as an isotropic mixture of two fluids (species) with different molecular weights, as encountered in non-premixed combustion. Using turbulence-resolving shock-capturing simulations, we have examined the turbulence statistics, turbulence budgets, conditional statistics and energy spectrum in multi-fluid STI and found that the nonlinear effects from the density variations significantly change the turbulence properties in both physical and spectral spaces. The relation between velocity and a passive scalar field has also been studied by Buttay, Lehnasch & Mura (2016) and Boukharfane, Bouali & Mura (2018). Other studies (Jin et al. 2015; Huete et al. 2017) used LIA and shock-capturing simulations to investigate the interaction of a reactive premixed mixture with shock and turbulence. These studies help in better understanding of complex STI problem. However, there still exist many gaps in our knowledge of the variable-density effects on the post-shock turbulence structure and flow topology.

In this study, we focus on the density effects on the post-shock turbulence structure by examining the velocity field. The properties of the velocity gradient tensor (VGT) determine a wide variety of turbulence characteristics, such as the flow topology, deformation of material volume, energy cascade and intermittency. Understanding both the VGT field immediately after the shock wave and its dynamics as the flow evolves away from the shock wave is also crucial to the development of subgrid-scale models that can accurately describe the shock interaction and return-to-isotropy effects. Perry & Chong (1987) and Chong, Perry & Cantwell (1990) proposed an approach to classify the local flow topology and structure using the invariants of the VGT. The dynamical behaviour of the VGT has been studied for incompressible flows using the Lagrangian evolution of the invariants along conditional mean trajectories (Meneveau 2011). The statistics regarding the invariants of the VGT and their Lagrangian dynamics have been used to understand the structure of turbulence in many canonical flows, such as isotropic turbulence, turbulent boundary layers and mixing layers (e.g. Chong et al. 1998; Ooi et al. 1999; Wang & Lu 2012; Bechlars & Sandberg 2017). Previous studies on single-fluid STI have examined the probability density function (PDF) of the VGT. Ryu & Livescu (2014) and Livescu & Ryu (2016) took a step further to investigate the turbulence structure and vorticity dynamics based on the examination of VGT invariants. By taking advantage of the Shock-LIA method, they extracted the statistics of the VGT and its invariants for a wide range of shock Mach numbers, even though the dynamics of the VGT as the turbulence evolves away from the shock wave could not be examined with the Shock-LIA method. Our earlier numerical studies of variable-density STI have revealed some important new features of velocity and scalar statistics in this set-up (Tian et al. 2017b ; Tian, Jaberi & Livescu 2018). However, these studies have not yet fully identified the variable-density effects on the post-shock turbulence/scalar structure.

This study uses the recently generated database of turbulence-resolving shock-capturing simulations of multi- and single-fluid STI to: (1) develop a better understanding of variable-density and shock effects on the turbulence structure immediately after the shock wave and (2) perform the first Lagrangian analysis of this flow configuration for better understanding of the dynamical behaviour of the VGT as the turbulence evolves away from the shock. While the compressibility effects are weak for the current parameter range and not discussed, variable-density effects are very significant and the focus of this study. The paper is organized as follows. Details of the simulations and the testing conducted to assess the accuracy of the Lagrangian and Eulerian analysis are discussed in § 2. Results are presented in § 3 and concluding remarks are made in § 4.

2 Numerical method and accuracy

In this section, we first briefly discuss the numerical approach used for shock-capturing turbulence-resolving simulations in our previous study (Tian et al. 2017a ), from which we have extracted the VGT statistics addressed in this paper. The extended variable-density STI configuration is described next, followed by a discussion of the new Lagrangian simulations used to examine the VGT dynamics away from the shock.

2.1 Governing equations and numerical approach

The conservative form of the dimensionless compressible Navier–Stokes equations for flows with two miscible species (i.e. continuity, momentum, energy and species mass fraction transport equations) has been solved numerically together with the perfect gas law using a high-order hybrid numerical method (Tian et al. 2017a ). The inviscid fluxes for the transport equations have been computed using the fifth-order monotonicity-preserving scheme, as described in Li & Jaberi (2012). The molecular transport terms have been calculated using the sixth-order compact scheme (Lele 1992). The third-order Runge–Kutta scheme has been used for time advancement.

Figure 1. Instantaneous contours of vorticity and shock surface in isotropic turbulence interacting with a Mach 2 shock. (a) Vortex structures are identified by the $Q$ criterion (i.e. iso-surface of the second invariant of VGT: $Q=2\langle Q_{w}\rangle$ , where $\langle Q_{w}\rangle$ is the averaged magnitude of the rotation tensor), coloured by the mole fraction of the heavy fluid. Fluid particles are initialized as a sheet that spans over the homogeneous directions at a given post-shock streamwise position and allowed to develop with the flow. (b) Visualized particle sheet, convected and distorted by the post-shock turbulence.

2.2 Numerical set-up

The physical domain for the simulations considered in this paper is a box that has a dimension of 4 $\unicode[STIX]{x03C0}$ in the streamwise direction (denoted as $x$ ) and $(2\unicode[STIX]{x03C0},2\unicode[STIX]{x03C0})$ in the transverse directions (denoted as $y$ and $z$ ), as shown in figure 1(a). The flow in this figure is visualized using the iso-surface of $Q$ , the second invariant of the VGT $\unicode[STIX]{x1D608}_{ij}=\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j}$ . The normal shock is located at $x=2\unicode[STIX]{x03C0}$ . A buffer layer is used at the end of the computational domain from $4\unicode[STIX]{x03C0}$ to $6\unicode[STIX]{x03C0}$ to eliminate reflecting waves. In the transverse directions, periodic boundary conditions are used as the flow is assumed to be periodic and homogeneous in these directions. To provide inflow turbulence, pre-generated decaying isotropic turbulence is superposed on the uniform mean flow with Mach number of 2.0 and convected into the domain using Taylor’s hypothesis. The inflow turbulent Mach number, Reynolds number and peak wavenumber are $M_{t}\approx 0.1$ , $Re_{\unicode[STIX]{x1D706}}\approx 45$ and $k_{0}=4$ , respectively. For this $M_{t}$ value, Taylor’s hypothesis is appropriate for approximating spatially developing turbulence with temporally developing turbulence (Lee, Lele & Moin 1992). The variable-density (multi-fluid) effects arise from compositional variations of a binary mixture of miscible fluids with different molar masses, which is generated by correlating the density to an isotropic scalar field representing the mole fraction of the heavy fluid. The scalar field is generated as a random field following a Gaussian spectrum with a peak at $k_{s}=8.0$ and has double-delta PDF distribution so that the scalar value initially is either 1.0 or 0. The initial scalar field is smoothed by solving a diffusion equation so that the scalar field can be fully resolved by the chosen mesh. The resulting scalar field is then allowed to decay in the fully developed isotropic turbulence set-up for one eddy turnover time as a passive scalar. The density field is then calculated by imposing $X=\unicode[STIX]{x1D719}$ (where $X$ is the mole fraction of the heavy fluid). The generated variable-density isotropic turbulence is then superposed onto the mean flow and allowed to develop into a more realistic state before reaching the shock wave. The Atwood number, $A_{t}=(MW_{2}-MW_{1})/(MW_{2}+MW_{1})$ , calculated from the molar masses of the two fluids, $MW_{1}$ and $MW_{2}$ , is 0.28. This value of the Atwood number was chosen such that the variable-density effects are non-negligible, yet the interaction with the shock wave is still in the wrinkled-shock regime. At larger Atwood numbers, the interaction enters the broken-shock regime, where more complicated dynamics exists. The extension of the current study to this regime poses significant challenges, which are beyond the goals of the current study. The Prandtl number, $Pr$ , and Schmidt number, $Sc$ , are the same and equal to 0.75. Immediately before the shock wave, $M_{t}$ and $Re_{\unicode[STIX]{x1D706}}$ reach around 0.09 and 42 due to turbulence decay. For these values, the nonlinear and viscous effects on turbulence passing through the shock wave are weak based on the results of LIA convergence tests done in our previous study Tian et al. (2017a ).

2.3 Numerical method for the Lagrangian study

For the current study, we have tracked more than 4.5 million particles that are initialized uniformly at various streamwise positions $\boldsymbol{x}_{\mathbf{0}}$ , and calculated various turbulence statistics following their trajectories. The aim is to understand the evolution of flow structures following fluid particles as the turbulence develops downstream of the shock. Figure 1(a) marks with red lines a typical streamwise plane where particles are initialized. Each set of particles is initialized uniformly in the spanwise directions at the same streamwise location, corresponding to a planar sheet. The spacing between the neighbouring particles in the spanwise directions is the same as the grid size $(2\unicode[STIX]{x03C0}/512)$ . We have uniformly sampled around 20 particle sets (sheets) for each cycle of the inflow turbulence box. The particles are then convected by the instantaneous turbulent velocity obtained by turbulence-resolving shock-capturing simulations and moved to a region marked by the blue lines. At this stage, the initially flat particle sheet is distorted by the turbulence as shown in figure 1(b).

The fluid particles are non-inertial and follow the local flow velocity. The corresponding transport equations for particle positions $x_{i}^{+}$ are

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}x_{i}^{+}(t\mid \boldsymbol{x}_{\boldsymbol{ 0}},t_{0})}{\text{d}t}=u_{i}^{+}(t\mid \boldsymbol{x}_{\mathbf{0}},t_{0}), & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle u_{i}^{+}(t\mid \boldsymbol{x}_{\mathbf{0}},t_{0})=u_{i}(x_{i}^{+},t), & \displaystyle\end{eqnarray}$$

where $x_{i}^{+}(t\mid \boldsymbol{x}_{\mathbf{0}},t_{0})$ represents the positions of the particles at time $t$ that are initialized at $\boldsymbol{x}_{\mathbf{0}}$ and time $t_{0}$ . The particle velocity $u_{i}^{+}(t\mid \boldsymbol{x}_{\mathbf{0}},t_{0})$ can be obtained from the Eulerian velocity field $u_{i}(x_{i}^{+},t)$ by interpolation. The interpolation is based on the cubic spline scheme, whose accuracy in predicting particle positions has been studied in Yeung & Pope (1988). The time-stepping scheme for Lagrangian particles is also the third-order Runge–Kutta scheme. Therefore, at each sub-time step, the particle velocity is interpolated from the Eulerian velocity field with the same sub-time step. In the STI configuration, there is a sharp change of the flow velocity at the shock, which reduces the interpolation accuracy. To achieve accurate interpolation of the particle velocity, the domain is partitioned into three different regions as shown in figure 1(a): pre-shock, shock and post-shock regions. The instantaneous shock surface is identified using the sensor: $s=-\unicode[STIX]{x1D703}/(|\unicode[STIX]{x1D703}|+\langle \unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{i}\rangle _{yz}^{0.5})>0.5$ (Larsson & Lele 2009), where $\unicode[STIX]{x1D703}=\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{i}$ is the dilatation, $\unicode[STIX]{x1D714}_{i}=\unicode[STIX]{x1D716}_{ijk}\unicode[STIX]{x2202}u_{k}/\unicode[STIX]{x2202}x_{j}$ is the vorticity and $\langle \rangle _{yz}$ represents the instantaneous average over the homogeneous directions. After the instantaneous shock region is identified, the pre- and post-shock turbulence fields can be separated for interpolation. Note that the cubic spline interpolation scheme requires information from neighbouring cells, so a buffer region (around three grid points) is added between the shock region and the post-shock region. Lagrangian dynamics of particles across the shock wave is not considered in this study because the shock profile is numerical and its thickness depends on the grid size. This introduces numerical artifacts when considering the particle dynamics across the shock wave.

2.4 Grid and statistical convergence

The accuracy of the numerical results is addressed in this subsection through a series of convergence tests. To ensure that all the turbulence length scales are well resolved, a grid convergence test was conducted in Tian et al. (2017a ). Here, we summarize these results for completeness, together with additional convergence results for small-scale quantities. Figure 2 shows the turbulence dissipation rate $\unicode[STIX]{x1D700}=-\overline{\unicode[STIX]{x1D70E}_{ij}(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})}$ , where $\unicode[STIX]{x1D70E}_{ij}$ is the viscous stress tensor, and scalar (mass fraction for the multi-fluid STI) dissipation rate $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D719}}=\overline{(\unicode[STIX]{x1D707}/(Re_{0}Sc))(\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}x_{j})(\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}x_{j})}$ as a function of the normalized streamwise direction $k_{0}x$ for a series of meshes. The grey regions in the following figures indicate the unsteady shock region, inside which the results are affected by the shock wrinkling and unsteady shock movement. As the grid is refined in all three directions, both quantities display convergence, proving the accuracy of the turbulence database. Another issue that needs to be considered is the scale separation between the numerical shock thickness $\unicode[STIX]{x1D6FF}_{n}$ and the Kolmogorov length scale $\unicode[STIX]{x1D702}$ as suggested in our previous study (Tian et al. 2017a ). Thickness $\unicode[STIX]{x1D6FF}_{n}$ is calculated as $(u_{1,u}-u_{1,d})/|\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}|_{max}$ , where $|\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{1}|_{max}$ denotes the maximum magnitude of streamwise velocity gradient. Grid numbers for grids 1 to 5 shown in figure 2 are $256\times 256\times 1024$ , $384\times 384\times 1024$ , $384\times 384\times 1536$ , $512\times 512\times 1536$ and $512\times 512\times 2048$ . With the finest mesh $(512\times 512\times 2048)$ , the scale separation ratio $\unicode[STIX]{x1D702}/\unicode[STIX]{x1D6FF}_{n}$ is around 1.9, which is sufficient for resolving the interaction between the numerical shock wave and small-scale turbulent motions. Therefore, in the current study, we have obtained all the statistics from the turbulence field based on the finest grid to ensure accuracy. Finally, LIA convergence tests were conducted in Tian et al. (2017a ) following Ryu & Livescu (2014) to show that the shock-capturing simulations can capture the correct limits. Turbulent Mach number $M_{t}$ and Taylor Reynolds number $Re_{\unicode[STIX]{x1D706}}$ were varied for the canonical single-fluid simulations, covering a wide range of parameter space. The shock-capturing simulation results do converge to LIA predictions for individual Reynolds stress components as long as certain conditions are satisfied (Tian et al. 2017a ). This was the first time that the asymptotic values for individual Reynolds stresses were approximated using shock-capturing simulations.

Statistical convergence is another important factor that needs to be examined. To reduce the statistical variability, all the results that are based on the Eulerian data are space-averaged over homogeneous directions and time-averaged for around two flow-through times. The averaging is performed after the flow has reached a statistically steady state to eliminate the effects of transient processes (Larsson et al. 2013). For the Lagrangian statistics, the number of fluid particles needs to be large enough for statistical convergence, especially for conditional averaged statistics. The conditional averaged value of $X$ , conditioned on the variables $A$ and $B$ , is defined as

(2.2) $$\begin{eqnarray}\displaystyle \langle X\mid (A=A_{0},B=B_{0})\rangle & = & \displaystyle \langle \!X\mid (A_{0}-{\textstyle \frac{1}{2}}\unicode[STIX]{x0394}A)\leqslant A<(A_{0}+{\textstyle \frac{1}{2}}\unicode[STIX]{x0394}A),\nonumber\\ \displaystyle & & \displaystyle (B_{0}-{\textstyle \frac{1}{2}}\unicode[STIX]{x0394}B)\leqslant B<(B_{0}+{\textstyle \frac{1}{2}}\unicode[STIX]{x0394}B)\!\rangle ,\end{eqnarray}$$

where $\unicode[STIX]{x0394}A$ and $\unicode[STIX]{x0394}B$ are the bin sizes. The conditional statistics are obtained by ensemble averaging (denoted by $\langle \rangle$ ) over all the fluid particles that fall into the bins. Figures 3 and 4 show the convergence of two important conditional Lagrangian statistics $\langle \text{D}Q/\text{D}t\rangle /\langle Q_{w}\rangle ^{3/2}$ and $\langle \text{D}R/\text{D}t\rangle /\langle Q_{w}\rangle ^{2}$ , and their standard deviation, depending on the number of particles in each bin. Here, $\text{D}Q/\text{D}t$ and $\text{D}R/\text{D}t$ represent the material derivative of the second invariant ( $Q$ ) and third invariant ( $R$ ) of the VGT. For the multi-fluid case, we note that the convergence of both conditional means and standard deviations can be achieved when using around 10 000 particles, larger than that needed for the canonical single-fluid case as shown in figure 4. This suggests that the variable-density effects make the simulations more computationally demanding. The effects of the bin sizes are also examined by comparing three different sets of bin numbers, $30\times 30$ (solid), $40\times 40$ (dashed) and $60\times 60$ (dotted), in the $(Q,R)$ phase plane at the same point $(3.0,3.0)$ . These bin numbers correspond to the following bin sizes: $(1.3,1.3)$ , $(1.0,1.0)$ and $(0.67,0.67)$ . Our analysis indicates that the statistics converge to almost the same values when the sample size is large enough. In the present study, we uniformly sampled more than 4.5 million particles and made sure that there are at least 10 000 particles in each sample bin with the number of bins being $40\times 40$ ( $(\unicode[STIX]{x0394}Q,\unicode[STIX]{x0394}R)=(1.0,1.0)$ ).

Figure 2. Results of multi-fluid grid convergence tests at $Re_{\unicode[STIX]{x1D706}}=45$ and $M_{t}=0.1$ . Streamwise development of (a) turbulent dissipation rate $\unicode[STIX]{x1D700}$ and (b) mass fraction dissipation rate $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D719}}$ . The region of unsteady shock movement is marked in grey.

Figure 3. The statistical convergence for (a)  $(\text{D}Q/\text{D}t)/\langle Q_{w}\rangle ^{3/2}$ and $(\text{D}R/\text{D}t)/\langle Q_{w}\rangle ^{2}$ and (b) their standard deviations conditioned at point $(3.0,3.0)$ in the $(Q,R)$ phase plane for the multi-fluid case.

Figure 4. The statistical convergence for (a)  $(\text{D}Q/\text{D}t)/\langle Q_{w}\rangle ^{3/2}$ and $(\text{D}R/\text{D}t)/\langle Q_{w}\rangle ^{2}$ and (b) their standard deviations conditioned at point $(3.0,3.0)$ in the $(Q,R)$ phase plane for the single-fluid case.

3 Results and discussion

The variable-density effects on the post-shock turbulence structure and dynamics are examined in this section. The results obtained from the multi-fluid STI simulation are compared with those of a reference single-fluid case and standard isotropic turbulence. First, the post-shock turbulence state and its evolution away from the shock wave are examined to identify the variable-density effects. The results are based on time- and space-averaged statistics obtained from the Eulerian data. The flow topology is studied next to further understand the post-shock turbulence evolution. The dynamics that dominates the transient evolution of post-shock turbulence structure is examined using the Lagrangian equation of the VGT and Lagrangian data collected for sample fluid particles.

3.1 Density effects on post-shock turbulence

3.1.1 Turbulence state immediately after the shock

In this subsection, the turbulence structure immediately after the shock wave is analysed to identify the different roles that density plays through the shock wave.

The PDFs of streamwise and spanwise longitudinal velocity derivatives for pre- and post-shock ( $k_{0}x=0.5$ ) multi-fluid turbulence are shown in figure 5(a) alongside the Gaussian distribution as a reference. The non-Gaussian nature of the velocity gradient PDFs and their connection to the energy cascade and intermittency are well documented in the turbulence literature. The PDFs of the pre-shock velocity derivatives are negatively skewed as expected. After passing the shock wave, they become closer to the Gaussian distribution, especially for the streamwise component. The PDFs for both single-fluid and multi-fluid post-shock turbulence are shown in figure 5(b). Here, we note that immediately after the shock wave, the PDF of the spanwise velocity gradient for both cases remains negatively skewed, as in isotropic turbulence. The streamwise component, however, becomes more symmetric and Gaussian-like due to the interaction with the shock wave. This indicates that the energy transfer to small scales is suppressed in the streamwise direction. We also note that the density has a relatively weak effect on the velocity derivative PDFs since the single-fluid and multi-fluid cases have similar PDFs.

Figure 5. Comparison of the PDFs of the normalized post-shock velocity derivatives with a Gaussian distribution. Comparison of (a) multi-fluid pre-shock with post-shock results and (b) multi-fluid with single-fluid post-shock results.

The preferential amplification of the transverse components of the rotation and strain rate tensors is an important effect in STI and has been extensively studied for the canonical single-fluid flows (Mahesh et al. 1997; Ryu & Livescu 2014; Livescu & Ryu 2016). This amplification can lead to an increase in the correlation between the two quantities. To better understand the variable-density effects on post-shock turbulence, the PDF of the strain-enstrophy angle, $\unicode[STIX]{x1D6F9}$ , is considered in figure 6. Angle $\unicode[STIX]{x1D6F9}$ is calculated using $\unicode[STIX]{x1D6F9}=$ tan $^{-1}(\unicode[STIX]{x1D61A}_{ij}\unicode[STIX]{x1D61A}_{ij}/(\unicode[STIX]{x1D61E}_{ij}\unicode[STIX]{x1D61E}_{ij}))$ , where $\unicode[STIX]{x1D61A}_{ij}=1/2(\unicode[STIX]{x1D608}_{ij}+\unicode[STIX]{x1D608}_{ji})$ and $\unicode[STIX]{x1D61E}_{ij}=1/2(\unicode[STIX]{x1D608}_{ij}-\unicode[STIX]{x1D608}_{ji})$ are the strain and rotation tensors. In isotropic turbulence, the PDF of $\unicode[STIX]{x1D6F9}$ peaks near $90^{\circ }$ (Jaberi, Livescu & Madnia 2000), indicating a strain-dominated flow. In single-fluid post-shock turbulence, the PDF of $\unicode[STIX]{x1D6F9}$ exhibits a shift of the peak from $90^{\circ }$ to around $45^{\circ }$ , as the shock Mach number increases. This has been observed by Livescu & Ryu (2016) and is interpreted as the increase in correlation of strain and rotation. However, in the multi-fluid case, the peak still occurs at relatively large angles and the increase in correlation is not as pronounced as that in the single-fluid case, at the same shock Mach number. Figure 6 implies that the rotation and strain are amplified differently by the shock when large density variations are present, which compromises the correlation between the two quantities.

Figure 6. The PDF of the strain-enstrophy angle $\unicode[STIX]{x1D6F9}$ in radians for post-shock turbulence.

Figure 7. Conditional expectation of the magnitude of strain rate tensor as a function of density after the shock wave.

The variable-density effects on strain and rotation tensors can be studied by examining the conditional expectations of their magnitudes as a function of density. It was shown in Tian et al. (2017a ) that through the shock wave, the amplification of vorticity is stronger in the mixed-fluid regions with near-average density, but weaker in the pure-fluid regions. This is not observed in the single-fluid simulation. One mechanism that might be responsible for this behaviour is the baroclinic torque $(\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \unicode[STIX]{x1D735}p)/\unicode[STIX]{x1D70C}^{2}$ in the vorticity transport equation. A strong pressure gradient $\unicode[STIX]{x1D735}p$ exists through the shock wave; at the same time, large density gradients $\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$ also exist, especially in the mixed-fluid regions. Since the pre-shock density field is isotropic, $\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$ and $\unicode[STIX]{x1D735}p$ can be locally misaligned, especially when the spanwise component of $\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$ is large, becoming a source of vorticity generation through the baroclinic torque. In addition, the generated vorticity field should be perpendicular to the spanwise density gradient. In the pure-fluid regions or single-fluid simulation, however, the density gradients are much smaller, so that the cross-product of $\unicode[STIX]{x1D735}p$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$ is also small. Note that the density gradient in the streamwise direction has no contribution, because it is aligned with the pressure gradient. To confirm this, the PDF of the angle between the spanwise component of density gradient and the vorticity vector is plotted in figure 8. After the shock wave, the multi-fluid case exhibits a stronger tendency of the vorticity vector being perpendicular to the density gradient. In contrast, this tendency is not observed in the single-fluid case. This provides evidence that density gradient and baroclinic torque play important roles in establishing the preferential deposition of vorticity across the shock wave.

Figure 9 can help visualize the changes in the flow structure across the shock wave. The vortex tubes are captured using the $Q$ -criterion and are coloured by their local density. Figure 9(a) shows the vortex structures for pre-shock multi-fluid isotropic turbulence. For the visualized vortex tubes, there are no identifiable effects from the density variations; the vortex tubes are not preferentially distributed due to the density effects. However, the interaction with the shock has a clear effect on the post-shock vortical structures (figure 9 b). Immediately behind the shock wave, vortex tubes are aligned in the spanwise direction, which has been observed in previous STI studies (Larsson et al. 2013; Boukharfane et al. 2018). More importantly, the vortex tubes also become aligned with the density iso-surfaces, meaning that the vorticity becomes perpendicular to the density gradient. This is consistent with the earlier analysis of the baroclinic torque. As a consequence, the post-shock vorticity field enhances the mixing between adjacent density regions. This coupling is further explored in the next section.

Figure 8. The PDF of the orientation between the vorticity vector and density gradient in $y$ – $z$ direction immediately after the shock wave.

For the strain rate tensor, figure 7 shows that its magnitude tends to be stronger in the heavy-fluid region and weaker in the light-fluid region. This trend is hypothesized to be related to the dependence of shock strength on the pre-shock density. Tian et al. (2017a ) showed that shock compression is stronger in the heavy-fluid region, while it is weaker in the smallest-density region, leading to the observed trend in the amplification of the magnitude of the strain rate tensor. This trend is different from that observed for the vorticity, which is explained above. As a result, the trend of the strain-enstrophy angle PDF peaking at around $45^{\circ }$ , observed in the single-fluid case at higher shock Mach numbers, is weakened in the multi-fluid case. Identifying the specific mechanisms behind variable-density turbulence interactions with shock waves, such as shock intensity dependence on density, density gradient effects, inertial effects and so on, can potentially be beneficial for modelling variable-density STI.

Figure 9. Vortex structures captured using the $Q$ -criterion, coloured by density, for multi-fluid (a) pre-shock turbulence and (b) post-shock turbulence.

3.1.2 Evolution of turbulence state downstream of the shock

The evolution of variable-density turbulence away from the shock wave involves many coupled nonlinear processes. In this subsection, the focus is on the evolution of turbulence structures.

Figure 10 shows the development of some of the fundamental turbulence statistics. The study of the evolution of these statistics helps in the understanding of the general characteristics of single- and multi-fluid STI. Figure 10(a) shows that with the introduction of strong density variations, the shock amplification of dissipation rate is stronger. Figure 10(b) shows the fluctuating pressure variance as a function of the streamwise position to highlight the development of the acoustic field. The amplification of the pressure fluctuations across the shock wave is noted, agreeing with Sethuraman et al. (2018). The acoustic wave is stronger in the multi-fluid case immediately after the shock wave. This is related to the shock intensity fluctuations induced by the strong density variations. As a result, the decay of the acoustic field is also faster for the multi-fluid case, causing a faster increase in turbulence kinetic energy (TKE). After the post-shock transient pressure adjustment, the multi-fluid case still exhibits larger absolute pressure fluctuations. However, after normalizing with $\unicode[STIX]{x1D70C}\overline{u^{\prime }u^{\prime }}$ , the pressure fluctuations become somewhat similar in magnitude in these two cases. In figure 10(c), the vortex stretching term $\unicode[STIX]{x1D6F4}=\overline{\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D714}_{j}(\unicode[STIX]{x2202}u_{i}/\unicode[STIX]{x2202}x_{j})}$ is decomposed into its streamwise $\unicode[STIX]{x1D6F4}_{x}=\overline{\unicode[STIX]{x1D714}_{1}\unicode[STIX]{x1D714}_{j}(\unicode[STIX]{x2202}u_{1}/\unicode[STIX]{x2202}x_{j})}$ and spanwise $\unicode[STIX]{x1D6F4}_{yz}=\overline{\unicode[STIX]{x1D714}_{2}\unicode[STIX]{x1D714}_{j}(\unicode[STIX]{x2202}u_{2}/\unicode[STIX]{x2202}x_{j})}$ components to explore the axisymmetric state and return-to-isotropy of post-shock turbulence. Previous studies (Livescu & Ryu 2016) have demonstrated that the normalized vortex stretching term reaches a low value after passing through the shock wave, indicating a tendency towards an axisymmetric state. Without normalization, figure 10(c) shows that the absolute values of the vortex stretching terms are magnified in both single- and multi-fluid cases, more so for the spanwise component. The two components then undergo a transient process, where they first increase and cross each other, before the flow returns to an isotropic state.

In order to quantitatively study the evolution of turbulence anisotropy, we consider here the Reynolds stress anisotropy tensor defined as $\unicode[STIX]{x1D623}_{ij}=\overline{u_{i}^{\prime }u_{j}^{\prime }}/\overline{u_{k}^{\prime }u_{k}^{\prime }}-\unicode[STIX]{x1D6FF}_{ij}/3$ . A similar anisotropy tensor, $\unicode[STIX]{x1D625}_{ij}$ , can also be defined for the vorticity field, as $\unicode[STIX]{x1D625}_{ij}=\overline{\unicode[STIX]{x1D714}_{i}^{\prime }\unicode[STIX]{x1D714}_{j}^{\prime }}/\overline{\unicode[STIX]{x1D714}_{k}^{\prime }\unicode[STIX]{x1D714}_{k}^{\prime }}-(\unicode[STIX]{x1D6FF}_{ij}/3)$ . Due to the homogeneity in spanwise directions, the diagonal components of the anisotropy tensor are related by $\unicode[STIX]{x1D623}_{22}=\unicode[STIX]{x1D623}_{33}=-0.5\unicode[STIX]{x1D623}_{11}$ , so only $\unicode[STIX]{x1D623}_{11}$ is discussed. The near-zero value of $\unicode[STIX]{x1D623}_{11}\approx 0$ is an indication that flow has reached an isotropic state, while $\unicode[STIX]{x1D623}_{11}\approx -1/3$ means that the turbulent field has a tendency towards a two-dimensional axisymmetric state. Figure 10(d) shows that $d_{11}$ , a small-scale turbulent variable, attains a value near $-0.3$ in the multi-fluid case, which is lower than that observed for the single-fluid case. This indicates that density intensifies the trend towards an axisymmetric state for small-scale turbulence. On the other hand, the stronger turbulent stretching mechanism as observed in figure 10(c) makes the return to isotropy much faster in the multi-fluid case as compared to that in the single-fluid case. For Reynolds stresses, large-scale turbulent variables, the multi-fluid flow reaches a quasi-isotropic state immediately after the shock wave ( $\unicode[STIX]{x1D623}_{11}\approx 0$ ), while single-fluid turbulence exhibits a tendency towards an axisymmetric state. This is in good agreement with Boukharfane et al. (2018). Evidently, the variable-density effects on the post-shock turbulence appear differently at small and large scales. Additionally, the quasi-isotropic state of the multi-fluid turbulence is not stable and is modified in the post-shock transition. Due to the energy transfer between the acoustic field and solenoidal turbulence field, $R_{11}$ quickly increases, causing $\unicode[STIX]{x1D623}_{11}$ to become larger than zero. The anisotropy reaches its maximum value at around the peak TKE position ( $k_{0}x\approx 2.0$ ) and then slowly decreases. For the single-fluid case, $\unicode[STIX]{x1D623}_{11}$ keeps increasing until $k_{0}x\approx 13.0$ , even though the acoustic effects almost vanish after peak TKE location of $k_{0}x\approx \unicode[STIX]{x03C0}$ .

Figure 10. Development of (a) turbulence dissipation rate, (b) pressure variance, (c) vortex stretching and (d) anisotropy ( $\unicode[STIX]{x1D623}_{11}$ ) of Reynolds stress and vorticity.

In figure 11, the developments of skewness and flatness of the longitudinal velocity gradients are examined before and after the flow interaction with the shock wave. They show how the non-Gaussian behaviours of the velocity field and specifically the VGT are affected by the combined shock and density effects. For isotropic turbulence, the skewness of the longitudinal velocity gradient should be around $-0.5$ , which is observed to be true in the pre-shock region for both single- and multi-fluid cases for both streamwise as well as spanwise components (figure 11 a). Immediately after the shock, different components of the derivative skewness tensor are shown to be modified in different ways. The streamwise component for both single-fluid and multi-fluid cases approaches values very close to 0, which is consistent with the tendency towards a two-dimensional axisymmetric state observed above. As the turbulence evolves away from the shock wave, the streamwise velocity derivative skewness decreases rapidly. Due to the strong density variations, the multi-fluid case exhibits a faster decrease in skewness before $k_{0}x=5.0$ , after which it slowly increases towards a value of $-0.54$ . The shock modification of the skewness of the transverse derivative is relatively small for the single-fluid case. For the multi-fluid case, the longitudinal transverse velocity derivative becomes less negatively skewed, with a value of around $-0.25$ . This difference can be attributed to stronger shock intensity variations and shock wrinkling in the multi-fluid case. Away from the shock wave, for both cases, the skewness of $\unicode[STIX]{x2202}v/\unicode[STIX]{x2202}y$ increases first until it reaches a peak and then slowly decreases. Comparably, the multi-fluid case exhibits a shorter but more intensified transition. At the end of the domain, however, the spanwise derivative skewness is still larger than $-0.5$ , as the flow is still anisotropic. Figure 11(b) shows the development of longitudinal velocity derivative flatness factor across and after the shock wave. Immediately after the shock, the flatness of the streamwise component decreases in value while that of the spanwise component increases. Similar to the skewness, the effect of density variations is relatively small on the flatness of the streamwise component for the Atwood number considered in this study. On the other hand, the density variations in the multi-fluid case make the increase in flatness of the transverse component less significant, with the pre- and post-shock values being almost the same. Away from the shock wave, the flatness of the longitudinal streamwise velocity derivative increases, returning to its pre-shock value, while the growth is much faster in the multi-fluid case. For the transverse longitudinal derivative component, the flatness slowly decreases after a small change.

Figure 11. Development of (a) skewness and (b) flatness of the streamwise and transverse components of velocity derivatives.

Figure 12. The PDF of the density gradient at different streamwise locations for (a) single-fluid case and (b) multi-fluid case.

From the results above, it can be stated that the variable-density effects are not strongly manifested immediately after the shock wave for some of the statistics, but they play an important role in the post-shock adjustment. It is possible for these statistics that the dominating effect across the shock is the shock compression. However, the density variations cause differences in the post-shock turbulence structure, which affect the turbulence development away from the shock wave. To get an insight into this behaviour, density gradient PDFs are examined at various streamwise positions in figure 12. Before the shock wave, the PDFs of the density gradients are symmetric in all three directions for both single- and multi-fluid cases (not shown). For the single-fluid case, after passing through the shock wave, the density gradient PDFs remain symmetric, but the streamwise component PDF becomes wider due to the shock compression (Boukharfane et al. 2018). As the turbulence develops away from the shock wave towards the peak TKE position, the density gradient PDFs still remain symmetric and become narrower, which is related to the fast decay of the acoustic field. For the multi-fluid case, the density gradient PDFs are strongly amplified through the shock wave, but the changes are relatively small far from the shock, because the density variations are controlled by the mixture composition instead of the acoustic field. More importantly, the streamwise component becomes negatively skewed.

To identify the mechanisms responsible for the skewness of the streamwise density gradient, we examine the orientation of the eigenvectors of strain rate tensor $\unicode[STIX]{x1D61A}_{ij}$ . The PDFs of the cosines of the angles between the three eigenvectors and the streamwise direction, conditioned on regions with positive or negative density gradients, are plotted in figure 13. The eigenvalues of the strain rate tensor are $\unicode[STIX]{x1D6FE}_{1}$ , $\unicode[STIX]{x1D6FE}_{2}$ and $\unicode[STIX]{x1D6FE}_{3}$ , where $\unicode[STIX]{x1D6FE}_{1}<\unicode[STIX]{x1D6FE}_{2}<\unicode[STIX]{x1D6FE}_{3}$ . The angles between these eigenvectors and the streamwise direction are denoted by $\unicode[STIX]{x1D712}_{1}$ , $\unicode[STIX]{x1D712}_{2}$ and $\unicode[STIX]{x1D712}_{3}$ . For the multi-fluid case, in the positive density gradient regions, the extensive ( $\unicode[STIX]{x1D6FE}_{3}$ ) eigenvector is more likely to be aligned with the streamwise direction (figure 13 a). The intermediate ( $\unicode[STIX]{x1D6FE}_{2}$ ) eigenvector is misaligned with the streamwise direction and the compressive ( $\unicode[STIX]{x1D6FE}_{1}$ ) eigenvector has no preferential alignment. This implies that the density field is generally being stretched in the streamwise direction, making the magnitude of the density gradient smaller. On the other hand, the alignment of the $\unicode[STIX]{x1D6FE}_{1}$ and $\unicode[STIX]{x1D6FE}_{3}$ eigenvectors with the streamwise directions is reversed in the negative density gradient regions as shown in figure 13(b). The density field is then compressed so that the magnitude of the density gradient is increased. This asymmetry in the alignment is caused by the nonlinear variable-density effects when the flow passes through the shock wave and explains the negatively skewed PDF of density gradient in the multi-fluid case. It is also interesting to note the different roles of density gradient across the shock wave: spanwise density gradients contribute to the generation of the vorticity field, while the streamwise component affects the strain field. For the single-fluid case, the asymmetry in the eigenvector behaviour is small and vanishes quickly away from the shock wave. This implies that even though density variations may not affect some of the turbulence statistics directly, they modify the topology and structure of turbulence immediately after the shock and continue to manifest their effects in the post-shock turbulence evolution.

Figure 13. The PDFs of the cosine angle between eigenvectors of the strain rate tensor and streamwise $(x)$  axis for regions with (a) $\text{d}\unicode[STIX]{x1D70C}/\text{d}x>0$ and (b) $\text{d}\unicode[STIX]{x1D70C}/\text{d}x<0$ and for multi-fluid (solid lines) and single-fluid (dashed lines) cases.

3.2 Topological analysis of the post-shock turbulence

To further characterize the turbulence structure behind the shock wave, we have analysed the invariant space of the VGT. The second and third invariants (denoted by $Q^{\ast }$ and $R^{\ast }$ ) of the anisotropic/deviatoric part of the VGT can reveal important features of the flow topology (Pirozzoli & Grasso 2004). In highly compressible turbulence, there exits a richer set of flow topologies due to the dilatational part of the VGT (Suman & Girimaji 2010). For the parameter range considered in this study, however, the compressibility effects are weak. This is demonstrated in figure 14, where the normalized PDFs of the dilatation and vorticity for pre-shock isotropic turbulence, single-fluid post-shock turbulence and multi-fluid post-shock turbulence are shown. The pre-shock isotropic turbulence has a very low magnitude of dilatation. The shock wave expectedly amplifies the dilatation magnitude, and more so when variable-density effects exist, but the dilatation values are still considerably lower than those studied in Suman & Girimaji (2010), Chu & Lu (2013) and Vaghefi & Madnia (2015). Considering that the focus of this study is on the variable-density effects, here we only present the topological structure of the anisotropic VGT. The anisotropic part of the VGT is calculated using the formula $\unicode[STIX]{x1D608}_{ij}^{\ast }=\unicode[STIX]{x1D608}_{ij}-\unicode[STIX]{x1D703}/3I$ . Correspondingly, the second and third invariants can be calculated from

(3.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle Q^{\ast }=-{\textstyle \frac{1}{2}}A_{ij}^{\ast }A_{ji}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle R^{\ast }=-{\textstyle \frac{1}{3}}A_{ij}^{\ast }A_{jk}^{\ast }A_{ki}^{\ast }. & \displaystyle\end{eqnarray}$$

Similarly, the invariants of the symmetric and skew-symmetric parts of the anisotropic VGT, $\unicode[STIX]{x1D61A}_{ij}^{\ast }$ and $\unicode[STIX]{x1D61E}_{ij}^{\ast }$ , can also be calculated using the corresponding form of (3.1). They are denoted here as ( $Q_{s}^{\ast }$ , $R_{s}^{\ast }$ ) and ( $Q_{w}^{\ast }$ , $R_{w}^{\ast }$ ). The following equations relate the above variables for the anisotropic part of the VGT (Ooi et al. 1999):

(3.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle Q^{\ast }=Q_{s}^{\ast }+Q_{w}^{\ast }, & \displaystyle\end{eqnarray}$$

(3.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle R^{\ast }=R_{s}^{\ast }-\unicode[STIX]{x1D714}_{i}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D714}_{j}^{\ast }, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D714}_{i}^{\ast }=\unicode[STIX]{x1D714}_{i}$ is the vorticity vector. The scalar variables $Q_{s}^{\ast }$ and $Q_{w}^{\ast }$ are related to the local dissipation rate ( $-Q_{s}^{\ast }=1/2\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }$ ) and enstrophy ( $Q_{w}^{\ast }=1/2\unicode[STIX]{x1D61E}_{ij}^{\ast }\unicode[STIX]{x1D61E}_{ij}^{\ast }$ ), respectively. For constant viscosity, $Q^{\ast }$ represents the difference between enstrophy and dissipation (Chu & Lu 2013). Similarly, $R_{s}^{\ast }=-1/3\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D61A}_{jk}^{\ast }\unicode[STIX]{x1D61A}_{ki}^{\ast }$ is related to the production of dissipation due to strain field and $\unicode[STIX]{x1D714}_{i}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }\unicode[STIX]{x1D714}_{j}^{\ast }$ is the vortex stretching contribution to the enstrophy. Therefore, for constant viscosity, $R^{\ast }$ represents the difference between enstrophy production and dissipation production. Based on the local values of $Q^{\ast }$ and $R^{\ast }$ , four types of local flow topologies can be identified: stable-focus/stretching (SFS), unstable-focus/contracting (UFC), stable-node/saddle/saddle (SN/S/S) and unstable-node/saddle/saddle (UN/S/S). For isotropic turbulence, the joint PDF of ( $Q^{\ast },R^{\ast }$ ) has a tear-drop shape. This has been further observed in other fully developed turbulent flows, such as boundary layers, mixing layers and channel flows (Pirozzoli & Grasso 2004; Wang et al. 2012). This type of distribution of $Q^{\ast }$ and $R^{\ast }$ is an indicator that the turbulence is more likely to have a local topology of SFS or UN/S/S. In figure 15(a), it is shown that the joint PDF of normalized second and third invariants, $Q^{\ast }/\langle Q_{w}\rangle$ and $R^{\ast }/\langle Q_{w}\rangle ^{3/2}$ , has the same tear-drop shape in the pre-shock flow. Using Shock-LIA and DNS data, Ryu & Livescu (2014) showed that for single-fluid STI, the ( $Q^{\ast },R^{\ast }$ ) distribution is significantly modified by the shock wave, with a tendency towards symmetrization of the joint PDF. This indicates that the regions with UFC and SN/S/S (first and third quadrants) are more likely to occur as turbulence develops a two-dimensional axisymmetric flow structure. To understand the variable-density effects on this shock-induced symmetrization, the joint PDFs of ( $Q^{\ast }$ , $R^{\ast }$ ) for both single-fluid and multi-fluid post-shock turbulence are compared in figure 15(b,c).

Figure 14. The PDFs of the normalized dilatation and vorticity for isotropic turbulence (IT), single-fluid post-shock turbulence (s) and multi-fluid post-shock turbulence (m).

Figure 15(b) shows the joint distribution for the single-fluid post-shock turbulence. The dashed lines denote the locus of zero discriminant of $A^{\ast }$ , where $Q^{\ast }$ and $R^{\ast }$ satisfy $27R^{\ast 2}/4+Q^{\ast 3}=0$ . Compared to the pre-shock joint PDF, there is a tendency towards symmetrization, with more points located in the first and third quadrants. Similar to single-fluid STI, multi-fluid STI demonstrates a tendency towards symmetrization of the $(Q^{\ast },R^{\ast })$ distribution. However, the multi-fluid distribution is slightly more symmetric and has a larger variance, with more points away from the axes. This implies that more extreme ‘events’ exist in the post-shock multi-fluid turbulence.

Figure 15. Iso-contour lines of joint PDFs of normalized second and third invariants of the anisotropic part of the VGT, ( $Q^{\ast }$ , $R^{\ast }$ ), for (a) pre-shock, (b) single-fluid post-shock turbulence and (c) multi-fluid post-shock turbulence. The lateral lines denote the locus of zero discriminant.

The density effects on the post-shock joint PDF of second and third invariants are further explored by comparing the conditional distribution, conditioned on regions with different densities, in figure 16(a–c). Figure 16(a) corresponds to regions with relatively high density ( $\unicode[STIX]{x1D70C}>(\overline{\unicode[STIX]{x1D70C}}+90\,\%\unicode[STIX]{x1D70C}_{rms}^{\prime })$ ), figure 16(b) to regions with density around the post-shock mean value and figure 16(c) to low-density regions ( $\unicode[STIX]{x1D70C}<(\overline{\unicode[STIX]{x1D70C}}-90\,\%\unicode[STIX]{x1D70C}_{rms}^{\prime })$ ). For consistency check, the joint PDFs corresponding to these regions are also computed for the pre-shock flow (not shown) and found to be close to the single-fluid PDFs. After the shock wave, the joint PDFs demonstrate significant differences between regions with different densities. In regions with density closer to that of the post-shock mean density, the distribution of invariants appears to be very similar to that shown in figure 15(c). But for regions with higher density (figure 16 a), the joint PDF becomes more symmetric compared to the overall flow for both the single-fluid and multi-fluid cases. There is a much larger portion of data points having a local topology of SN/S/S, and fewer data points fall into the first and second quadrants, indicating larger strain-dominated regions. On the other hand, the post-shock regions with low density values (figure 16 c) exhibit features similar to that of isotropic turbulence, with almost the same tear-drop shape, only with a larger variance or a wider distribution. The quantitative difference is hypothesized to be related to the higher shock strength variation in the multi-fluid case. It was observed in our previous studies (Tian, Jaberi & Livescu 2019) that the local shock strength is positively correlated with the pre-shock density. With a stronger shock, the two-dimensionalization effect on the post-shock turbulence should also appear stronger in the high-density regions (Livescu & Ryu 2016). For low-density regions, the smaller two-dimensionalization effect reduces the symmetrization trend. Moreover, the relatively lower inertia in these regions leads to a faster response to the local strain field (Livescu et al. 2010), which could make a faster return to isotropic turbulence. The different characteristics of the ( $Q^{\ast },R^{\ast }$ ) joint PDF in regions with different densities provide additional evidence for the previous argument made about the role of density in the preferential amplification of the strain and rotation tensors.

Figure 16. Iso-contour lines of post-shock ( $k_{0}x\approx 0.44$ ) joint PDF of second and third invariants of the anisotropic part of the VGT, ( $Q^{\ast }$ , $R^{\ast }$ ), in regions with different densities. (a) Regions with high density values, $\unicode[STIX]{x1D70C}>(\overline{\unicode[STIX]{x1D70C}}+90\,\%\unicode[STIX]{x1D70C}_{rms}^{\prime })$ , (b) regions with density around the post-shock mean value and (c) regions with low density values, $\unicode[STIX]{x1D70C}<(\overline{\unicode[STIX]{x1D70C}}-90\,\%\unicode[STIX]{x1D70C}_{rms}^{\prime })$ .

Figure 17. Colour illustration of the flow topology for multi-fluid STI. The flow topology is represented by the quadrants (denoted as $Q_{i}$ ) of the joint PDF of ( $Q^{\ast }$ , $R^{\ast }$ ). (a) Two-dimensional colour contours in the $x$ – $z$ plane at $y=3.14$ (half $y$ -domain). The shock wave is located in the middle of the domain at $k_{0}x\approx 0$ . The ratio of the fluid volume in different quadrants in the pre-shock region is $Q_{1}:Q_{2}:Q_{3}:Q_{4}=26.7\,\%:38.7\,\%:7.8\,\%:26.8\,\%$ . The two-dimensional colour contours in the homogeneous $y$ – $z$ plane at streamwise locations of (b)  $k_{0}x\approx 0.2$ $(28.7\,\%:34.4\,\%:14.3\,\%:22.6\,\%)$ , (c)  $k_{0}x\approx 2.0$ , peak TKE location $(26.7\,\%:36.9\,\%:11.2\,\%:25.2\,\%)$ and (d)  $k_{0}x\approx 4.0$ $(26.3\,\%:37.9\,\%:9.3\,\%:26.2\,\%)$ .

In figure 17, the planar distributions of the flow topologies are shown. Here, $Q_{i}$ refers to the quadrants on the joint PDF of ( $Q^{\ast }$ , $R^{\ast }$ ), which amounts to a representation of the local flow topology. Figure 17(a) presents the two-dimensional visualization of the flow topology in a typical $x$ – $z$ plane. The regions occupied by different quadrants are marked using different colours. Evidently, the vorticity-dominated regions ( $Q_{1}$ and $Q_{2}$ ) cover a large portion of the flow and have more compact shapes. These regions are connected by UN/S/S areas ( $Q_{4}$ ), which are more elongated. The SN/S/S ( $Q_{3}$ ) areas can be located either at the edge of $Q_{4}$ regions or between different $Q_{4}$ regions. These strain-dominated regions seem to be more fragmented than the compact vorticity-dominated regions. Figure 17(b–d) shows the two-dimensional contours of $Q_{i}$ in the homogeneous ( $y$ – $z$ ) planes at different streamwise locations after the shock. These locations are also marked on figure 17(a). Immediately after passing through the shock wave, the volume fractions of different quadrants are calculated to be $Q_{1}:Q_{2}:Q_{3}:Q_{4}=28.7\,\%:34.4\,\%:14.3\,\%:22.6\,\%$ , indicating a trend towards symmetrization in the joint PDF. This can also be observed in the two-dimensional contours in figure 17(b). Moreover, the characteristic length scales associated with the regions occupied by different quadrants are decreased across the shock wave. As the flow evolves away from the shock wave, the distribution slowly changes back to the pre-shock shape but still with smaller turbulence length scales. Most of fluid in different quadrants returns to the pre-shock values at $k_{0}x=4.0$ . The reorientation of the flow structures into the streamwise direction is also noted in figure 17(a), consistent with the return to isotropy of the flow. However, the rates at which different flow features return to an isotropic state are slightly different. The dynamics of flow and the return-to-isotropic turbulence process are examined in detail in the next section using Lagrangian statistics.

The quasi-axisymmetric state immediately after the shock wave, identified above based on the joint PDF of ( $Q^{\ast }$ , $R^{\ast }$ ), is further explored below by considering the vortex stretching rate and other flow topological features.

The rate at which the vorticity is stretched or contracted, i.e. the normalized vortex stretching rate, can be calculated based on the VGT invariants using the formula $\unicode[STIX]{x1D6F4}^{\ast }=w_{i}^{\ast }\unicode[STIX]{x1D61A}_{ij}^{\ast }w_{j}^{\ast }=(R_{s}^{\ast }-R^{\ast })/Q_{w}^{\ast }$ (Ooi et al. 1999). In figure 18, the joint PDF of ( $-Q_{s}^{\ast }$ , $\unicode[STIX]{x1D6F4}^{\ast }$ ) is plotted to investigate the effects of the strain field on the vortex stretching rate. Positive and negative $\unicode[STIX]{x1D6F4}^{\ast }$ values correspond to the vortex being stretched or contracted. Figure 18(a) shows the joint PDF of ( $-Q_{s}^{\ast }$ , $\unicode[STIX]{x1D6F4}^{\ast }$ ) for the isotropic turbulence. The results agree very well with those of Ooi et al. (1999), which indicates that the flow favours positive $\unicode[STIX]{x1D6F4}^{\ast }$ values or an overall vortex stretching, especially in the strong strain-dominated regions. Here, we compare the results from isotropic turbulence to those from single-fluid and multi-fluid post-shock turbulence to understand the shock and variable-density effects. We note that in figure 18(b), the joint PDF becomes more symmetric around $\unicode[STIX]{x1D6F4}^{\ast }=0$ after passing through the shock wave. For the multi-fluid case, as shown in figure 18(c), the joint PDF becomes almost fully symmetric, especially at lower $-Q_{s}^{\ast }$ values. This symmetry has a strong effect on the overall vortex stretching rate for the multi-fluid post-shock turbulence because the positive and negative $\unicode[STIX]{x1D6F4}^{\ast }$ values tend to cancel each other through averaging. Moreover, the variances of the stretching term are almost the same for single- and multi-fluid cases, meaning that the lower stretching rate is mainly due to changes in the turbulence structure (especially in more negative $\unicode[STIX]{x1D6F4}^{\ast }$ regions), and not simply to the decrease in the magnitude of $\unicode[STIX]{x1D6F4}^{\ast }$ .

To understand the contribution from different topological states to the vortex stretching, the joint PDF of ( $-Q_{s}^{\ast }$ , $\unicode[STIX]{x1D6F4}^{\ast }$ ) is conditioned on different quadrants for the multi-fluid case. Figure 19(a,b) shows the joint PDF for $Q_{1}$ and $Q_{2}$ regions, or areas with a local topology of SFS and UFC. It can be observed that in these rotation-dominated regions, the magnitude of the vortex stretching rate $\unicode[STIX]{x1D6F4}^{\ast }$ is relatively small. Moreover, $Q_{1}$ is dominated by areas with negative $\unicode[STIX]{x1D6F4}^{\ast }$ and $Q_{2}$ is dominated by positive $\unicode[STIX]{x1D6F4}^{\ast }$ areas, which can be inferred from their corresponding topologies. However, for $Q_{3}$ and $Q_{4}$ (figure 19 c,d), the magnitude of the vortex stretching rate is larger than that in the rotation-dominated regions (figure 19 a,b). In $Q_{3}$ , the joint PDF is relatively symmetric and seems to be slightly biased towards negative vortex stretching rate values. Quadrant $Q_{4}$ , on the other hand, is dominated by positive vortex stretching. Overall, the results explain the lower averaged vortex stretching rate values in the multi-fluid case caused by the enhancement of $Q_{1}$ and $Q_{3}$ regions.

Figure 18. Iso-contour lines of joint PDF of ( $-Q_{s}^{\ast }$ , $\unicode[STIX]{x1D6F4}^{\ast }$ ) for (a) isotropic box turbulence and (b) single-fluid and (c) multi-fluid turbulence at post-shock position of $k_{0}x\approx 0.44$ .

Figure 19. Iso-contour lines of joint PDF of ( $-Q_{s}^{\ast }$ , $\unicode[STIX]{x1D6F4}^{\ast }$ ) for different quadrants right after the shock wave: (a)  $Q_{2}$ , (b)  $Q_{1}$ , (c)  $Q_{3}$ and (d)  $Q_{4}$ .

3.3 Lagrangian dynamics of VGT

In this section, we use non-inertial Lagrangian particles/tracers that move with the local flow velocity because their statistics can reflect important transient turbulent dynamics, which is difficult to study using Eulerian data (Yeung 2002). More importantly, in the context of variable-density flows, the Lagrangian statistics enable us to differentiate among particles with different densities and investigate their dynamics separately. Lagrangian data are used here to perform an analysis of the post-shock turbulence structure and VGT dynamics with and without significant density fluctuations.

The first result considered is the time scale of particle motions related to different flow topologies. In figure 20, the percentages of fluid particles that remain in their starting quadrants are plotted over time so that we can identify the residence time of particles for different turbulence structures. In figure 20(a), the percentages of fluid particles are plotted for decaying isotropic turbulence as a reference. It is noted that $Q_{3}$ and $Q_{4}$ , which are the strain-dominated regions, have the smallest associated residence times among the four quadrants, with $Q_{3}$ time being the smaller of the two. For the rotation-dominated regions, the residence times are expectedly longer, especially for $Q_{2}$ . The residence times for single-fluid and multi-fluid simulations can be inferred from figure 20(b,c). For both cases, $Q_{3}$ always has the least residence time and $Q_{2}$ has the largest one. Comparing all three cases, the particles in the multi-fluid and single-fluid post-shock turbulence are shown to evolve faster away from the original quadrant than particles in isotropic turbulence, indicating smaller time scales of the flow topologies. Between the two post-shock turbulence fields, the multi-fluid case presents shorter residence times.

Figure 21 presents an example of the temporal development of the above-mentioned structures. The evolution of a vortex tube in the post-shock turbulence is tracked and visualized as it moves away from the shock wave in figure 21(a). As expected, the depicted vortical structure maintains its shape, except that it is being stretched and reoriented by the local flow field. Moreover, the vortex tube surface is almost parallel with the iso-surface of the density field, i.e. perpendicular to the density gradient. This is consistent with the discussion in § 3.1.1 regarding the bulk of vorticity generation across the shock wave. As the vortex tube evolves away from the shock wave, the reorientation of the density gradient by the vortex is also observed. In figure 21(b), a strain-dominated structure is visualized using the iso-surface of negative $Q^{\ast }$ . It can be clearly seen that such structures lack temporal coherency since they tend to be become fragmented as they evolve.

Figure 20. Percentage of fluid particles that stay in each quadrant following particles initialized uniformly in (a) isotropic turbulence and (b) single-fluid and (c) multi-fluid turbulence at post-shock position of $k_{0}x=0.44$ .

Figure 21. Visualization of the temporal development (left to right) of the turbulence structure using iso-surfaces of $Q^{\ast }$ coloured by density for multi-fluid post-shock turbulence. These structures are captured immediately after the shock wave. (a) Vorticity-dominated structure and (b) strain-dominated structure.

Figure 22. Contributions to the vortex stretching rate from particles starting in each quadrant. The particles are initialized uniformly at the post-shock position $k_{0}x\approx 0.44$ and traced downstream until the vorticity returns to an isotropic state. (a) Single-fluid case and (b) multi-fluid case.

In figure 22, the contributions to the normalized vortex stretching rate from particles that are initialized in each of the four quadrants are plotted following these particles. As expected, at $t=0$ , particles from $Q_{2}$ and $Q_{4}$ add positively to the vortex stretching rate, while those from $Q_{1}$ have a negative vortex stretching rate contribution on average. This is in good agreement with the joint PDFs of ( $-Q_{s}^{\ast }$ , $\unicode[STIX]{x1D6F4}^{\ast }$ ), shown in figure 19. For $Q_{3}$ , the initial contribution is close to zero. As the fluid particles move with the turbulent flow, their contributions to the vortex stretching also change. It can be seen that the fast increase in vortex stretching can be mainly attributed to the particles originating in $Q_{1}$ and $Q_{3}$ . Particles starting in $Q_{4}$ have an increasing vortex stretching contribution for a short period before their combined/average contributed value starts to decrease. The behaviour is qualitatively similar for the single-fluid case, but the changes are smaller in this case. For both cases, the vortex stretching contribution from the initial $Q_{2}$ particles decreases in time.

To further understand this behaviour, the Lagrangian equations of the VGT and its invariants are considered. The time evolution of $\unicode[STIX]{x1D608}_{ij}$ for fluid particles can be obtained by taking the spatial derivatives of the Navier–Stokes equations. In dimensionless form, it can be written as (Chu & Lu 2013)

(3.3a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D608}_{ij}}{\unicode[STIX]{x2202}t}+u_{k}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D608}_{ij}}{\unicode[STIX]{x2202}x_{k}}+A_{ik}A_{kj}=-H_{ij}+{\mathcal{T}}_{ij}, & \displaystyle\end{eqnarray}$$

(3.3b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\unicode[STIX]{x1D608}_{ij}}{\text{D}t}=-A_{ik}A_{kj}-H_{ij}+{\mathcal{T}}_{ij}, & \displaystyle\end{eqnarray}$$

with

(3.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle H_{ij}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}\right)=-\frac{1}{\unicode[STIX]{x1D70C}^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}x_{j}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}p^{2}}{\unicode[STIX]{x2202}x_{i}\unicode[STIX]{x2202}x_{j}}=H_{ij}^{b}+H_{ij}^{p}, & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{T}}_{ij}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ik}}{\unicode[STIX]{x2202}x_{k}}\right), & \displaystyle\end{eqnarray}$$

(3.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}_{ij}=\frac{\unicode[STIX]{x1D707}}{Re_{0}}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}-\frac{2}{3}\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{k}}\unicode[STIX]{x1D6FF}_{ij}\right), & \displaystyle\end{eqnarray}$$

where $Re_{0}$ is the reference Reynolds number. From here, the dynamic equations for the three invariants of the VGT, $P$ , $Q$ and $R$ , can be derived in the following form (Chu & Lu 2013):

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}P}{\text{D}t}=(P^{2}-2Q)+H_{ii}^{p}+H_{ii}^{b}-{\mathcal{T}}_{ii}, & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}Q}{\text{D}t}=(PQ-3R)+(PH_{ii}^{p}+\unicode[STIX]{x1D608}_{ij}H_{ji}^{p})+(PH_{ii}^{b}+\unicode[STIX]{x1D608}_{ij}H_{ji}^{b})+(-P{\mathcal{T}}_{ii}-\unicode[STIX]{x1D608}_{ij}{\mathcal{T}}_{ji}), & \displaystyle\end{eqnarray}$$

(3.5c ) $$\begin{eqnarray}\displaystyle \frac{\text{D}R}{\text{D}t} & = & \displaystyle PR+(QH_{ii}^{p}+P\unicode[STIX]{x1D608}_{ij}H_{ji}^{p}+\unicode[STIX]{x1D608}_{ij}\unicode[STIX]{x1D608}_{jk}H_{ki}^{p})\nonumber\\ \displaystyle & & \displaystyle +\,(QH_{ii}^{b}+P\unicode[STIX]{x1D608}_{ij}H_{ji}^{b}+\unicode[STIX]{x1D608}_{ij}\unicode[STIX]{x1D608}_{jk}H_{ki}^{b})+(-Q{\mathcal{T}}_{ii}-P\unicode[STIX]{x1D608}_{ij}{\mathcal{T}}_{ji}-\unicode[STIX]{x1D608}_{ij}\unicode[STIX]{x1D608}_{jk}{\mathcal{T}}_{ki}),\end{eqnarray}$$

where the three invariants of the VGT are defined as

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle P=-\text{tr}(\unicode[STIX]{x1D608}_{ij}), & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle Q={\textstyle \frac{1}{2}}(\text{tr}(\unicode[STIX]{x1D608}_{ij})^{2}-\text{tr}(\unicode[STIX]{x1D608}_{ij}\unicode[STIX]{x1D608}_{jk})), & \displaystyle\end{eqnarray}$$

(3.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle R=-\text{det}(\unicode[STIX]{x1D608}_{ij}). & \displaystyle\end{eqnarray}$$

Here, $\text{tr}(\unicode[STIX]{x1D608}_{ij})$ and $\text{det}(\unicode[STIX]{x1D608}_{ij})$ denote the trace and determinant of a tensor. Note that instead of the deviatoric part of the VGT, the dynamic equations for the full VGT are considered. The reason is that due to the variable-density effects and shock compression, the incompressibility condition is not satisfied, especially when $M_{t}$ and $A_{t}$ become large. Even though $M_{t}$ and $A_{t}$ in this study are small, we still consider the full equations for any future comparisons. The contributions from the dilatational part of the VGT and their coupling with the variable-density effects in highly compressible turbulence are still unknown and need to be explored for STI in future studies.

The dynamical equations can be divided into contributions of four different parts: (1) mutual interaction among invariants, (2) pressure Hessian $H_{ij}^{p}$ , (3) baroclinic $H_{ij}^{b}$ and (4) viscous term ${\mathcal{T}}_{ij}$ . The statistics regarding these terms can be extracted from the Lagrangian data.

Figure 23. The PDFs of (a) $(\text{D}Q/\text{D}t)/\langle Q_{w}\rangle ^{3/2}$ and (b) $(\text{D}R/\text{D}t)/\langle Q_{w}\rangle ^{2}$ for fluid particles with different densities at streamwise location of $k_{0}x\approx 0.5$ .

Some general features of the Lagrangian dynamics of the VGT invariants are examined through the PDFs of their material derivatives. The variable-density effects can be identified by comparing the PDFs corresponding to regions with different densities (figure 23). In the light-fluid regions, the PDFs of $\text{D}Q/\text{D}t$ and $\text{D}R/\text{D}t$ have narrower tails, while the tails are wider in the heavy-fluid regions. Another important observation is that the skewness of $\text{D}Q/\text{D}t$ is different in the light- and heavy-fluid regions. Heavy-fluid particles have a positively skewed PDF, similar to the overall flow. On the other hand, the $\text{D}Q/\text{D}t$ skewness resulting from light-fluid particles is negative. This implies that heavy-fluid particles are more likely to move towards rotation-dominated regions and vice versa. These differences can be attributed to differences in the return-to-isotropy, experienced by fluid particles with different densities.

Figure 24. Conditional mean rate of change vectors of ( $\text{D}Q/\text{D}t/\langle Q_{w}\rangle ^{3/2}$ , $\text{D}R/\text{D}t/\langle Q_{w}\rangle ^{2}$ ) in the ( $Q,R$ ) plane for (a) isotropic turbulence, (b) single-fluid post-shock turbulence and (c) multi-fluid post-shock turbulence at streamwise location of $k_{0}x\approx 0.5$ . To ensure that the vectors can be properly visualized, their sizes are re-scaled by multiplying by a constant of 0.3. This applies to all the following vector plots.

The Lagrangian dynamics of the turbulence and the evolution of flow topology are further examined here by considering the conditional mean rate of change of $Q$ and $R$ in the invariants plane (Ooi et al. 1999). The rates of change are used to form a vector at each point in the invariants plane. The trajectories implied by these vectors can be followed to understand the return-to-isotropy process. In fully compressible turbulence, the $(P,Q,R)$ invariant space becomes three-dimensional (Suman & Girimaji 2010; Chu & Lu 2013; Vaghefi & Madnia 2015) and there exists an out-of-plane $(Q,R)$ component of the trajectory due to the contribution from the compressibility ( $P$ ) effect. Due to the low compressibility effect in this work, however, it would be more appropriate to consider only the in-plane $(Q,R)$ dynamics and leave the compressibility effects for future study. Therefore, the results presented below correspond to the data points with small magnitude of $P$ ( $P/\langle Q_{w}\rangle ^{0.5}<$ 0.1) for the relatively ‘incompressible’ region of the flow. These points comprise approximately $60\,\%$ of the flow.

The procedure used to obtain the conditional mean vectors (CMVs) in this study is similar to that in Ooi et al. (1999). Based on the conditional averages introduced in (2.2), $X(Q,R)$ represents a statistical quantity that is conditioned on $Q$ and $R$ . The statistical convergence concerning the bin sizes and the number of samples in each bin has been discussed in § 2.4.

The normalized CMVs ( $\text{D}Q/\text{D}t/\langle Q_{w}\rangle ^{3/2}$ , $\text{D}R/\text{D}t/\langle Q_{w}\rangle ^{2}$ ) for different flows are shown in figure 24. The vectors obtained from isotropic turbulence data are shown in figure 24(a) for reference. It can be seen that the CMVs exhibit a circulating behaviour in the $(Q,R)$ plot around the origin in the clockwise direction, indicating that the flow evolves from SFS to UFC, UN/S/S, SN/S/S then back to SFS on average. This circulating behaviour represents the Lagrangian dynamics in fully developed turbulence that maintains the tear-drop shape of the ( $Q,R$ ) distribution. This has been observed in many incompressible/compressible canonical turbulent flows (Ooi et al. 1999; Chu & Lu 2013). The CMVs for single-fluid and multi-fluid post-shock turbulence are shown in figure 24(b,c). Evidently, the joint PDF of ( $Q,R$ ) becomes more symmetric due to shock compression. From the Lagrangian point of view, the circulating behaviour as seen in figure 24(a) for isotropic turbulence is weakened. The particles in $Q_{2}$ tend to have an increasing $Q$ and decreasing $R$ , resulting in an overall trend of getting away from the original point, instead of circulating and then moving towards $Q_{1}$ . This trend in the second quadrant represents an increase of enstrophy. The particles in $Q_{1}$ have similar dynamics as in isotropic turbulence and tend to move downward in the ( $Q,R$ ) plane towards the zero discriminant curve. The particles in $Q_{3}$ are more likely to move straight up towards $Q_{2}$ , while those in $Q_{4}$ are likely to move away from the original point following the direction of the zero discriminant line and then circulate back to $Q_{3}$ . The overall behaviour exhibited by these particles demonstrates the return-to-isotropy process, with an enlarging head in the second quadrant and elongating tail in the fourth quadrant, anticipating the formation of the classic tear-drop shape.

Figure 25. The CMVs in the ( $Q,R$ ) invariants plane for (a) light fluid, (b) medium-density fluid and (c) heavy fluid at streamwise location of $k_{0}x\approx 0.5$ .

The density effects can be further examined by conditioning the ( $\text{D}Q/\text{D}t,\text{D}R/\text{D}t$ ) vector field on the local density. Figure 25(a) shows the CMVs for the light-fluid regions. The light-fluid particles retain the circulating motion, except that the particles in $Q_{3}$ and $Q_{4}$ are likely to go straight left instead of following the zero discriminant line. In general, the flow dynamics in the light-fluid regions is less affected by the shock wave. For the medium-density-fluid regions (figure 25 b), the circulating motion disappears. On the right-hand side of the ( $Q,R$ ) plane ( $R>0$ ), which is the strong dissipation-production region based on (3.2), the fluid particles tend to move downward, resulting in lower $Q$ values. On the left-hand side of the ( $Q,R$ ) plane ( $R<0$ ), which is the enstrophy-production-dominated region, the fluid particles tend to move to the left, indicating an increased enstrophy production. The overall downward-moving behaviour of the medium-density-fluid particles is indicative of decreasing vorticity. This is possibly due to the fact that vorticity is preferentially amplified in the medium-density region across the shock wave. After passing the shock wave, the vorticity will decrease as the correlation between density and vorticity vanishes. Figure 25(c) shows the CMVs for the heavy-fluid regions. Interestingly, the heavy-fluid particles exhibit counterclockwise motion. The heavy-fluid particles start from $Q_{3}$ and move to $Q_{4}$ , $Q_{1}$ , and finally to $Q_{2}$ . This implies that they become vorticity-dominated due to the fast depletion of strain. Evidently, density plays an important role in the development of the flow topology in the post-shock region, so special attention should be paid to the modelling of variable-density STI.

Figure 26. The PDFs of the normalized magnitude of the different contributions from Lagrangian dynamics for (a) isotropic turbulence, (b) single-fluid post-shock turbulence and (c) multi-fluid post-shock turbulence.

To better understand the underlying mechanisms that cause the behaviour highlighted above, the dynamic equations (3.5c ) governing the vector ( $\text{D}Q/\text{D}t,\text{D}R/\text{D}t$ ) are examined. In figure 26, PDFs of the normalized magnitude of the different contributions from Lagrangian equations are shown to study the relative importance of different dynamics. The normalization used here for the vectors is the same as that used in figure 24. Figure 26(a) shows that for isotropic turbulence, the pressure Hessian term has the largest magnitude and the baroclinic contribution is the smallest. Mutual interaction and viscous terms have almost the same magnitude and distribution. After interacting with the shock wave, the magnitude of the baroclinic term is amplified for both single- and multi-fluid turbulence, but still remains the smallest comparing to the other contributions. The mutual interaction term becomes less important due to its reduced magnitude for both cases. The viscous term, however, exhibits different behaviour between single- and multi-fluid cases: it is amplified in the single-fluid case and reduced in the multi-fluid case. The pressure Hessian term is also amplified and remains the largest among all the terms. The percentages of contributions, using the means, indicate that the percentage of pressure Hessian contribution increases from 61.3 % to 74.9 % (single-fluid case) and to 73.9 % (multi-fluid case) across the shock wave.

Figure 27. Contributions to the transport equations of the VGT invariants by different terms for isotropic turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term and (d) viscous term.

The Lagrangian dynamics of the flow can be understood better by considering the CMVs of different terms in the $(Q,R)$ plane. As a reference, these terms are shown in figure 27 for isotropic turbulence. The variable $Q$ tends to be amplified in the enstrophy-production-dominated region due to the effects of the vortex stretching mechanism and is decreased in the dissipation-production-dominated region due to self-amplification of the strain rate tensor. On the other hand, the mutual effects on $R$ are small because the first invariant $P$ is usually small and the positive and negative values are likely to cancel each other. The contributions from the pressure Hessian (figure 27 b) tend to move the particles away from an asymptotic line, ending up amplifying the magnitude of $R$ . This result agrees well with that observed in turbulent boundary layers (Chu & Lu 2013). For the current simulation, the asymptotic line starts from $Q_{2}$ and ends in $Q_{4}$ with a slope of around $-2.5$ . The baroclinic contributions are very small in the post-shock turbulence as shown in figure 27(c). The viscous effects as shown in figure 27(d) and as expected reduce the magnitudes of $Q$ and $R$ and push the particles towards the origin. This has been observed in various types of turbulence (Ooi et al. 1999; Chu & Lu 2013). The combined effects from the four above mechanisms determine the circulating behaviour of the conditional mean of ( $\text{D}Q/\text{D}t,\text{D}R/\text{D}t$ ) vectors.

Figure 28. Contributions to the transport equations of the VGT invariants by different terms for single-fluid post-shock turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term and (d) viscous term.

After interaction with the shock wave, the CMVs in the ( $Q,R$ ) plane are different from those in the pre-shock isotropic turbulence. Figure 28 shows the results for the single-fluid case. By comparing it with figure 27, we note that even though the conditional means of ( $\text{D}Q/\text{D}t,\text{D}R/\text{D}t$ ) vectors are different, the contributions from mutual interaction (figure 28 a), baroclinic term (figure 28 c) and viscous term (figure 28 d) are very similar. The only term that is qualitatively different in post-shock turbulence and isotropic turbulence is the pressure Hessian term (figure 28 b). The importance of the pressure Hessian term is also reflected in the dynamical contributions in the $(Q,R)$ plane. In the post-shock single-fluid turbulence, the asymptotic line that separates the vectors into two regions with different behaviours disappears. Instead, the pressure Hessian term tends to move the particles away from the origin in $Q_{1}$ and $Q_{2}$ , thus increasing Q and $|R|$ values of the particles. In $Q_{3}$ and $Q_{4}$ , the vectors are parallel to the left-hand zero discriminant line, making the particles move from $Q_{3}$ to $Q_{4}$ , and then to $Q_{1}$ .

For multi-fluid post-shock turbulence, the pressure Hessian term is also the only term that is qualitatively different from that in isotropic turbulence (figure 29). Despite the increased density and pressure gradient in the multi-fluid case, the baroclinic term is still considerably smaller than all the other terms. In $Q_{2}$ and $Q_{3}$ , an asymptotic line similar to that in isotropic turbulence seems to exist, which ‘repels’ the vectors away from it, causing an increase in $|R|$ values. In $Q_{1}$ and $Q_{4}$ , the magnitude of the pressure Hessian term becomes much smaller. The further conditioned pressure Hessian term based on the local densities in figure 30 indicates that fluid particles with different densities have very different behaviours with respect to pressure Hessian dynamics. Specifically, the pressure Hessian generally moves the heavy-fluid particles towards the regions with larger $Q$ values. In $Q_{3}$ and $Q_{4}$ , it also moves the heavy-fluid particles towards the $R>0$ plane. For the light-fluid particles, the pressure Hessian term tends to make them move towards regions with larger $|R|$ values in the first and second quadrants. In $Q_{3}$ and $Q_{4}$ , the fluid particles move from $Q_{4}$ to $Q_{3}$ . Last but not least, the fluid particles with medium density seem to exhibit behaviour similar to that of light-fluid particles, except in $Q_{1}$ , where the pressure Hessian contribution is moving the fluid particles towards the regions with large $Q$ values. Examining figures 25 and 30 together, we observe that the differences in particle dynamics in the ( $Q,R$ ) plane in regions with different densities are mainly due to differences in the pressure Hessian contributions.

Figure 29. Contributions to the dynamics of the VGT invariants by different terms for multi-fluid post-shock turbulence. (a) Mutual interaction among invariants, (b) pressure Hessian term, (c) baroclinic term and (d) viscous term.

Figure 30. Contributions from pressure Hessian to the dynamics of the VGT invariants in (a) light-fluid region, (b) medium-density-fluid region and (c) heavy-fluid region.

4 Conclusions

Accurate shock-capturing turbulence-resolving simulations are used together with Eulerian and Lagrangian particle tracking post-processing methods to investigate the interaction of an isotropic turbulence with a normal shock wave for both a single fluid and a binary mixture of fluids (species) with different densities. The main objective is to develop a better understanding of the variable-density effects on the post-shock turbulence structure and its evolution away from the shock. Grid convergence tests are conducted to establish the numerical accuracy of the simulated data. The results show that the turbulence statistics are grid-converged, indicating good accuracy of the current computational method. Statistical convergence is also conducted for Lagrangian data.

The analysis is restricted here to an Atwood number of 0.28, based on the molar masses of the two fluids. At this Atwood number value, the variable-density effects introduce important changes in the turbulence structure, while the shock remains in the wrinkled regime for the shock Mach number considered. Similarly, the turbulent Mach number is also small enough that the multi-fluid case does not transition to the broken shock regime and the post-shock compressibility effects are weak. On the other hand, the Reynolds number is large enough so that the viscous effects stay small through the interaction with the shock and the corresponding single-fluid simulation satisfies the LIA limit. As the flow transitions to the broken shock regime due to larger turbulent Mach number and/or Atwood number, additional effects appear. These are left for future studies.

The density effects on the post-shock turbulence structure are first studied using Eulerian data. The different roles that the variable-density field could play through the shock interaction are identified for some of important statistics. The non-Gaussianity of the VGT is studied by examining the PDFs of velocity gradient components. The preferential amplification of rotation and strain rate tensors is found to be affected by the density variations, leading to a weaker correlation between the two tensors in the multi-fluid case. This is shown to be caused by different roles that density plays in the modification of rotation and strain rate tensors across the shock wave. The skewness and flatness of VGT components before and after the shock wave are then examined to study the evolution of the VGT. It is shown that density effects are weak across the shock, but are stronger in the post-shock development. The density variations are also shown to cause preferential alignment between eigenvectors of the strain rate tensor and density gradient vector, which then modifies the skewness of the velocity gradient and density gradient PDFs.

The density effects on the flow topology are then examined by first comparing the joint PDF of the second and third invariants of the deviatoric part of the VGT. The pre-shock joint PDF has the classic tear-drop shape. However, after the shock wave, a tendency towards symmetrization of the joint PDF, as in single-fluid STI, is observed for the multi-fluid case, with more data points falling into the first and third quadrants. After conditioning the joint PDFs based on fluid density, large differences among heavy-, medium- and light-fluid regions are observed. In the heavy-fluid regions, the joint PDF becomes almost completely symmetric with an increasing portion of data falling in the third quadrant. In contrast, the majority of the light-fluid data points have a similar distribution to that of isotropic turbulence. A connection between low vortex stretching and the joint $Q^{\ast }$ , $R^{\ast }$ statistics is established for the post-shock turbulence, by considering the contribution to vortex stretching rate from each quadrant.

Furthermore, Lagrangian fluid particles are used to track the development of the turbulence and VGT after the interaction with the shock. The Lagrangian dynamics of the VGT is also examined by using the conditional mean rate of change of the invariants of the VGT. For the parameter range considered, the results show that particles in $Q_{3}$ have the least residence times, while those in $Q_{2}$ have the longest residence times. The residence times are smaller than those in isotropic turbulence, especially in the multi-fluid case. It is also shown that particles starting in quadrants $Q_{1}$ and $Q_{3}$ play an important role in recovery of the vortex stretching term. After interacting with the shock wave, the ‘clockwise circulating’ behaviour (as observed in isotropic turbulence) disappears in both single- and multi-fluid cases. Our analysis highlights the mechanisms through which post-shock turbulence recovers the classical tear-drop shape, with an enlarging head in the second quadrant and elongating tail in the fourth quadrant. The contributions from different terms in the dynamic equations of VGT invariants, compared with isotropic turbulence, show that the pressure Hessian term is critical to the topological evolution of turbulence. The relative magnitude of the pressure Hessian term is increased and its dynamical contributions in the $(Q,R)$ plane are modified across the shock wave. The pressure Hessian term is also shown to be strongly dependent on the local density in the multi-fluid case, resulting in completely different dynamics in regions with different densities. In this work, the out-of-plane $(Q,R)$ compressibility effects are not considered due to the relatively low $M_{t}$ and $A_{t}$ . The compressibility effects and their coupling with the variable density effects will be investigated in more detail in future studies.