Effect of adiabatic index on Richtmyer–Meshkov flows induced by strong shocks

Richtmyer–Meshkov instability (RMI) is a phenomenon in fluid mechanics that describes the evolution of an interface between two fluids of distinct acoustic impedance and distinct densities when a shock wave impacts the interface. The flow evolution is illustrated by Fig. 1, where, after the shock passage, the light fluid (red) and the heavy fluid (blue) and the fluid interface (light green) travel as whole unit, and RMI develops. (Richtmyer, 1960; Meshkov, 1972). RMI may develop when the interface between the fluids is given an initial perturbation of magnitude a₀ and wavelength λ (seen on the right in Fig. 1) the interface amplitude increases over time with growth rate v₀ and the interface evolves into a large-scale coherent structure of bubbles and spikes Fig. 1) (Abarzhi, 2010; Abarzhi, 2008).

RMI appears in a variety of processes in high energy density plasmas, controlling material transformation under the high strain rate, governing the formation of hot spots in inertial confinement fusion, determining energy transport in a core-collapse supernova, and strongly influencing the evolution of blast waves and explosions (Meshkov, 1972; Richtmyer, 1960). RMI forms in situations characterized by strong shocks, sharply and quickly changing flow fields, and by small effects of dissipation and diffusion (Abarzhi, 2010). Interaction of a shock wave with a density discontinuity may result in the development of RMI and in extensive interfacial mixing (Meshkov, 1972; Richtmyer, 1960). Since RMI plays such a large role in these applications, the ability to understand and control this instability is important.

To the best of the author's knowledge, previous numerical simulations investigating the dynamics of RMI flows have only contained ideal monotonic gases. This study aims to investigate RMI flows for gases with more than one atom per molecule, e.g., diatomic or triatomic gases such as $O2$ or $H2O$ or more complex molecules such as $SF6$⁠. This is significant for investigation, as it can enable RMI studies to have a wider field of application to more accurately model situations with high-speed non-monatomic gases, and to aid in control of the instability through the initial parameter set-up. The applications include rocket thrust flow (such as those in scramjets), inertial confinement fusion, and explosions of supernovae (Drake, 2009; Bodner et al., 1998; Zel'dovich, 1967).

