2.1 Initial condition creation and characterization

Figure 1.

Multi-mode initial condition creation. (a) Sample corrected multi-mode initial condition image where nitrogen is light and carbon dioxide is dark. (b) Simplified schematic showing the general method of interface creation by outflowing the mixture of gases at the interface, shown as red arrows, and flowing pure nitrogen and carbon dioxide below and above the interface, respectively, to create shear and buoyancy effects. In the schematic, the angle of the interface is slightly exaggerated for clarity. In both depictions gravity acts normal to the interface. The actual inlets and outlets consist of small perforations to prevent disturbance of the shock.

Figure 2.

IC characterization: (a) Sample multi-mode IC (image intensity corresponds to mole fraction of nitrogen). (b) Ensemble average of mole fraction for 200 multi-mode IC images. (c) Standard deviation of same 200 images. (d) Multi-mode IC after $10^{\circ }$ rotation with respect to inclined shock tube ( $80^{\circ }$ ). (e) Average energy spectra for increasing number of ensemble-averaged images for rotated IC. (f) Energy spectra calculated from a rotated IC in (d) within a concentration between 0.18 and 0.3. Different slopes are shown by red lines.

One of the main goals of this paper is to understand the role of the initial condition in shock-driven, variable-density flows. In order to achieve this goal, we have implemented two sets of initial conditions which vary in their modal composition. The first case is an inclined interface (tube angle of $80^{\circ }$ relative to horizon) which is a predominantly single-mode interface. The single-mode interface is created following the process outlined previously (McFarland et al. Reference McFarland, Reilly, Creel, McDonald, Finn and Ranjan2014b ; Reilly et al. Reference Reilly, McFarland, Mohaghar and Ranjan2015). It was determined that this method creates a relatively flat interface with a diffusion thickness ${\approx}0.96~\text{cm}$ . For insight into the mathematical representation of the inclined interface, which is effectively half of a period of a triangular wave, Mikaelian (Reference Mikaelian2016) suggested that the inclined interface can be closely represented by half of a cosine wave which has a Fourier coefficient carrying 81 % of the Fourier expansion. The multi-mode interface is created by injecting light gas into the heavy gas just below the interface creation plane and by injecting heavy gas into the light gas just above the interface creation plane. The light gas then rises from its injection point in the heavy gas section while the heavy gas falls from its injection point in the light gas section to create a multi-mode interface between the two larger sections of gas. This configuration, combined with the suction of a vacuum pump at the interface outflow slots, creates a counter flow condition with a combination of shear and buoyancy instabilities along the interface. A sample multi-mode initial condition PLIF image corrected for variations in laser intensity and background noise effects is shown in figure 1(a), where the intensity of the image concentration corresponds to the mole fraction of the light gas (white), and the multi-mode interface creation is described in figure 1(b). The initial condition characterization was done by taking 200 PLIF images of the multi-mode interface and then using energy spectra of density and Proper Orthogonal Decomposition (POD) to provide a statistical analysis of the interface.

Figure 3.

Dominant modes of the POD for 200 ensemble-averaged multi-mode initial conditions and their respective kinetic energy: (a) Mode 0 (Ensemble average), (b) Mode 1 (4.4 % of total kinetic energy), (c) Mode 2 (4.1 % of total kinetic energy). (d) Modal and cumulative fractional kinetic energy for the first 100 (90 %) modes. (e) Modal fraction of kinetic energy for all 200 modes.

Figures 2(b) and 2(c) show the ensemble average and standard deviation, respectively, of the nitrogen mole fraction for 200 PLIF images. As seen in figure 2(c), the largest fluctuation occurred near the bottom surface (the left side of the shock tube) where the mole fraction standard deviation is 1.2 %, which shows the repeatability of our initial condition. In order to isolate the effect of the secondary modes from the large-wavelength inclined perturbation on the energy spectra, images of the inclined multi-mode interface are rotated $10^{\circ }$ . One such rotated interface sample is shown in figure 2(d). In a method similar to Weber et al. (Reference Weber, Haehn, Oakley, Rothamer and Bonazza2014), the energy spectrum is calculated along the perturbed region within the concentration between 0.18 and 0.3 (small to mid-range scales), then averaged in the streamwise direction for each realization and then ensemble-averaged for the total number of realizations. Figure 2(e) shows the ensemble average of density energy spectra for an increasing number of interface images, indicating that the spectrum converges after approximately 20 images.

Using 100 images, the ensemble average of energy spectra in the perturbed region is shown in figure 2(f). Three ranges of different slopes can be distinguished from the analysis of the rotated image in figure 2(f). According to figure 2(f), the primary peak wave number is $k\approx 40~\text{m}^{-1}$ , and the spectral range that has the highest energy ( $k<40$ ) is related to the shock tube dimension. The spectra undergo a change in slope at $k\approx 200~\text{m}^{-1}$ , highlighting an intermediate range from $40<k<200$ . The scales in this range relate to shear and buoyancy effects in the initial condition. Finally, the smallest scales, which are related to the diffusion thickness, are seen in a third range ( $k>200$ ), which has the steepest slope of the spectra. Rather than a fit to the data, the figure shows fiducials for the three ranges.

To compute the POD of the initial condition, the same set of 200 PLIF realizations used in the density energy spectra investigation were utilized. In POD, each mode represents a coherent flow structure based on its reoccurrence in all images. Furthermore, to find the kinetic energy of a particular mode, one must look at the associated eigenvalue (Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993). In the POD analysis, mode 0 represents the ensemble average of the 200 density fields, as shown in figure 3(a). Additionally, when comparing figure 3(a) with figure 2(b), it can be shown that the images are similar, and thus indicates the accuracy of the POD code used (both of them show ensemble average of flow fields). Although the total, 199, modes in this decomposition were analyzed, to achieve better visualization only the first 90 % of the the modal fraction and cumulative fractional energy of modes are considered, as is shown in figure 3(d). Figures 3(b) and 3(c), respectively, represent modes 1 and 2 of the decomposition, containing 4.4 % and 4.1 % of the total energy, respectively. Figure 3(d) shows that the first 18 modes contain 50 % of the perturbations’ kinetic energy, which suggests multi-mode perturbations are occurring on the interface. To further clarify this point, the modal fraction of kinetic energy is shown in figure 3(e), which illustrates that the first 18 modes comprise a majority of the flow perturbations’ kinetic energy. The connection between POD and energy spectra of density can be observed by considering the mid-range scales in modes 1 and 2 of POD analysis which are approximately $2{-}5~\text{cm}$ and are in good agreement with the scales of the intermediate range of energy spectra.

4.1 Analysis of the density field

We will now probe the state of the flow by analyzing quantitative concentration fields (correction method described in § 2). Values of density calculated from 1D gas dynamics can be applied to the known values of concentration. From the density fields we can extract statistics that describe the state of mixing and begin to understand the flow field evolution.

The total width of the mixing region is measured as the distance in the direction of shock propagation between 5 % and 95 % values of the concentration profile (averaged in the direction normal to shock propagation). Figure 6(a) shows the mixed width versus time. The error bars represent plus or minus one standard deviation. The data points without error bars are from a previous study (Reilly Reference Reilly2015), with the same conditions as the single-mode case S presented here, but with only PLIF measurements. The data from the current study is consistent with the previous data as expected, but the multi-mode case displays a slightly larger mixed width at 2.64 ms (note that error bars overlap). The larger mixed width could be due to the inflow of gases slightly above and below the interface, as detailed in § 2, combined with the increased thickness of the mixing layer which is superposed on the inclined interface. The growth of the interface follows a fairly linear trend at the earliest times, as expected for an amplitude-to-wavelength ratio of the inclined perturbation $\unicode[STIX]{x1D6FC}=0.088$ , and maintains a trend which is fairly close to linear even until reshock. The linear growth corresponds to the strong penetration seen in the density fields. The overall growth rate of the mixing layer calculated from experiments can be compared to analytical results. The growth rate calculated from experiments is approximately $17.2~\text{m}~\text{s}^{-1}$ just after shock passage (initial growth rate is very similar for both cases), for which the linear theory of Richtmyer predicts a growth rate of $16.9~\text{m}~\text{s}^{-1}$ (Richtmyer Reference Richtmyer1960). Using the model of Weber, Cook & Bonazza (Reference Weber, Cook and Bonazza2013) to extract the growth rate by modeling vorticity deposition and applying the Biot-Savart law yields a growth rate of $17.3~\text{m}~\text{s}^{-1}$ . As secondary instabilities begin to have larger effects, a change in growth rate from linear to non-linear is expected as the penetrating mass spreads through the mixing region more rapidly and homogeneously. Although a strongly nonlinear trend is not obvious, it is apparent that secondary instabilities play a role in the mixing process; we see significant spreading of mass throughout the mixing region in the multi-mode case before and after reshock. Additionally, the post-reshock, single-mode case appears to break up due to secondary instabilities, although the linear perturbation at the largest scale is somewhat separated from the scales at which this mixing is occurring.

Figure 6.

(a) Total mixing width at different times for both single-mode (case 1) and multi-mode (case 2), and (b) total mixedness at different times for both single-mode (S) and multi-mode (M) cases.

Mixedness profiles are a measure of how well-entrained fluid is molecularly mixed. If none of the entrained fluid is mixed molecularly (immiscible case) the mixedness is zero, but if all the entrained fluid is molecularly mixed (a pure substance, locally homogeneous) then mixedness is one. The mixedness was calculated using $\unicode[STIX]{x1D703}=\overline{\unicode[STIX]{x1D709}_{1}\unicode[STIX]{x1D709}_{2}}/\overline{\unicode[STIX]{x1D709}_{1}}\times \overline{\unicode[STIX]{x1D709}_{2}}$ , where $\unicode[STIX]{x1D709}$ represents mole fraction and subscripts 1 and 2 refers to light and heavy gases.

Mixedness can be calculated by averaging over the entire mixing region, as shown in a temporal progression in figure 6(b). The time plot shows that from 2.64 ms to 5 ms the penetration of additional fluid as the interface stretches decreases mixedness due to gradient intensification. Case M shows higher mixedness at times before and just after reshock. After reshock, both cases undergo an increase in mixedness because the mixed material is confined to a smaller region during the inversion. From 6.4 ms to 9 ms the mixedness decreases again due to intensification of the gradients as the interface grows again. Although the mixedness is higher at all times for the multi-mode case, at 9 ms it is the closest between the cases. As the interface stretches, which causes additional entrainment but also increased turbulent mixing, the mixedness remains close to 0.8. Similar to the observations of Schilling et al. (Reference Schilling, Latini and Don2007), Lombardini, Pullin & Meiron (Reference Lombardini, Pullin and Meiron2012) and Tritschler et al. (Reference Tritschler, Olson, Lele, Hickel, Hu and Adams2014) the mixedness approches 0.85 asymptotically, as increased mass is injected into the mixing region and entrained fluid continues to become homogenized by turbulent mixing. It appears that the interface is approaching this point; however, given the extent of inhomogeneity (which remains throughout the mixing region at 9 ms), it is likely that the interface is still approaching the asymptotic state. Moreover, observation of figure 7 gives insight into the how well-mixed fluid is in different regions. Mixing in the bubble region (heavy gas side) is enhanced in the multi-mode case before reshock, while the bubble region of the single-mode case remains sharper before reshock. The effect of the rollup manifests itself in case S, at 5 ms, as a slight bump in the mixedness followed by a sharp upturn to one as as it approaches the pure nitrogen region. The multi-mode case is less mixed in the bubble region at 9 ms, which is evidence of continued turbulence development at smaller scales.

Figure 7.

Temporal evolution of mixedness for both (a) single-mode (S) and (b) multi-mode (M).