Having started from the seminal works of Robert Richtmyer and Evgeny Meshkov (Richtmyer, 1960; Meshkov, 1972, Abarzhi and Sreenivasan, 2020), RMI remains the subject of active research in fundamentals and in applications (Abarzhi et al., 2019). Recent developments include studies of strong-shock driven RMI using particle methods (Huang et al., 2020; Sun et al., 2020; Dell et al., 2017; [Dell et al., 2015; Stanic et al., 2012); numerical modeling of the effect of turbulence, presuming it may develop, on Richtmyer–Meshkov mixing (Thornber et al., 2019); studies of RMI with magnetohydrodynamics (Li et al., 2018); simulations of Richtmyer–Meshkov dynamics with re-shocks and in complex geometry (Li et al., 2019; Guan et al., 2020). This indicates a need in a systematic study of the effect of the adiabatic index on RMI evolution in order to provide reliable benchmarks for experiments, and to better understand a broad range of processes in nature and technology to which RMI is relevant (Abarzhi et al., 2019; Abarzhi, 2010).

The novelty of our work is in the systematic and scrupulous study of the effect of the adiabatic index on the evolution and the growth rate of RMI driven by strong shocks in fluids with contrast densities with intermediate and large amplitudes of the initial perturbations of the interface, and in the quantification of the properties of Richtmyer–Meshkov dynamics in a broad range of parameters, including the shock strength, the density ratio of the fluids, the amplitude of the initial perturbation, and the fluids' adiabatic indexes.

Sometimes, in applications such as in scramjets, plasma thursters, or inertial confinement fusion, the effects of the RMI and the RMI evolution are necessary to mitigate and control (Lindl et al., 2004). Some methods of achieving this include controlling RMI through adjusting the initial parameters of the system (Anisimov et al., 2013; Abarzhi, 2010; Demskoĭ et al., 2006). To do this, information is required about the effect the initial parameters have on the development of the instability. Of particular interest in controlling the development of the instability is the amount of energy available for interfacial mixing deposited into the interface by the shock wave. This is one of the aims of this study.

The volume of physical experimental data of RMI produced by strong shocks is sparse as the experiments have tight and challenging to achieve requirements on the flow control, implementation, and diagnostics. (Motl et al., 2009; Orlicz et al., 2009; Jacobs and Krivets, 2005). Therefore, numerical modeling of RMI is a powerful tool to aid in designing and building systems in which RMI is present. However, the dynamics of RMI are complex and a numerical model should be able to manage numerous competing requirements, such as shock capturing, interface tracking, and accurate accounting for the dissipation processes (Stanic et al., 2012; McFarland et al., 2011; Herrmann et al., 2008; Dimonte et al., 2004).

Here, by applying the Smoothed Particle Hydrodynamics simulations and the theory, we systematically study the hydrodynamic problem of RMI evolution driven by strong shocks in a broad parameter regime. We focus on the effect of the adiabatic index on the fraction of energy available for interfacial mixing in RM flow. To do this, we will obtain data on two quantities of the flow- the interface speed and the interface growth rate. These quantities inform us to how much mixing of the interface occurs, and how well the numerical simulations capture small-scale structures, which the simulations must do well in order to obtain data on the interfacial mixing. The results of each of these quantities are compared whenever possible with rigorous theoretical theories, finding good quantitative and qualitative agreement.

RMI develops when a shock impacts an interface between two fluids of differing densities and the energy is distributed throughout the fluids (Aleshin et al., 1990). This paper will only focus on two dimensional RMI case, and with the shock propagating from the light fluid into the heavy fluid. When the shock hits an ideally planar interface (without any perturbation), it splits into a reflected shock traveling back through the light fluid and a transmitted shock traveling through the heavy fluid (Stanic et al., 2012; Herrmann et al., 2008; Demskoĭ et al., 2006). The bulk of the fluid influenced by the shock impact (the bulk between the reflected shock and the transmitted shock, including the fluid interface) all moves with the same velocity which we call the background velocity and which has magnitude $v∞$⁠, as seen in Fig. 1.

The speed $v∞$ quantifies the amount of energy transferred into the fluid bulk by the shock and is a function of the shock strength and fluid properties (Stanic et al., 2012; Richtmyer, 1960).

If the interface between the two fluids is perturbed (no longer planar), the bulk fluid containing the interface still moves with velocity $v∞$⁠, yet now the amplitude of the interface perturbation increases with the growth-rate v₀ due to the impulsive acceleration induced by the shock (Meshkov, 1972; Richtmyer, 1960). In Fig. 1, arrows mark the direction of fluid motion at the tip of the bubble (right) and spike (left). Eventually, the bubble and spike growth decelerates, and small-scale structures appear on the sides of the developing spikes and in the bulk (Stanic et al., 2012).

The Mach number, denoted M, is defined as the ratio of the speed of a shock to the speed of sound in the light fluid c_l (Richtmyer, 1960). For weak shocks, $M≈1$⁠, the background motion (the motion of the whole body of fluid) is subsonic relative to the light fluid, $v∞/cl<1$⁠, where c_l is the speed of sound in the light fluid. For $M≈3$ the ratio is $v∞/cl≈1$⁠. For shocks with $M≈5$ the background motion is supersonic, $v∞/cl>1$⁠, and for shocks with M > 7, the background motion may be hypersonic, $v∞/cl≫1$ (Stanic et al., 2013; Stanic et al., 2012; Herrmann et al., 2008). In this work, we consider Mach numbers 3 and 5, which are usually challenging to model accurately in the simulations.

The Atwood number, denoted A, describes the density difference between two adjacent fluids with a common interface:

where $ρh,ρl$ are the densities of the heavy and light fluids, respectively. In this work we consider Atwood numbers 0.6, 0.8, and 0.95, the latter of which is a very extreme density contrast of 40 to 1.

The adiabatic index of a substance, denoted γ, can be understood from three perspectives. From a thermodynamic point of view, the adiabatic index gives an important relation for an adiabatic process of an ideal gas:

where P, V, and γ is the pressure, volume, and adiabatic index, respectively, of the fluid.

It can also be understood as the ratio of the heat capacity at constant pressure C_P to the heat capacity at constant volume C_V,

From a molecular dynamics point of view, the adiabatic index can also be related to the degrees of freedom f of a molecule as

To the best of the author's knowledge, most of the simulation analyses of the dynamics of RMI flows have been conducted with an adiabatic index of $γ=5/3$⁠, which is the value for ideal monotonic gases. Its value decreases for gases with more than one atom per molecule, e.g., for diatomic gases $γ=7/5$⁠. While γ is known to have a strong influence on the flow dynamics, no systematic study has been undertaken on the effect of γ on the dynamics. The gases analyzed in this study are theoretical ones, where we vary the adiabatic index of the heavy and light fluids systematically instead of choosing particular gases.

The initial amplitude of the perturbed interface is the important parameter of RM dynamics, since for steady shocks RMI develops only when the interface is initially perturbed (Dell et al., 2017; Stanic et al., 2012). In this work, we consider the single-mode sinusoidal initial perturbations of the interface with amplitude a₀ and wavelength λ. For given λ, we systematically vary the ratio $a0/λ$ from 0 (initially flat interface) to 1 (interface with large initial perturbation) (Dell et al., 2017; Dell et al., 2015; Stanic et al., 2012). This systematic variation is needed because even in ideal fluids with adiabatic indexes $γh=γl=5/3$ and with given values of the perturbation wavelength λ, the Mach number M and the Atwood number A, the initial growth rate of the RMI v₀ has a complex dependence on the initial perturbation amplitude a₀ (Dell et al., 2017; Velikovich et al., 2014; Nishihara et al., 2010; Wouchuk, 2001; Richtmyer, 1960).

After the incident shock completely passes the perturbed interface, the RM interfacial growth begins, and the structure of bubbles and spikes appears. The RMI growth rate is a function of time and of the initial perturbation amplitude, its accurate quantification is a challenge. Initially (i.e., immediately after the shock passage), the amplitude of the interface perturbation grows linearly with time, and the RMI initial growth rate v₀ can be defined as the slope of the growth of the amplitude of RM unstable interface. For very small initial amplitudes, $a0/λ≤10−2$⁠, the initial growth rate v₀ increases linearly with the initial amplitude a₀. For relatively small initial amplitudes, $(a0/λ)≈10−2$ to $10−1$⁠, the initial growth rate increases with $a0/λ$ slower than linear. For relatively large initial amplitudes, $(a0/λ)≈O(10−1)$ the initial growth rate v₀ achieves a maximum value. For very large initial amplitudes $(a0/λ)≈1$ the initial growth rate of RMI is negligible (Dell et al., 2017).

The recent study (Dell et al., 2017) finds that for ideal fluids with adiabatic indexes $γh=γl=5/3$ the maximum initial growth rate of RMI is achieved for $(a0/λ) ∼0.3$ for a broad range of the Mach and Atwood numbers, in agreement with experiments (Dell et al., 2017). In this work we systematically investigate the effect of the fluid's adiabatic indexes on the RMI initial growth rate, the value of the maximum initial growth rate, the initial amplitude at which the maximum growth rate is achieved, and the amount of energy available for interfacial mixing over a broad range of the Mach and Atwood numbers, the fluids' adiabatic indexes and amplitudes of the initial perturbation.

The parameter regime we investigate is for Mach numbers 3 and 5, and for Atwood numbers 0.6, 0.8, and 0.95. We vary the adiabatic index of the heavy and light fluid $(γl,γh)=(1.2,1.3,…,1.6)$⁠. In order to explore the effects of the perturbation of the interface on the dynamics, we vary the amplitude of the initial perturbation of the interface from 0% to 100% of the interface wavelength for each pair of γ_l and γ_h, i.e., $a0/λ=(0,0.1,0.2,…,1)$ where a₀ is the amplitude of the sinisoidal initial perturbation, and λ is the wavelength of the perturbation. As discussed before, there is a well-studied regime for $γl=γh=5/3$⁠, documented well in the past, allowing us to compare our results with previous studies (Dell et al., 2015; Stanic et al., 2012; Stanic et al., 2013).

Our parameter regime is very broad and gives us a good physics picture of the behavior of RM dynamics. Our parameter regime contains 1650 cases. We ran a numerical simulation for each case. The average simulation takes 36 h to run, making a total of about 60 000 h. The simulations were run on three Windows laptop computers with i7 processors and with 8 GB, 12 GB, and 16 GB of RAM.

In this work, we compare the results from our numerical simulations to analytical solutions. This analysis takes different forms depending on the progression of the instability. In the initial linear regime of RMI, the interface perturbation grows at a constant rate v₀, which is a function of the amplitude a₀ and wavelength λ of the initial perturbation (Richtmyer, 1960). In the following nonlinear regime, the interface perturbation growth rate decreases, and a large coherent structure of spikes and bubbles appears (Abarzhi, 2010; Abarzhi, 2008). The heavy fluid penetrates the light fluid in spikes as seen on the bottom right of Fig. 1. As they travel, the spikes decelerate and small-scale structures form on the sides of the spikes and in the bulk (Stanic et al., 2012).

For ideal fluids, RM dynamics is governed by the conservation of mass, momentum, and energy as

with the spatial coordinates x_i, time t, and the fields of density ρ, velocity v, pressure P, and energy density $E=ρ(e+v2)/2$⁠, where e is the specific internal energy (Abarzhi, 2010). These equations are augmented with the standard Rankine–Hugoniot conditions defining the post-shock flow fields, and with the boundary value problem, including the boundary conditions at the interface (i.e., the continuity of normal component of velocity and pressure at the interface), and the conditions at outside boundaries of the domain. The initial conditions prescribe the flow fields at the interface and in the bulk. Rigorous theories consider RM dynamics in certain limiting cases and find diagnostic benchmarks for RMI (Richtmyer, 1960; Nishihara et al., 2010; Abarzhi, 2010; Dell et al., 2017).

An important parameter of RMI dynamics is $v∞$⁠, the magnitude of the velocity (hereafter: velocity) of the bulk, or the background motion. This value can be precisely calculated by zeroth-order theory from the conditions of the conservation of mass, momentum, and energy, and the equations of state of the fluids (Richtmyer, 1960). For ideal gases, it is a function of the initial shock's Mach number, the adiabatic index of the fluids $γh(l)$⁠, and the Atwood number, $v∞=v∞(M,A,γh(l))$⁠, and it is useful because it quantifies the amount of energy deposited by the shock into the fluid bulk (Stanic et al., 2012). The analysis of the numerical simulations becomes much simpler once $v∞$ is obtained and is used to set the characteristic timescale. It also defines the inertial frame of reference moving at speed $v∞$ (Dell et al., 2015).

Another important parameter in RMI dynamics is the initial growth rate v₀ of the interface. It is a function of M, A, $γh(l)$⁠, and the initial perturbation amplitude and wavelength, $v0=v0(M,A,γh(l),a0,λ)$ (Stanic et al., 2012; Nishihara et al., 2010; Holmes et al., 1999). For very small amplitude (⁠$a0/λ=10−2$ or smaller), v₀ is precisely calculated by linear theory, and grows linearly with a₀, $v0≈a0/λ(Mcl)$ (Nishihara et al., 2010; Wouchuk, 2001; Richtmyer, 1960).

For moderately small values of a₀, the growth rate v₀ becomes non-linear and is calculated in weakly nonlinear theories Velikovich et al., 2014 and Nishihara et al., 2010. For larger values of a₀, the rate v₀ may grow with a₀ even slower than the linear and weakly nonlinear theory predict: it achieves a pick value for some initial amplitude and then exponentially decays (Stanic et al., 2012; Dell et al., 2015).

The late stages of Richtmyer–Meshkov dynamics have been calculated with a group theory approach (Abarzhi, 2010; Abarzhi et al., 2003; Abarzhi, 2002). At this late stage, the bubbles decelerate and flatten, and extensive interfacial mixing occurs (Stanic et al., 2013; Stanic et al., 2012; Herrmann et al., 2008).

Numerical modeling of RMI in the extreme conditions of strong shocks, contrast fluid densities and large initial perturbation amplitudes is a difficult task because the method should be able to accurately handle large speeds, strong shocks, and preserve small-scale structures with high precision and accuracy. These small-scale structures are embedded in large-scale dynamics, moving at high speeds, so the order of precision required is very large (Abarzhi et al., 2019; Anisimov et al., 2013). To model these complex dynamics we have employed the Smoothed Particle Hydrodynamics code (SPHC), which is an open-source code developed by Dr. Stellingwerf; it has free access to a complete set of validation test cases (Stellingwerf, 1991). This code has been widely used and tested on a broad variety of shock and flow problems; particularly it has been used by NASA to investigate the Space Shuttle Columbia incident (Stellingwerf et al., 2004).

SPH is a grid-free method that represents a continuous fluid with fixed-mass SPH particles, which are each represented by a mathematical basis function (or kernel) (Stellingwerf, 1990). In essence, SPHC keeps track of a large array of particles and for each time step, and calculates interactions between all particles. This enables an accessible treatment of the governing equations – which are the partial differential equations representing the conservation of mass, momentum, and energy (Stanic et al., 2012, Stellingwerf, 1990). Similarly to the Eulerian finite-difference methods approximating the flow fields at a finite number of grid points, SPH applies the particle and the kernel approximations to reduce the governing PDEs to a set of ordinary differential equations. In particle approximation in the SPH, a fluid quantity is presented as $f(x)=∫Ωf(x′)δ(x′−x)dx′,$ where x is the position vectors of the particles, Ω is the integration volume and δ is the Dirac delta-function. The Dirac delta function is replaced by a kernel W, such that

where h is a smoothing length and K is a coefficient determining the support domain boundary. The kernels are usually symmetric cubic B-spline functions resembling a Gaussian shape. This leads to $⟨f(x)⟩=∫Ωf(x′)W(x′−x)dx′$ and $⟨∇f(x)⟩=∫Ωf(x′)∇W(x′−x)dx′$ and upon discretization to

where $⟨…⟩$ represents the kernel approximation operator, and m_j, ρ_j are the mass and the mass density of the j particle. These expressions for the fluid quantities are substituted to the governing equations, and their numerical solution is sought. The SPH approach fully employs the advantages of Lagrangian method and is relatively easy to implement (Dell et al., 2015; Stanic et al., 2012; Monaghan, 2005; Stellingwerf, 1990).

Here we employ the Smoothed Particle Hydrodynamics code (SPHC) which handles well the SPH technique and is suitable for high-performance simulations, and which is diverse, proven, and robust (Stellingwerf, 1990). SPHC was successfully tested and validated on the strong shock-driven RMI, the Noh problem, and other shock and flow problems (Dell et al., 2017; Stanic et al., 2012; Stellingwerf, 1990). It was also applied for modeling of reactive fluids, material ablation in hypersonic flows, electrodynamics problems, and gravitational problems (Stellingwerf, 1990).

In this paper, the SPHC is employed in a hydrodynamic approximation. The SPHC consists of a standard SPH technique with Monaghan viscosity to handle shocks (Monaghan, 2005; Monaghan, 1988) and with Balsara weighting to reduce shear forces (Balsara, 1995). SPHC conserves mass, energy, momentum, as well as angular momentum globally and locally. In the SPHC, the interface is tracked by the particle ordering and is free from the needs for sub-grid reconstruction and grid refinement.

In the simulations, the particle size is varied in regions of varying density to maintain contact between adjacent particles, leading to an automatic increase in the resolution in the regions of higher density, e.g., near the interface, and to track the interface with high resolution. A small amount of conservative smoothing is applied to thermal energy to avoid artificial effects caused by rapid shock heating. This correction is substantially weaker than the “wall-heating” thermal conduction correction (Monaghan and Gingold, 1983]. The time integration in SPHC uses a Runge–Kutta algorithm (Benz, 1988; Benz, 1990) and automatically provides for fluid quantities the estimate of a computational error due to a finite time interval. In case the error exceeds a specified level, the step is recomputed with a smaller time interval; otherwise, the time step gradually increases for computational efficiency. In the SPHC simulations, the boundaries are consist of so-called “reflect” conditions which produce correct solutions near the boundary (Stanic et al., 2012).

In this study, the computational setup is similar to that in (Stanic et al., 2012) in order to extend and compare with the results of previous studies.

The RMI SPHC simulation is a two-dimensional array of particles in three regions, as shown in Fig. 1. The region on the left consists of the low-density gas with density ρ_l and the region on the right consists of the high-density gas with density ρ_h. The region in the middle is the interface between the two fluids. Post-shock, after the shock refracts the interface and splits into transmitted and reflected shocks, the regions of the light gas and the heavy gas (between the reflected and the transmitted shocks) and the interface move in the direction of the transmitted shock with the velocity of the background motion. The top, bottom, and right edge boundaries of the simulation are “reflect” boundaries, where particles bounce off elastically (Stanic et al., 2013; Stanic et al., 2012). The orientation of our RMI simulations is horizontal, where the shock moves from left to right. However, in order to be consistent with the theoretical analysis, we define the horizontal and vertical axis as the z- and x-axis, respectively. The interface is given a perturbation which is sinusoidal in nature, defined by the equation $z*=a0 sin (kx)$ where $k=2π/λ$ is the wave-vector, z and x are the Cartesian coordinates in the directions along with and normal to the shock direction, respectively, and t is time. The number of particles in each simulation is $2×105$⁠.

The initial pressure and temperature is $P0=2494×103$ Pa and $T0=3×102$ K, respectively. The gas constants $Rh(l)$ and densities $ρh(l)$ are $Rl=8.314$ J/kmol, $ρl=1×10−3kg/m3$⁠, and $Rh=2.079,9.238×10−1,2.132×10−1$ J/kmol, $ρh=4×10−3,9×10−3,3.9×10−2$⁠. Their energy per proton is high, ranging from $≈192$ eV to $≈1.03$ MeV initially, and reaching values from $≈1.1$ keV to $≈2.9$ MeV after the shock passage (Dell et al., 2015).

Our simulations are run in the two-dimensional domain $1·10−2m×3·10−2$ m. The initial perturbation wavelength is $λ=3.3×10−3$ m. The Mach number is defined relative to the speed of sound of the light fluid with $cl=2.039×103$ m/s (Pandian et al., 2017).

The characteristic length scale in our simulations is set by the wavelength of the initial perturbation λ. The characteristic velocity scale is chosen as $v∞$ (Stanic et al., 2012; Herrmann et al., 2008). The characteristic timescale is given by $τ=λ/v∞$⁠.

For each simulation, the total run time is $3τ$⁠, which ranges from $2.5×10−6$ s for $M=5,A=0.8$ to $1.9×10−6$ s for $M=10,A=0.95$⁠. This run time is sufficient for the reflected shock to reach the end of the computational domain, and it is sufficient to study the initial growth rate of RMI (Dell et al., 2015).

Figure 2 represents the diagnostics of the interface evolution caused by the shock-interface interaction (Richtmyer, 1960; Nishihara et al., 2010; Dell et al., 2017). The initial amplitude of the interface perturbation is set at the start of every run by the initial conditions. The shock propagates from the light to the heavy fluid and hits the interface. Then, the shock refracts the interface and splits on the transmitted shock and reflected shock, whereas the interface is compressed and its amplitude decreases. The RMI starts to develop only after the shock completely passes the interface, the transmitted and reflected shocks are fully formed at both sides of the interface, and the interface amplitude begins to increase. The shock-interface interaction is a complex process, and the accurate quantitative description of this process is challenging. Qualitatively, both the “impulsive” acceleration of the interface by the shock and the shock-induced vorticity deposition at the interface lead to the development of RMI and define the value of RMI initial growth rate v₀ (Nishihara et al., 2010; Wouchuk, 2001).

In our simulations, we want to obtain the initial growth rate of the RM unstable interface. Hence, as seen in Fig. 2 representing the interface evolution, we locate the origin $(0,0)$ at the first minimum of the interface amplitude, such that for t > 0 the amplitude increases from a minimum value of the already compressed amplitude achieved at the corresponding discrete instance of time. We then derive the RMI initial growth rate v₀ for t > 0 over some time interval, and we take v₀ to be the slope of the linear regression line from the first minimum amplitude over the next data points. Note that v₀ is defined as the time derivative of the difference of the initial positions of the bubble and spike, which is twice the amplitude, as in Stanic et al., 2012.

We note that an issue arises when choosing what data points to include, because the value of v₀ may vary depending on where the cutoff for the initial growth is defined. To solve this problem, we notice that the shape of the amplitude-time curve for early time changes depending on the value of the initial perturbation, marked in orange in Fig. 2. For $a0/λ=0.05$ the curve is concave, for $a0/λ=0.15$ the curve is linear, and for $a0/λ=0.5$ the curve is convex, forming a bubble like shape. Each of these shapes has clearly defined endpoints, marked on the diagram by the first blue points: the first finishes at the inflection point (the transition from concave to convex), the second where the growth starts to become nonlinear, and the third finishes at the end of the convex “bubble” (and where the second bubble begins). These endpoints for all three shapes all finish at $0.8τ$⁠, where $τ=v∞/λ$⁠, giving $0.8τ=6.67×10−7$ s. These patterns are consistent across all different values of $γh(l)$ and $a0/λ$ considered, and this method produces good results, which indicates that it is a good choice. While the time-interval $0.8τ$ which we choose in the present work is bigger than the time-interval $0.3τ$ chosen in our earlier works, in the appropriate cases the results in the paper are consistent the earlier results (Dell et al., 2017; Dell et al., 2015; Stanic et al., 2012).

To find $v∞$⁠, we measure the speed of the interface in the ideally planar interface case. The position of the interface is taken to be the leftmost position of the interface, as in Dell et al., 2015. The planar interface position motion is almost perfectly linear in time, allowing very accurate calculation of $v∞$ by taking the slope of the linear regression line of the interface position over time. When the interface is not planar, the interface may interfere with the dynamics of the bulk. In our study, the values of $v∞$ calculated for the planar case are taken to be the same for cases where $a0/λ≠0$ (Dell et al., 2015).

Using the SPHC simulations and the method described above, we calculated the value of $v∞$ for the cases $M={3,5},A={0.6,0.8,0.95}$⁠, with γ_h and γ_l ranging from $1.2,1.3,…,$ 1.6. We look at each pair of Mach and Atwood numbers and investigate how the background motion changes with the adiabatic index. Choosing two particular pairs of the Mach and Atwood numbers as an example in Fig. 3, we see that the speed of the background motion is negatively correlated with the adiabatic indexes, sometimes in a non-linear fashion (right) and sometimes in a linear fashion (left).

The speed of sound in fluids is given by the Newton–Laplace equation

where c is the speed of sound, K_s is a coefficient of stiffness, the isentropic bulk modulus, and ρ is the density. In an adiabatic process, $Ks=γP$ where P is the pressure and γ is the adiabatic index. In an adiabatic process, we also have

where V is volume, the constant is $a=P0ρ0−γ$⁠, and $P0,ρ0$ are some “initial” pressure and density of the gas.

Putting these together, we have

Figure 4 displays the values of $v∞$ scaled by the speed of sound in the light fluid (c_l) over our regime of γ_h and γ_l. Since the values of $γh(l)$ vary within 10%–30% in our analysis, $γh(l)∈[1.2,1.6]$⁠, then across all cases the dependence of the quantity $v∞$ on $γh(l)$ is approximately linear. Hence, in order to analyze how the Mach and Atwood numbers affect the speed of the background motion, we apply a plane of best fit using least squares regression. We also find that the average sum of the residuals in the fit is insignificant. The speed of the background motion increases as the adiabatic indexes of both fluids decrease, i.e., the post collision fluid moves faster when the molecules are more complex. More complex molecules tend to be stiffer, which allows for shock to be transmitted quicker.

Table I represents the dependence of $v∞/cl$ on $γh(l)$ across each Mach and Atwood number. The function $f∞=v∞/cl$ is given for variables $γ̃h(l)=5/3−γh(l)$⁠, in order to better see the variation in $v∞/cl$⁠, as we move away from the adiabatic index's well-studied value of $γh(l)=5/3$⁠. We see from these results that the dependence of $v∞/cl$ on $γ̃h(l)$ and hence $γh(l)$ is non-symmetric. This may owe to the direction of propagation of the initial shock from the light to the heavy fluid.

Finding the planes of best fit for each pair of Mach and Atwood numbers, we can investigate the effect these have on the rate of change of $v∞$⁠. Figure 5 displays the rate of change of $v∞$ with respect to γ_l (left) and γ_h (middle), as well as the value of $v∞/cl$ when $γl,γh=0$ (right), which correspond with the values of a, b, and c, respectively, in the equation of a plane $z=ax+by+c$⁠.

We find that as the Atwood number increases, the rate of change of $v∞$ over γ_h increases in an apparently linear fashion, getting closer to 0 (Fig. 5, middle) From our equation above, $Ks=γP$⁠, we know that the adiabatic index is proportional to the stiffness of a fluid, and so a possible interpretation of this result is that at higher density contrasts, the stiffness of the impacted fluid has less of an effect on the shock speed than at lower density contrasts. This makes intuitive sense since at higher densities, the fluid will be under greater pressure and will already be quite stiff, and hence less able to be affected by the adiabatic index. We also notice that the initial shock speed (the Mach number) has little effect upon the rate of change of $v∞$ over γ_h, which is in stark contrast with γ_l, where the Mach number has a large impact on the rate of change of $v∞$ over γ_l (Fig. 5, left). A possible explanation for this is that a low density, non-stiff fluid traveling at high speed is able to lose a more of its energy on impact due to absorption as the speed compresses it than the same fluid traveling at a lower speed.

From the plot of $v∞/cl$ against the Atwood number (Fig. 5, right), we observe that the shock speed decreases as the Atwood number increases, which may be because as the density contrast increases, the lighter fluid is unable to move the increasingly heavier fluid as easily.

We compare these results to the analytical calculations produced by zero-order theory (Richtmyer, 1960), where $v∞$ can be calculated precisely from zero-order theory in the planar interface case (when the initial perturbation amplitude is 0) (Stanic et al., 2012; Richtmyer, 1960). The percentage error between the simulation results and the zero-order theory for the case Mach 5, Atwood 0.8 is shown in Table II.

From this table we can see that when γ_h and γ_l are close, simulation agreement with the linear series is close, <7%. Values with errors above these thresholds are marked in bold. In particular, when $γl=γh$⁠, the agreement is excellent, $>99%$⁠.

We found that if γ_l and γ_h are too different with $γh(l)≥γl(h)+0.2$ for $γl(h)=(1.2,1.3)$ or $γh(l)≥γl(h)+0.3$ for $γl(h)=(1.4,1.5,1.6)$⁠, the simulations don't handle those situations as easily. To the best of the author's knowledge, this observation has not been made before, and may prove useful in the understanding of scenarios with varying adiabatic indexes.

From Table III we see that the average errors for each Mach and Atwood pair are very good, but increase as the Mach and Atwood numbers increase. This makes sense because higher Mach and Atwood numbers create more and more extreme environments with extreme densities, speeds, and shocks. This is harder for simulations to model accurately.

Table IV shows the average value of $v∞$ scaled by c_l, the speed of sound in the light fluid, and the Mach number M across all our Mach and Atwood numbers. By dividing the interface speed by the shock speed, this gives us a representation of the proportion of energy deposited into the interface. We see that the velocity of the background motion is only a fraction of the shock velocity, ranging from 20% to 50%.

After the shock passes through the interface, the interface amplitude grows approximately linearly with time, as we see in Fig. 2. When this value is calculated from the simulations, we scale it by the velocity of the bulk fluid, $v∞$⁠, which is done to compare the interface growth rate to the background motion. This value $v0/v∞$ quantifies the distribution of the energy imparted by the shock wave between the interfacial fluid and the bulk fluid (Dell et al., 2015).

We ran simulations for $M=3,5,A=0.6,0.8,0.95,γl,γh={1.2,1.3,…,1.6}$⁠, and $a0/λ={0,0.1,0.2…,1}$⁠. Determining the interface amplitude growth rate v₀ and bulk velocity $v∞$ from the simulations using the method as described in Sec. II C 2, we plot their ratio $v0/v∞$ for each value of $a0/λ$⁠, displayed in Fig. 6. A well-defined shape is formed, which increases linearly for early time, becomes nonlinear and eventually peaks and decreases for late time, asymptotically approaching 0. As in Abarzhi et al., 2019; Dell et al., 2017, and Dell et al., 2015, a function satisfying these criteria is

This depends on the Atwood number and two constants, C₁ and C₂. Our objective is to ascertain how well the simulation data fits this curve and to compare these constants with linear theory to assess the accuracy of our simulations.

Letting $x=a0λ$ and $y=v0v∞·1A$⁠, we have

As we see in Fig. 6, the data fits the large hump of the equation, and we would expect values of $a0/λ>0.5$ to fit the tail, which is confirmed in previous studies (Abarzhi et al., 2019; Dell et al., 2017, and Dell et al., 2015).

To determine the constants C₁ and C₂ of the equation of the curve and to quantify the strength of the relationship between the data and the curve, we rearrange Eq. (10) to form a linear relationship so that a linear regression line can be fitted to the data:

Letting $ŷ=ln (yx),m=−C2,n=ln (C1)$⁠, we have the linear relation

We plot $ŷ$ against x for all chosen values of $γl,γh$ (excluding those found to produce high errors) and for $a0/λ=(0.1,0.2,…,1)$⁠. The plot of all the values of γ_h and γ_l combined is shown in Fig. 7, where the relationship can be seen to be strongly negative and linear. The value of C₁ is calculated by taking the exponential of the y-intercept of the linear regression line, m, and C₂ is the negative of the slope of the line, n, which follows from the setup of Eq. (10).

Table V represents the dependence of the parameters $C1(2)$ on $γh(l)$ from the planes of best fit across each Mach and Atwood number. The plane function $f1(2)$ is given for variables $γ̃h(l)=5/3−γh(l)$⁠, in order to better see the variation in $v∞/cl$⁠, as we move away from the adiabatic index's well studied value of $γh(l)=5/3$⁠. We see that the dependence of $C1(2)$ on $γ̃h(l)$ and hence $γh(l)$ is non-symmetric. This may be associated with the direction of propagation of the initial shock. Figure 8 displays all the values of C₁ and C₂ across our entire regime of adiabatic indexes, Atwood and Mach numbers. We notice that each plot increases in value as the adiabatic indexes increase. (Dell et al., 2017) found that the values of C₁ and C₂ are approximately constant across all Atwood numbers for $γl,γh=5/3$⁠. This value is just outside of the space that we have analyzed, but extrapolating from our values using the equations in Table V we agree with this conclusion, with an average error of 9.53% which is good for extrapolation. However when $γh,γl$ are further away from 5/3, we find that the values of C₁ and C₂ are not as consistent. The data points are still clustered around the same region, but we find that for Atwood numbers of 0.95, C₁ actually decreases as γ_l increases for both Mach numbers, which is opposite to every other case.

In order to analyze some trends in the data and to make the data set more manageable, we take the average of the three data points for each Atwood number, since for the most part they are reasonably close for a particular Mach number. To better understand a two-dimensional image of a three-dimensional plot, we display the data in a graph with the adiabatic indexes arranged diagonally from smallest to largest in order to get a clearer picture of the data (Fig. 9).

From this graph, we can observe that both constants are higher for Mach 3. From our equation above, we know that the interface growth rate is a non-monotone function of the initial perturbation amplitude.

Observing this we see that a higher value of C₁ means a faster interfacial growth rate as it is directly proportional to v₀, and a higher value of C₂ means the maximum growth rate is achieved at smaller perturbations as well as a slower interfacial growth rate, because it increases the decay of the exponential.

To find the maximum value of the interfacial growth rate, we let $x=a0λ$ and $y=v0v∞·1A$⁠.

Substituting that into our original equation, we obtain

Note that in the case when C₁ and C₂ both change at the same rate, multiplying C₁ and C₂ by the same constant has no effect, and v₀ does not change.

Figure 10 displays the initial interface perturbation size $(a0/λ)$ at which the maximum interfacial growth rate is reached, and how this varies as we change the adiabatic indexes. We find that the perturbation decreases as the adiabatic indexes increases, across all Atwood numbers, in full consistency with results in Fig. 9.

Physically, this means that stiffer fluids spread out more quickly with smaller perturbations. The result may be of interest in experiments and applications that want to reduce the effects of RMI, as they need to be more careful with any perturbations between their stiffer fluids.

Figure 11 shows how the maximum interfacial growth rate changes as the adiabatic indexes change. We see that for Atwood numbers 0.6 and 0.8, the maximum growth rate gradually increases as the adiabatic indexes increase, but stay relatively constant for Atwood number 0.95. In consistency with the results presented in Fig. 9, when the values of C₁ and C₂ both change, the parameter C₁ is stronger influenced by the adiabatic indexes when compared to C₂ for Atwood numbers 0.6 and 0.8. We also notice that Atwood 0.95 cases have a much higher interfacial growth rate than the other two Atwood numbers. A possible explanation for this is that higher density contrasts lead to a higher compression of the interface during the initial impact because the lighter fluid is unable to move the heavier fluid as easily, which could lead to a more rapid “rebound” effect as the interface develops. From our simulation results presented in Figs. 10 and 11, the experiments can learn that stiffer fluids spread out faster when the density contrasts are not extremely high.

Of interest in experiments is the amount of energy that can be deposited into the interface by the shock, which determines the amount of energy available for interfacial mixing, discussed in Sec. IV. The maximum scaled interface growth rate $v0/v∞$ quantifies this amount.

Our results agree with experiments (Glendinning et al., 2003; Dimonte et al., 1996; Dimonte and Remington, 1993; Aleshin et al., 1990). According to our results, for any shock strength, density ratio, adiabatic indexes, and initial perturbation wavelength, there is a characteristic amplitude scale $[a0/λ]max$ at which the maximum growth rate is achieved. The existence of the amplitude scale in RMI suggests a possible mechanism of saturation of the growth of the mode with the wavevector $k=2π/λ$ at some initial amplitudes, which was observed in experiments (Glendinning et al., 2003; Dimonte et al., 1996; Dimonte and Remington, 1993; Aleshin et al., 1990). In our SPHC simulations, for ideal gases with $γh(l)=5/3$⁠, this characteristic scale is $[ka0]max=2.4$⁠. For $γh(l)<5/3$ its value overall increases (since the value of C₂ overall decreases) and leads to $[ka0]max∈[2.4,3.3]$ for a broad range of adiabatic indexes $γh(l)∈[1.2,1.6]$⁠.

The value $v0/v∞$ quantifies the fraction of energy imparted into the interface by the shock, and these results show that on average, ∼55% of the bulk velocity is available for interfacial mixing. From Table IV we know that at on average ∼40% of the shock velocity is imparted into the bulk motion, meaning that on average, only ∼25% of the shock velocity (and only ∼6% of the shock kinetic energy) is available for interfacial mixing. Because the area we are interested in is embedded in large, high-speed flows and there is only a small proportion of the total energy available for our interfacial mixing, it is quite difficult to capture these dynamics accurately. Because of this, it is beneficial to compare our results with theoretical analysis.

For small $a0/λ, a0/λ≤0.1$⁠, the initial growth rate v₀ linearly depends on a₀. Linear theory can predict this growth for these small amplitudes, as seen in Fig. 12 (Dell et al., 2015). We compare this linear theory with our simulation results. Restating our equation relating interface growth rate and initial perturbation amplitude, we have

The linear theory finds values of $v0/(v∞·a0λ)$ for small $a0/λ$⁠, so we choose the smallest initial amplitude of our data set, $a0/λ=0.1≪1$⁠. We rearrange Eq. (15) to get

where $[·]T$ is the value obtained from the linear theory. In order to investigate the linear approximation more closely, we ran simulations over a finer increment of $a0/λ$⁠, increasing by 0.05 instead of 0.1, as displayed in Fig. 12.

The approximations for linear theory assume the perturbation amplitude to be very small, $ka0≪1$ or $a0/λ<0.05$⁠, but we see that the approximation holds for larger amplitudes, up to $a0/λ<0.1$⁠, which is a result consistent with previous studies (Dell et al., 2017; Dell et al., 2015).

We compare the values of the quantity $(v0/v∞)/(a0/λ)$ and the quantity C₁ A exp(−0.1C₂) in the linear theory and in the simulations. Table VI displays the percentage error between the theory and the simulations, and Table VII presents the averaged values of $(v0/v∞)/(a0/λ)$⁠. The simulation results are in good agreement with the theoretical values, with the smallest error being along the diagonal where $γl=γh$⁠, which is the same observation we made with the error for $v∞$⁠. Figure 13 displays the average error for each Mach and Atwood number, and we can clearly see that the error increases as the Mach and Atwood numbers increase. This is the same observation we made with the errors for $v∞$⁠, because higher Mach and Atwood numbers contain more extreme densities, speeds, and shocks and are harder for simulations to model accurately. Notably, the average error for all cases that do not have an Atwood number of 0.95 is under 10%, which is very good. The Atwood number 0.95 is a very extreme case, with density contrasts of 40:1 which are a challenge to model accurately.

Despite the very small scale of interfacial mixing that we are measuring, our results are in good agreement with the theory, demonstrating the ability of SPHC to capture small-scale dynamics embedded in large-scale dynamics very accurately.

In this study, through use of numerical simulations we have systematically studied the effect of a previously unconsidered regime of adiabatic indexes of the fluids on the early time dynamics of RMI, specifically the extent of the interfacial mixing. The key properties of the dynamics we have analyzed are the interface velocity and the interface growth rate. Our regime included the Mach numbers 3 and 5, the Atwood numbers 0.6, 0.8, and 0.95, a range of adiabatic indexes $γl,γh=(1.2,1.3,…,1.6)$⁠, and a range of initial perturbations from 0% to 100% of the wavelength. In this regime, the simulations are repeatable and qualitatively similar, and we found good quantitative and qualitative agreement between the simulation results and the theory, demonstrating the accuracy at which SPHC simulations are able to capture small-scale dynamics embedded within large-scale dynamics in extreme and challenging situations.

In order to find the appropriate scale for the amount of energy deposited into the interface by the shock, we obtained the velocity of the background motion. In experiments, the background motion makes reliable diagnostics of RMI challenging because flow measurements must be taken from a quickly moving interface. Numerical simulations have the ability to scale the dynamics by the bulk velocity via the use of a Galilean transformation to a moving frame of reference. In order to scale our results, we obtained the bulk velocity $v∞$ in Sec. III B, and we found that the more challenging cases for the simulations to model occur when the adiabatic indexes of the two fluids are further apart (Table II), which to the best of the author's knowledge is the first time this observation has been made.

In order to quantify the amount of energy deposited into the interface by the shock, we found the initial growth rate of the interface v₀, and scaled it by the velocity of the background motion $v∞$⁠. The speed at which the interface spreads out indicates how much energy the interface received from the initial shock. We found that the initial growth rate is a non-monotone function of the initial perturbation amplitude, and that this relationship holds across all adiabatic indexes of the fluids (Fig. 6). For each pair of adiabatic indexes, we found the maximum values of the growth rate scaled by the background motion, $v0/v∞$⁠, and at which values of $a0/λ$ this occurred. We found that this maximum generally increases as γ_l and γ_h increase (Fig. 11), and the perturbation size at which this maximum occurs decreases as the adiabatic indexes increase (Fig. 10). We found that on average, only ∼25% of the shock velocity (only ∼6% of the shock kinetic energy) is transferred to the growth of the interface, indicating that only a fraction of the energy from the shock is deposited into the interface and is hence available for interfacial mixing. We compared our simulation results to the theory, and found the average error to be less than 7% for the bulk velocity (Table III) and less than 10% for the interface growth rate (Table VI) when the Atwood numbers were less than 0.95 (which is an extreme case). This indicates that the SPHC simulations can handle small-scale dynamics embedded in large-scale dynamics with high accuracy.

We also found a relationship between C₁, C₂, and γ_l, γ_h. This is useful in a variety of applications because the adiabatic index of the fluids is often hard to ascertain due to the highly complex molecules that form during the high-energy environment. Knowing the speed and amplitude of the interface, the constants C₁ and C₂ can be found, and using our results the adiabatic indexes can now be deduced, informing scientists of the molecules present.

Our results provide good benchmarks for further studies and experiments and open up further avenues of investigation for non-ideal adiabatic indexes. Our results also have implications for hydrodynamic instabilities and mixing in inertial confinement fusion (ICF). To achieve ICF ignition, the ability to avoid or control the RMI that forms during the implosion process is necessary (Lindl et al., 2004). One method is to fully suppress the development of RMI, which is based on traditional scenarios of RMI that suggest that the development of RMI may produce uncontrolled growth of small-scale imperfections and lead to disordered mixing that is similar to canonical turbulence (Dimonte et al., 2004). However, research conducted in studies like this through simulations and theoretical analysis suggests that the interfacial mixing may keep a significant degree of order, shown by the ability to accurately predict the evolution of the interface through theoretical analysis. These findings suggest that the dynamics can be controlled through the initial perturbation, so that interfacial mixing may be prevented without having to completely suppress RMI, which may be easier to implement (Dell et al., 2015; Anisimov et al., 2013; Abarzhi, 2010; Demskoĭ et al., 2006). This study also suggests there is reason to have confidence in the ability of the numerical simulations produced by SPHC in accurately capturing small-scale dynamics embedded in large-scale structures despite its difficulty.

The work is supported by the National Science Foundation in the USA (Award/Contract No. 1404449), by the University of Western Australia in Australia (Award/Contract No. 10101047), and by the Australian Government Research Training Program Scholarship in Australia. The authors deeply thank Dr. Robert F. Stellingwerf for valuable comments and discussions on the implementation and use of the Smoothed Particle Hydrodynamics code.

The methods, results and data in this work are available to interested researchers.