The autocorrelation of density fluctuations with specific volume fluctuations, also known as the density self-correlation (DSC) or the $b$ parameter as introduced by Besnard et al. (Reference Besnard, Harlow, Rauenzahn and Zemach1992) in the Besnard, Harlow, and Rauenzahn, or BHR model, helps describe instantaneous mixing occurring in the flow:

The averaging in (4.1) has typically been performed in two different ways by the RMI community: spanwise averaging, when there are not enough realizations to form an ensemble, and true ensemble averaging when there are enough realizations. The effect of averaging method on the value of $b$ was demonstrated by Balakumar et al. (Reference Balakumar, Orlicz, Ristorcelli, Balasubramanian, Prestridge and Tomkins2012). The DSC, averaged over the spanwise direction, is zero for pure fluids, including homogeneous mixtures, where there is no spatial variance in the density. Thus, the spanwise averaging method describes the level of inhomogeneity in the spanwise direction. If averaged over an ensemble, the DSC is zero where there are no fluctuations from the ensemble. This averaging method produces values for $b$ which correspond to the density fluctuation intensity, different information than mixedness provides. Livescu & Ristorcelli (Reference Livescu and Ristorcelli2008) shows the relationship between mixedness and the DSC, and if both are calculated from spanwise averaging they both describe homogeneity (inhomogeneity) in the flow. In the present work, ensemble averaging was used to produce a DSC field which was then averaged in the spanwise direction to obtain profiles along the streamwise direction. Because of the coupling between the density and velocity fields, density fluctuations are selected to the flux of mass (the velocity of the mass flux will be described in detail later) which drives the production of Reynolds stress in the BHR model. The production term of the Reynolds stress from the BHR model (Schwarzkopf et al. Reference Schwarzkopf, Livescu, Baltzer, Gore and Ristorcelli2016) is

where $a$ is the mass flux.

According to figure 8, the profiles of the DSC are similar at 2.64 ms because the interface has just been compressed by the incident shock. This compression confines the area over which the density can fluctuate while the interface stretching is sharpening the gradients. However, the mixing is slightly higher in the multi-mode case. At 5 ms the profiles are the most different, as the single-mode case shows localized peaks in mixing at shear-driven roll-ups, whereas the multi-mode case shows mixing throughout the mixing layer due to increased mode merging effects. Both cases show similarly high levels of mixing at 6.4 ms due to the recent forcing of the reshock wave. However, at 9 ms, the DSC is higher for case M in the spike tip and bubble regions, which is evidence that turbulent mixing mechanisms in those regions persist even after the forcing of reshock has ceased. The mixing intensity at 9 ms after reshock is commensurate with the asymptotic theory of mixedness in that as more fluid is injected into the mixing region, which is growing nonlinearly, density fluctuations act in synergy with strong velocity fluctuations, which describes the motivation for (4.2).

Figure 8.

Temporal evolution of the combined density self-correlation for both (a) single-mode (S) and (b) multi-mode (M) cases.

To conclude analysis of the density fields, and preface subsequent discussion on the velocity fields, the density gradients are analyzed in conjunction with vorticity. The PLIF resolution is not sufficient to show the exact scalar dissipation rate with $Sc\approx 1$ , the Batchelor scale, $\unicode[STIX]{x1D706}_{b}\approx \unicode[STIX]{x1D706}_{k}$ ( $\text{resolution}\approx 7\unicode[STIX]{x1D706}_{b}$ ); however, the dissipation rate filtered by the resolution scale is provided in appendix C. An alternative method of visualizing density gradients is a simulated schlieren technique (Quirk & Karni Reference Quirk and Karni1996) which is shown in figure 10 with the vorticity overlaid. This result from simultaneous measurements provides a description of how well the small-scale features of the flow (scalar and velocity gradients) are connected. The simulated schlieren field is derived from the density field by

where $\unicode[STIX]{x1D709}$ is mole fraction. Sharper gradients appear darker in the images. Most of the mixing region in case S has sharper gradients at each experimental time. The gradients appear sharper towards the spike side (light gas side, left) of the interface at all times as the pure heavy gas is injected into the light gas and the resultant shearing strips off mass as it mixes. The multi-mode case appears to break down to smaller scales much more quickly than the single-mode case, but both appear to have become much more disorganized at 9 ms. The vorticity is discussed in further detail in the subsequent section.

4.2 Analysis of velocity field

In this section, the acquired velocity fields are analyzed and discussed through the use of vorticity fields (pre- and post-reshock) and the effect of the velocity-averaging method on Reynolds stress profiles is discussed. Probability density functions (PDFs) of velocity fluctuations (Orlicz et al. Reference Orlicz, Balasubramanian, Vorobieff and Prestridge2015), and the positive, negative and total circulation are explained in appendix C. In order to reliably perform velocity field analysis, the velocity fluctuations must be calculated correctly, from spanwise or ensemble averaging. To calculate velocity fluctuations using spanwise averaging, the bulk velocity $u_{b,i}$ was defined as the streamwise and spanwise spatial average of the $i$ th component of the velocity. The relative velocity was then defined as the difference between the instantaneous and bulk velocities,

Next, the relative velocity was spatially averaged over the spanwise or $y$ -direction. Then the velocity fluctuation, $u^{\prime }$ , was defined as the relative velocity minus the spanwise-averaged velocity,

Figure 9.

Comparison of Favre-averaged Reynolds stress components in the multi-mode case at 9 ms based on two methods of velocity fluctuation calculation: (a) spatial versus (b) ensemble averaging.

Figure 10.

Vorticity evolution overlaid on simulated schlieren images for both cases. S refers to single-mode case and M refers to multi-mode case.

Velocity fluctuations found using ensemble averaging were obtained by taking the difference between the instantaneous velocity and an ensemble average. Equations (4.8) and (4.9) show the ensemble average, which is defined as the velocity at each $x,y$ location averaged over all realizations at each experimental time:

The two methods yield different results. To highlight the effect of ensemble averaging versus the use of traditional spanwise averaging, Favre-averaged Reynolds stresses for case M at 9 ms (late time, post-reshock) are shown in figure 9. For the Favre-averaged Reynolds stress, velocity fluctuation components come from mass-weighted average quantities: $\tilde{u_{i}}=\overline{\unicode[STIX]{x1D70C}u_{i}}/\overline{\unicode[STIX]{x1D70C}}$ . Figure 9 shows that the most significant difference between the two methods is the magnitude of the $u^{\prime \prime }u^{\prime \prime }$ term, which is significantly higher when using spanwise averaging, due to the large fluctuations in the streamwise velocity across the interface. The $v^{\prime \prime }v^{\prime \prime }$ and $u^{\prime \prime }v^{\prime \prime }$ terms are similar in magnitude but show more variation across the interface in the ensemble averaged case. We expect this variation would decrease as more experiments are executed. In a statistically repeatable, homogeneous turbulent flow field it would be expected that the ensemble and spatial averages would yield the same result. Comparatively, in our flow field the spanwise velocity is relatively insensitive to the averaging method, suggesting it is homogeneous in the spanwise direction. However, the large difference in the streamwise direction highlights the need for ensemble averaging due to the inhomogeneity of the streamwise velocity across the interface. Consequently, velocity statistics were calculated using ensemble averaging.

Resolving the vorticity field and flow circulation are two key aspects of a velocity field analysis. Using the PIV methods outlined earlier, the velocity fields were resolved using TSI INSIGHT 4G™ in order to provide a quantitative analysis of the flow field. Any holes in the velocity fluctuation field were filled via interpolation using nearby values. The field was subsequently filtered using a median filter with a bin of 11-by-11 vectors in order to smooth spurious noise induced by small gradients in the fields and allow for a more accurate analysis. The two methods that have been used to compute vorticity are explained in appendix C.

The evolution of the vorticity field is shown in figure 10, where red and blue represents positive and negative vorticity, respectively. Due to the misalignment of the pressure and density gradients, the baroclinic term results in vorticity deposition across the flow field, which in turn results in mixing and deformation of the interface. For the single-mode case, shown starting with figure 10(a), the incident shock wave passing from light to heavy gas creates strong negative vorticity across the interface. The magnitude of the vorticity immediately after incident shock compares well with the analytical model of Weber et al. (Reference Weber, Cook and Bonazza2013), which gives approximately $0.9\times 10^{4}~\text{s}^{-1}$ . At 5 ms after incident shock, the vorticity deposited on the interface has caused the initial perturbations to grow in amplitude and form larger vortical structures, as seen in figure 10(b). As the amplitude of the perturbation is larger than the IC and the density gradients are enhanced, it is expected that a larger amount of vorticity will be deposited as the reshock wave passes back from the heavy to light gas regions, leading to much higher levels of vorticity (Weber et al. Reference Weber, Haehn, Oakley, Rothamer and Bonazza2012). This can be seen in figure 10(c), where as the reshock wave passes through the interface, it deposits such a large amount of vorticity that the majority of the field becomes positive. This reversal of the vorticity field is a result of the application of additional opposite-signed (now positive) baroclinic torque, and the interface itself begins to rotate in the opposite direction immediately after reshock. Furthermore, small pockets of negative vorticity are preserved even after reshock. The presence of positive and negative vorticities after reshock, seen in figure 10(d), is a result of vorticity deposition on an interface with a much less unidirectional gradient. It is also apparent that, due to the more mixed state of the reshock initial condition (5 ms), there is a larger area of the flow over which vorticity is deposited.

For the multi-mode case, shown as a progression starting from figure 10(e), the general physics is similar but key differences should be highlighted. In figure 10(e), immediately following the incident shock wave, the complex nature of the multi-mode interface has already resulted in the deposition of positive vorticity even before reshock. In figure 10(f), the vorticity has grown again, feeding the growth of the perturbations that serve as the initial conditions for the multi-mode reshock case. Figure 10(g) shows the vorticity field immediately prior to reshock and displays the anticipated deposition of positive vorticity due to the application of baroclinic torque. As expected, the deposition of vorticity in the multi-mode case results in higher-amplitude positive vorticity post-reshock than in the single-mode case, due to the presence of positive vorticity pre-reshock. However, this figure displays a notable difference between the two cases, as although the amount of vorticity deposition is equivalent, its distribution is markedly different. In the single-mode case, the vorticity deposition occurs along the interface, as it closely follows its shape. Whereas in the multi-mode case the distribution of vorticity deposited is broader and more diffuse. Lastly, it can again be seen from figure 10(h) that by late time the vorticity becomes more widely deposited throughout the flow field, its amplitude grows, and the level of mixing increases.

4.3 Analysis of simultaneous velocity and density field statistics

Figure 11 shows velocity vector fields which have been registered to density fields as discussed in § 2. This visualization shows the connection between the velocity field and the density field, and highlights that there is a velocity gradient across the density gradients, especially at the edge of the mixing layer where the density gradient is strong. However, a more chaotic state is present in the bubble region, which is generally more mixed, but still contains a level of local inhomogeneity.

The alignment of vorticity with concentration gradients is just one example of the utility gained by registering the simultaneous density and velocity fields. Density-velocity statistics are critical to validation of turbulence models for variable-density flows and give insight into the physics of multi-component mixing in the developing RMI flow.

Figure 11.

Registered velocity and concentration fields for 9 ms multi-mode and single-mode inclined cases. (a,b) Include $1/4$ of the velocity field resolution and indicate the regions to be zoomed in (c,d), which show $1/2$ of the velocity field resolution for clarity.

The mass flux due to turbulent fluctuations found by correlating the density and velocity fluctuations, is $a_{i}=\overline{\unicode[STIX]{x1D70C}^{\prime }u_{i}^{\prime }}/\overline{\unicode[STIX]{x1D70C}}$ . The overbar represents an ensemble average which yields a field of values. To obtain a profile, the fields are averaged in the spanwise direction.

The streamwise velocity mass flux is shown in figure 12. Based on the choice of axes, a negative value denotes that, on average, positive fluctuations in density correspond to a velocity fluctuation in the streamwise-direction toward the region of pure light gas while negative fluctuations in density correspond to a velocity fluctuation in the opposite direction. This leads to turbulent interpenetration.

Figure 12.

Temporal evolution of the combined turbulent mass flux in the streamwise direction of the velocity fluctuations. Single-mode (S) and multi-mode (M) cases.

It is obvious that at all times case M has more fluctuating mass flux across the majority of the mixing layer. At the earliest time, 2.64 ms, the two cases appear very similar in shape, but with slightly more flux in the multi-mode case. At 5 ms, two small negative peaks can be seen in the single-mode case where shear-driven rollups cause increased fluctuations, while the multi-mode case has a maximum near the centre of mass where POD analysis of the initial condition suggests there is significant modal content. Reshock amplifies the flux in both cases as more energy is delivered to the mixing region. At 6.4 ms, the single-mode case shows peaks at the edges of the mixing region where there is intense penetration. The multi-mode case shows a similar level of flux near the edges, but a maximum near the centre of mass. At the latest time (9 ms) the bubble region (the right side of the image) is becoming more homogenized, and therefore has less mass flux than the spike region, where a sharper contrast in density yields more intense flux of mass. Similar to 6.4 ms, the multi-mode case has greater mass flux over almost the entire domain, reinforcing the assumption that the multi-mode case is in a state of more active turbulent development.

4.4 Reynolds stresses and anisotropy

In this section we examine the anisotropy of the flow field through the Reynolds stresses, quantified by the normalized anisotropy tensor:

The components of the tensor are bounded by $-1/3$ to $2/3$ , where a value of $-1/3$ indicates the component has no energy while a value of $2/3$ indicates the component contains all of the energy. A value of zero is approached as the Reynolds stresses become equal, indicating an isotropic turbulent energy. For our case, only the streamwise and spanwise components of velocity are available and so we assume that the fluctuations in the out-of-plane direction are equal to those of the spanwise direction. This sets $\unicode[STIX]{x1D6FD}_{22}$ equal to $\unicode[STIX]{x1D6FD}_{33}$ . This assumption is likely less accurate before reshock, but is a necessary limitation of 2D PIV measurements.

Figure 13 shows the components of the spanwise average of the normalized anisotropy tensor for each case at two times, immediately before reshock (5 ms), and late after reshock (9 ms). From this figure it is evident that the anisotropy decreases after reshock for both cases, with large regions of nearly equal components near the centre of mass in the flow. Before reshock, both cases show a peak in anisotropy near the centre of mass, though the single-mode case shows a higher anisotropy. Both cases at 9 ms show strong non-zero values for the $\unicode[STIX]{x1D6FD}_{11}$ and $\unicode[STIX]{x1D6FD}_{22}$ components in the upstream flow field, indicating that the former, pre-reshock bubble region now contains the majority of anisotropy. In the downstream region of the post reshock flows we see a much more isotropic region, with the $\unicode[STIX]{x1D6FD}_{22}$ term taking a positive value.

Figure 13.

The normalized anisotropy tensor components at (a,b) 5 and (c,d) 9 ms for both (a,c) single-mode and (b,d) multi-mode cases.

Additional insight into the flow can be gained by examining the spatial distribution of the anisotropy components. Figure 14 shows the spatial distribution of $\unicode[STIX]{x1D6FD}_{11}$ . It is noticeable that the anisotropy does not drop to zero away from the interface, due to a small amount of noise which manifests predominantly in the streamwise direction. Although the magnitude of some of the velocity fluctuations away from the interface is small compared to the fluctuations at the interface, the normalization of anisotropy does not typically weight the magnitude of fluctuations because it is classically used in homogeneous flows. Nevertheless, it is customary to present information about inhomogeneous flow in this context while understanding the caveats involved. Before reshock, $\unicode[STIX]{x1D6FD}_{11}$ is highest upstream of the interface in the bubble-front region and near the spike tip. The single-mode inclined interface case shows a less well organized field due in part to its less turbulent flow field. After reshock, these regions persist and are transported by the overturning process. The downstream region now contains the high $\unicode[STIX]{x1D6FD}_{11}$ found in the pre-reshock spike while the anisotropy of the pre-reshock bubble front is pushed upstream by the emerging post-reshock spike structure. In order to consider the trend of anisotropy over time, where much lower anisotropy for the multi-mode case before reshock can be seen more clearly, the spanwise average of normalized anisotropy tensor magnitude is shown and explained in appendix C.

Figure 14.

$\unicode[STIX]{x1D6FD}_{11}(x,y)$ for each case at (a,b) 5 ms (before reshock) and (b,d) 9 ms (late time post reshock). (a,c) Single-mode and (b,d) multi-mode cases.

5.1 Turbulent length scales and Reynolds numbers

The Taylor microscale can be defined using two methods. The first method is based on fitting a quadratic function to the peak and two points, both on the right side of the peak, of the velocity or density fluctuation autocorrelations. The slope of the resulting curvature at $f(0)$ is 0 (Pope Reference Pope2000; Ramaprabhu & Andrews Reference Ramaprabhu and Andrews2004). The average of the correlation functions for each field can be defined as

where $F$ is density or velocity fluctuation (streamwise or spanwise), and $\langle F(x,y)F(x+r,y)\rangle$ is the spatial average (over $x$ and $y$ ) and the ensemble average over all the fields. Correlation functions can also be used to calculate the integral scale by finding the area under the curve of $f(r)$ ,

Using correlation functions to calculate Taylor microscales and integral scales is useful for estimating the inertial range in energy spectra. Consider an example using the longitudinal correlation function of streamwise velocity fluctuations at 9 ms, as is shown for the single-mode and multi-mode cases in figures 16(a) and 16(b), respectively. To find the slope of the energy spectra, one must calculate the least-squares fit of the points between the Taylor micro and integral scales in the energy spectra. The error bars in figure 16 represent plus or minus one standard deviation. Inspection of figure 16 suggests that the separation between micro and integral scales, which is also an indication of the inertial range in the energy spectra, will increase after reshock. Furthermore, at 6.4 ms the highest $L/\unicode[STIX]{x1D706}=23/2.4\approx 10$ is observed in the single-mode case.

Figure 16.

Correlation functions of streamwise velocity fluctuation for both cases at 9 ms (a,b) and the temporal evolution of the corresponding Taylor microscales (c) and integral scales (d) calculated from the correlation functions.

An alternative method for calculating the Taylor microscale is by computing the variance and gradient of the velocity or density fluctuations (Pope Reference Pope2000). Thus the Taylor microscale can be calculated as

where $\unicode[STIX]{x1D702}$ is the velocity or density field. According to figure 17, both methods show similar trends for Taylor microscales given by the velocity fluctuations, as both increase after reshock. Figure 17 also shows that the Taylor microscales from streamwise and spanwise velocity fluctuations are higher in the multi-mode case compared to the single-mode case due to higher fluctuations. Taylor microscales derived from streamwise density fluctuations calculated using the spatial gradient method, are also shown in figure 18. This figure shows a trend similar to the streamwise velocity fluctuation Taylor microscales from the correlation function method. However, the magnitude of Taylor microscales based on velocity fluctuation is between 1.5 and 3 mm, notably different from those given by density fluctuations, which are around 0.4 to 1 mm. This difference is thought to be due to resolution effects and will be discussed later by comparing Batchelor and Taylor microscales from different Reynolds number calculations.

Figure 17.

Comparison between the use of spatial gradients and correlation functions to find the Taylor microscales derived from (a) streamwise and (b) spanwise velocity fluctuations.

Figure 18.

Temporal evolution of the streamwise density fluctuation based Taylor microscale, found using spatial gradients.

The Taylor microscales computed above were then used as a comparison for those found using the Reynolds number. When using the Reynolds number to characterize turbulence, several calculation methods must be considered because of the unsteady nature of the flow and the disagreement in the literature about which method is best for characterizing transitioning, variable-density flow. Therefore, three different methods are considered to determine the relative applicability and insight they can provide for this particular experiment.

The first method uses the total mixing width growth rate, $Re_{h}=h{\dot{h}}/\unicode[STIX]{x1D708}_{mix}$ , where $h$ is the total mixing layer width, ${\dot{h}}$ is the mixing width growth rate, and $\unicode[STIX]{x1D708}_{mix}$ is the mixed kinematic viscosity of the two gases. The mixed kinematic viscosity can be found using $(\unicode[STIX]{x1D707}_{1}+\unicode[STIX]{x1D707}_{2})/(\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2})$ , where 1 and 2 denote the light and heavy gases, respectively. The mixing width and the mixing width growth rate data were extracted from figure 6(a), which was previously explained in the analysis of density statistics. The mixing width growth rate was then determined by fitting a linear trend line to the mixing width data shown in figure 6(a). Due to the changes imparted by the reshock and overturning of the interface, the mixing width undergoes a drastic change. For this reason, two different values are also found for the pre- and post-reshock mixed kinematic viscosity by using one-dimensional gas dynamics and property relations; these values are approximated to be $6.27\times 10^{-6}$ and $4.08\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , respectively.

The circulation Reynolds number was the second method investigated, $Re_{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D6E4}/\unicode[STIX]{x1D708}_{mix}$ , where $\unicode[STIX]{x1D6E4}$ is the absolute value of the average total circulation, calculated using the methods outlined in § 4.2.

The final method, the turbulent kinetic energy based Reynolds number is $Re_{K}=\sqrt{\overline{K}}h/\unicode[STIX]{x1D708}_{mix}$ , where the turbulent kinetic energy is $K(x,y)=(u^{\prime }(x,y)^{2}+v^{\prime }(x,y)^{2})/2$ , and $\overline{K}$ is the spatial average of the $K(x,y)$ field on the whole domain.

Figure 19 shows the time dependence of each Reynolds number for both single-mode and multi-mode cases. All three methods show a similar trend, as they increase with time at approximately the same rate. This similar trend was expected for $Re_{h}$ and $Re_{\unicode[STIX]{x1D6E4}}$ because they have similar origins (Jacobs & Sheeley Reference Jacobs and Sheeley1996). Furthermore, Tomkins et al. (Reference Tomkins, Balakumar, Orlicz, Prestridge and Ristorcelli2013) found similar agreement, within a factor of two, between the two methods. Comparatively, $Re_{K}$ shows lower magnitudes than the other two, yet displays the same overall trend. As stated by Dimotakis (Reference Dimotakis2000), achieving $Re\approx 10^{4}$ is conjectured to be necessary to observe a mixing transition and is a more restrictive condition beyond conditions for turbulent transition. Employing this criterion showed that the Reynolds numbers obtained via the above-mentioned methods satisfy the necessary condition for transition to turbulence at all times.

Figure 19.

Temporal Reynolds number values found using three different methods: Mixing width growth rate, total circulation, and total kinetic energy.

Further information about the energy spectra can be extracted through the Taylor microscale and Kolmogrov scale. These scales can be estimated with a simplifying assumption of local isotropy and knowledge of the Reynolds number values. The Kolmogrov length scale is calculated from the largest length scale in the flow, the mixing width $h$ , and the Reynolds number using $\unicode[STIX]{x1D706}_{k}\approx hRe^{-3/4}$ . Similarly, the second useful length scale, the Taylor microscale, is $\unicode[STIX]{x1D706}_{T}\approx \sqrt{15}hRe^{-1/2}$ (Pope Reference Pope2000; Tomkins et al. Reference Tomkins, Balakumar, Orlicz, Prestridge and Ristorcelli2013; Orlicz et al. Reference Orlicz, Balasubramanian, Vorobieff and Prestridge2015).

With information provided by the Kolmogrov scale, one can better estimate in which scales the viscosity is dominant at higher wave numbers. The Taylor microscale gives insight into the inertial sub-range, and allows one to delineate the maximum size of scales which can be assumed to be isotropic and free of influence from large-scale flow features (Robey et al. Reference Robey, Zhou, Buckingham, Keiter, Remington and Drake2003).

By comparing the Taylor microscales derived from Reynolds numbers to those based upon correlation functions, the effects of resolution can be investigated. According to figures 17 and 20, the Taylor microscales measured from velocity correlation functions and spatial gradients are larger by a factor of 3–4 compared to those of $Re_{h}$ and $Re_{\unicode[STIX]{x1D6E4}}$ , and by a factor of 1.5–2 compared to that of $Re_{K}$ . The differences between these methods implies that the resolution effects in velocity could be responsible for the differences. However, comparing Taylor microscales given from the Reynolds number and density correlation function methods (figure 18) shows almost identical results, indicating that the PLIF resolution is fine enough to capture Taylor microscales from density fields.

Other important scales which have been used to show turbulent transition criteria are the Liepmann–Taylor scale ( $\unicode[STIX]{x1D706}_{L}\approx 5hRe^{-1/2}$ ) and the inner viscous scale ( $\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D708}}\approx 50\unicode[STIX]{x1D706}_{k}$ ). The threshold for transition, and therefore the possible presence of an inertial range in the spectra, based on these scales is ( $\unicode[STIX]{x1D706}_{L}/\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D708}}>1$ ), where both scales can be calculated from Reynolds numbers (Dimotakis Reference Dimotakis2000).

Figure 20.

Taylor microscale (a) and Kolomogrov scale (b) progressions found with Reynolds numbers derived using three different methods.

Figure 21.

Progression of the ratio of Liepmann–Taylor scale to the inner viscous scale provided by Reynolds number values from three different methods.

Values of $\unicode[STIX]{x1D706}_{L}/\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D708}}$ computed from the three different Reynolds numbers methods are shown in figure 21. For each of the three calculations, since the fraction is above 1, an inertial range is present, especially after reshock. However, work by Lombardini et al. (Reference Lombardini, Pullin and Meiron2012) suggests that for fully developed turbulence $\unicode[STIX]{x1D706}_{L}/\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D708}}$ should be above 10, which is not achieved in the present case. Therefore, based on all discussed above (Reynolds number calculations, the fraction of integral scale to correlation function Taylor microscale, and the fraction of Liepmann–Taylor scale to the inner viscous scale), one can conclude that while a mixing transition is occurring in the present case, especially after reshock, fully developed turbulence is likely not reached.

5.2 Time evolution of energy spectra and structure functions

In order to understand the growth of different scales and to quantify the inertial range, the time evolution of fluctuating kinetic energy, fluctuating enstrophy and the density variance spectra before and after reshock are presented. The general procedure used to compute the spectrum of a quantity $\unicode[STIX]{x1D711}(x,y,t)$ (velocity, density or enstrophy fluctuation) is described by Reilly et al. (Reference Reilly, McFarland, Mohaghar and Ranjan2015). Thus, the one dimensional energy spectrum of any quantity is

where $FFT(\unicode[STIX]{x1D711}^{\prime })$ is the Fourier coefficient of the $2\unicode[STIX]{x03C0}$ periodic function $\unicode[STIX]{x1D711}(x,y,t)$ in the spanwise direction and $\text{conj}(FFT(\unicode[STIX]{x1D711}^{\prime }))$ is its complex conjugate.

The temporal development of the total fluctuating kinetic energy spectra is computed by adding the $u^{\prime }$ and $v^{\prime }$ components of the energy spectra for the $i$ th column:

Figure 22.

Temporal evolution of the combined total fluctuation kinetic energy spectra (a,b), a comparison of compensated spectra at early and late reshock times (c,d), and the normalized cross-correlation spectrum (dashed line fitted with a solid line) (e,f).

The progression of the kinetic energy spectra without $w^{\prime }$ component and a comparison of the compensated spectra for the two experimental cases at each time is shown in figure 22. In figure 22(a,b), one can observe oscillations before reshock, especially in the single-mode case at early time (2.64 ms), which are then smoothed by reshock. The spectra increase after reshock (6.4 ms) for both cases, as the second shock adds more energy to the mixing region. Conversely, energy decays after incident shock at 5 ms and at late time after reshock (9 ms). Compensated spectra are also shown in figure 22(c,d), where the compensation exponent, $n$ , is shown in the legend. Because the out-of-plane component of velocity is unknown, the spectra is compared to that of the inertial sub-range scaling of both ( $-5/3$ ) and ( $-3$ ) predicted by 3D and 2D classical turbulence, respectively. It is important to note the assumptions of classical spectra are that the flow is both isotropic and incompressible. The flow is fairly incompressible after passage of shock, but is anisotropic and inhomogeneous. However, similar to Schilling et al. (Reference Schilling, Latini and Don2007) and Orlicz et al. (Reference Orlicz, Balasubramanian, Vorobieff and Prestridge2015), the scalings associated with classical turbulence of expected power laws of $5/3$ , for the small, energy-containing wave numbers, and 3, for larger wave numbers are used to compensate the velocity fluctuation spectra. Once compensated, a full decade ( $k_{y}\approx 200$ –2000) of the inertial range approximately follows the $-5/3$ law after reshock, at 6.4 and 9 ms. However, less than half a decade of spectra in the large wave numbers is observed to approximately follow the $k^{-3}$ power law. This is likely an effect of resolution limitation and not a physical dissipation range. The compensated energy spectra is used to emphasize the separation between large and small scales which is another indication of transition to turbulence (Zhou Reference Zhou2017).

In order to investigate the level of anisotropy and memory of initial conditions across different scales with time for both single and multi-mode cases, the normalized cross-correlation spectrum between $u^{\prime }$ and $v^{\prime }$ , $|E_{u^{\prime }v^{\prime }}(k_{y})|/[E_{u^{\prime }u^{\prime }}(k_{y})E_{v^{\prime }v^{\prime }}(k_{y})]^{1/2}$ , is considered in figure 22(e,f) (Tavoularis & Corrsin Reference Tavoularis and Corrsin1981a ,Reference Tavoularis and Corrsin b ; Shen & Warhaft Reference Shen and Warhaft2000). Anisotropy is clearly present before reshock (2.64 and 5 ms) in the single-mode case even at smaller scales; however, in the multi-mode case, the tendency toward isotropy of small scales is evident as the wave number increases. This local isotropy in the multi-mode case also suggests that memory of the initial condition is maintained before reshock (both 2.64 and 5 ms) because imposed perturbations in the IC causes a tendency toward isotropy in mid-range scales. After reshock (6.4 and 9 ms), the cross-correlation spectrum decays more rapidly in both single- and multi-mode cases, which indicates that local isotropy is reached at lower wave numbers than before reshock. The decrease in normalized cross-correlation spectrum from 6.4 ms to 9 ms in both cases at larger scales is an indicator of breakdown of large scales to small scales between 6.4 and 9 ms. Comparison of single- and multi-mode cases at 6.4 and 9 ms shows the similarity between imposed scales in the initial condition and the scales which result from breakdown of large to small scales upon reshock. Thus, a significant range of scales exists which are locally isotropic in both cases, and the anisotropy is primarily contained in the large scale which is common to the two cases.

The fluctuating enstrophy is defined as $(\unicode[STIX]{x1D714}^{\prime })^{2}/2$ and its time evolution is seen in figure 23(a,b), which also shows the compensated spectra for both of the cases after reshock (figure 23(c,d)), similar to the velocity energy spectra in figure 22. While the slope of fluctuating enstrophy spectra is lower than velocity spectra, it is similar to the velocity energy spectra in that it is again amplified by reshock. An additional similarity is that the enstrophy spectra also decays at late time after the incident shock (5 ms).

Figure 23.

Temporal evolution of the fluctuation enstrophy spectra (a,b) and a comparison of compensated spectra at early and late reshock times (c,d).

The fluctuating enstrophy spectra in homogeneous, isotropic flow (classical turbulence) can be related to velocity spectra via $E_{\unicode[STIX]{x1D714}^{\prime }\unicode[STIX]{x1D714}^{\prime }}(k,t)=k^{2}E(k,t)$ . The expected slope of the inertial range in the fluctuating enstrophy spectra will be $1/3$ at small wave numbers and $-1$ at large wave numbers for 3D turbulence. Moreover, according to Batchelor (Reference Batchelor1969) and Kraichnan (Reference Kraichnan1967), the energy spectra of classical 2D turbulence scales with $k^{-3}$ , which yields $-$ 1 scaling for enstrophy spectra at inertial range scales. Therefore, $-$ 1 scaling in the enstrophy spectrum can be compared with both classical 3D turbulence at larger wave numbers and classical 2D turbulence at smaller wave numbers (Kraichnan Reference Kraichnan1967; Batchelor Reference Batchelor1969). Figure 23(c,d) shows approximately $k^{-1}$ scaling in the compensated spectra for less than one decade after reshock; however, the largest scales show scaling close to $k^{1/3}$ for almost one decade ( $k_{y}\approx 20$ –200).

Structure functions were calculated in order to find the inertial ranges and slopes of spectra so they could be compared with those from different spectra. The second-order longitudinal structure function is $S_{2}(r)=\langle (\unicode[STIX]{x1D702}(x+r)-\unicode[STIX]{x1D702}(x))^{2}\rangle$ , where $\unicode[STIX]{x1D702}$ is the velocity or density fluctuation in the direction of the separation length, $r$ , and brackets $\langle \rangle$ denote both spatial and ensemble averaging.

Figure 24.

Temporal evolution of the velocity structure function (a,b) and a comparison of compensated structure function at early and late reshock times (c,d).

Figure 24 shows the time evolution of the velocity structure functions for each case (a,b), as well as the comparison between the compensated structure functions after reshock (c,d). The decay of the small-scale contribution can be seen during both pre- (2.64–5 ms) and post-reshock (6.4–9 ms) for the single- and multi-mode cases. In addition, the large scale peak is similarly dampened at 5 and 9 ms for both cases. Analyzing the compensated structure functions after reshock shows a slope close to that of the $2/3$ power law of Kolomogorov turbulence (classical turbulence) for almost an entire decade.

The temporal development of the density fluctuation spectrum and the comparison of compensated spectra for each case after reshock are shown in figure 25(a,b) and (c,d), respectively. Unlike the velocity fluctuation spectra, there is a small increase in the density spectra after incident shock, between 2.64 and 5 ms, and after reshock, between 6.4 and 9 ms, especially at small scales. Yet, similar to the velocity spectra, reshock sharply increased the spectra magnitude, which suggest that density fluctuations are stronger post-reshock. Classical scaling of $-5/3$ is again of interest, even for density spectra. Therefore, the assumption of isotropy at small scales should hold true in order to have the same scaling as a classical spectra. Some researchers have found similar scaling for density spectra even with inhomogeneous and unsteady flow, which is not assumed in the classical case (Vorobieff, Rightley & Benjamin Reference Vorobieff, Rightley and Benjamin1998). In the current case, the density fluctuation spectrum shows an approximate $k^{-5/3}$ scaling for half a decade at intermediate wave numbers after reshock, as shown by the compensated spectra. There is a small increase in the inertial range which can be seen from 6.4 ms (early reshock) to 9 ms (late time after reshock). Furthermore, these results are verified by density fluctuation structure functions which are shown in figure 26. Post-reshock, density fluctuation structure functions also show close to $2/3$ scaling for intermediate $r$ and for less than a decade.

Figure 25.

Temporal evolution of the density spectra (a,b) and a comparison of compensated spectra at early and late reshock times (c,d).

Figure 26.

Temporal evolution of the density structure function (a,b) and a comparison of compensated structure function at early and late reshock times (c,d).