The effect of compressibility on the single-mode Rayleigh–Taylor instability is examined using two (2D) and three-dimensional (3D) direct numerical simulations. To isolate compressibility from background stratification effects, this work employs a constant density profile on each side of the interface. The numerical simulations are performed at various Reynolds numbers using the gas kinetic method for static Mach numbers up to *M* = 0.4. The most important finding is that compressibility acting in isolation enhances the instability and perturbations grows faster with increasing Mach number, unlike previous results with background isothermal state, which show suppression of the instability at higher static Mach numbers. In addition, compressibility is also shown to increase the bubble-spike asymmetry. While the instability grows faster for the 3D case, the findings are qualitatively similar in 2D and 3D. The dynamical reasons underlying the effect of compressibility are elucidated by examining the evolution of vorticity and turbulent kinetic energy transport equations.

## I. INTRODUCTION

The Rayleigh Taylor instability (RTI) develops at the interface between a heavier fluid lying on top of a lighter fluid.^{1,2} This phenomenon involves many complex flow features as the heavy fluid accelerates toward the lighter one. The light fluid ascends in the shape of rounded structures called bubbles, and the heavy fluid descends through elongated structures called spikes. These bubble and spikes interact non-linearly and eventually transition to turbulence.^{3,4} Such instabilities are of much interest in engineering applications (inertial and magnetic confinement fusion^{5–7}) and natural phenomena (Type Ia supernovae and x-ray burst^{8–11} and convective flows in the atmosphere and oceans^{12,13}).

While most studies to date focus on the incompressible or nearly incompressible limit,^{3} in many cases, RTI occurs in regimes where the fluids exhibit compressibility effects. There are a growing number of studies addressing the influence of compressibility on RTI. These include single-mode^{14–17} and multi-mode numerical simulation.^{18} Multi-mode experimental studies^{19,20} have also been performed to understand the mixing efficiency of RTI for linear density stratifications. Owing to their simplicity, we focus on single-mode RTI in this study. Unlike the incompressible limit, there are significantly more parameters influencing compressible RTI flow physics.^{21} In particular, the type and strength of background stratification can play significant roles and obscure the intrinsic compressibility effects associated with acoustic or entropic modes and their coupling back to vorticity.^{4,22} In addition, compressible RTI can occur for both two-fluid configurations, where the fluids have different molar masses (e.g., at early inertial confinement fusion (ICF) stages between the shell and the interior gas), as well as single-fluid configurations, where density variations are connected to temperature changes through the equation of state (e.g., at late ICF stages between the hot spot and colder surrounding hydrogen plasma). Two-fluid configurations allow for more complexity in the RTI development, as material properties (e.g., ratio of specific heats^{21} or viscosities^{23}) can be different. Most studies to date of the two-fluid case start from thermal equilibrium, which, coupled with hydrostatic equilibrium, leads to an initial exponential variation of the background density in each of the two fluid regions. Weiland *et al.*^{14} performed a comprehensive analysis of the background isothermal stratification strength effects on single mode RTI. The stratification strength was characterized by a static Mach number, $ M = g \lambda a 0$, defined as the ratio of the free fall velocity (over the distance *λ*) to the speed of sound. Here, *g* is the gravitational acceleration, *λ* is the wavelength of the initial perturbation, and *a*_{0} is the isothermal speed of sound. For the isothermal case, *M* defines the exponent of the background density variation.^{21} Weiland *et al.*^{14} showed that increasing *M* can result in full instability suppression. Using the vorticity transport equation, they were able to correlate the suppression of the instability to a decrease in incompressible baroclinic torque with increasing stratification. Nevertheless, out-of-thermal-equilibrium initial background stratifications, e.g., isochoric or isentropic background thermodynamic variations, do not necessarily lead to instability suppression and may have different effects at early and late times.^{22} In particular, isochoric initial stratification may also offer a better way of separating intrinsic compressibility effects from those due to stratification.

To distinguish the effects of different stratification types, Weiland *et al.*^{14} focused on secondary instabilities on the sides of the bubbles and spikes and considered the behavior of a self-propagating vortex pair due to the induced vortical velocity in different density fields. The motion of the vortex pair can predict the suppression or acceleration of the instability with *M* for different types of density backgrounds. These simulations were performed mainly at low Atwood number, where $ A = ( \rho h \u2212 \rho l ) / ( \rho h + \rho l )$, with *ρ _{l}* and

*ρ*the densities of the light and heavy fluids, respectively, when the instability is mostly symmetrical between the bubble and spike sides. Nevertheless, the instability becomes more asymmetric at higher

_{h}*M*values for the isothermal background case.

^{14}Several questions are still left unanswered regarding the intrinsic compressibility effects that can best be evaluated using isochoric initial stratification, e.g., with respect to the acceleration or suppression of the instability, asymmetry between the two sides, and differences between the two- and three-dimensional (3D) cases.

Another line of research with respect to the single mode RTI case has been the long time behavior and its dependence on Reynolds and Atwood numbers.^{15,24–27} While the potential flow theory and buoyancy-drag models indicate that the instability may reach a constant (terminal) velocity at late times, it is well known that the spike side does not follow this behavior at high enough Atwood numbers and even the bubble side can undergo a re-acceleration at high enough Reynolds numbers.^{24,28} Wei and Livescu^{24} showed, as a counterexample to the terminal velocity assumption, that two-dimensional (2D), low Atwood number incompressible single-mode RTI reaches mean quadratic growth at late time, if the perturbation Reynolds number, *Re _{p}* is larger than about 10 000. In this case, the vortex pairs generated on the sides of the bubbles and spikes self-propagate toward the edges of the layer and add bursts of acceleration leading to a fluctuating growth behavior with a quadratic mean, similar to the multi-mode case. Such growth mechanism is different than simple “bubble mergers” or “competition,” thought to describe RTI growth

^{29}and relies on the induced vortical velocity pointing in the vertical direction. At low Reynolds numbers, the vortices dissipate before reaching the edges of the layer; however, under these conditions, a terminal velocity is not maintained at long times.

Extensions to the compressible case were investigated by Reckinger *et al.*^{30} and Bian *et al.*^{15} for the two-fluid and single-fluid cases, respectively. At high Atwood numbers, the faster growth of the spike side tends to inhibit the vortices from reaching the edge of the layer on the bubble side by pulling them to the opposite side and, thus, a constant velocity growth on the bubble side appears more likely. This argument is consistent with the spike velocity reaching free fall descent as Atwood number tends to unity. The vortices have to travel longer distance before entering the bubble tip region. Hence, the region of largest vorticity production occurs near the spike tip, away from the bubble tip.^{15} However, as *Re _{p}* approaches 20 000, sustained quadratic growth on the bubble side was achieved at moderate Atwood numbers and significant re-acceleration effects were seen even at Atwood number of 0.8. Thus, the results do not preclude the possibility of quadratic growth on the bubble side at high enough

*Re*. That is, the viscous dissipation does not overcome vortex acceleration due to induced vortical velocity, at any given Atwood number. In addition, Bian

_{p}*et al.*

^{15}also showed that 3D cases grow faster than 2D cases and likely reach quadratic growth on the bubble side at lower Reynolds numbers, which is consistent with the faster self-propagation of vortex rings (for the 3D case) compared to vortex pairs (for the 2D case). Thus, sustained quadratic growth was seen at Atwood number of 0.8 for $ R e p = 8000$ in the 3D case. The single-fluid simulations of Bian

*et al.*

^{15}considered an initial isochoric background state; however, the effect of the Mach number on the results was not explored.

The focus of the present study is to understand the intrinsic compressibility effects on the single-fluid RTI development by considering the growth of single-mode RTI with an initial uniform density on each side of the interface. We investigate the effect of Mach number on the growth of RTI by performing direct numerical simulations (DNS) using gas kinetic method (GKM). Both two-dimensional and three-dimensional simulations are performed to establish the fundamental and 3D effects. The simulations are further analyzed using compressible vorticity, and turbulent kinetic energy transport equations and the effect of various flow processes pertaining to growth of vorticity and turbulent kinetic energy are identified.

The paper is organized as follows. Section II discusses the governing equations and the initial background state, investigated in this study. Section III describes the kinetic (as opposed to continuum) based computational approach used to solve the governing equations and the important flow processes describing the instability. Section IV describes the results of the simulations and goes into details about the effects of compressibility on vortex dynamics and kinetic energy interchanges during the evolution of the perturbation. Section V concludes the paper with a summary of the results.

## II. FLUID GOVERNING EQUATIONS

*i*by $ u i *$, temperature by $ T *$, and pressure by $ p *$. $ g i * = ( \u2212 g , 0 , 0 )$ is the gravitational acceleration, and $ e *$ is the total (kinetic + internal) energy. Furthermore, $ R *$ is the specific gas constant and $ c v *$ is the specific heat capacity at constant volume. The viscous stress tensor ( $ \tau i j *$) and the heat flux vector ( $ q i *$) are given by

*δ*is the Kronecker delta function.

_{ij}### A. Background state

^{22}As mentioned in Introduction, extensive studies have been conducted of RTI with initial background state being in thermal equilibrium.

^{14,18,21,30}In this work, to isolate compressibility effects, we consider single-fluid RTI with an initial isochoric background state, similar to Bian

*et al.*

^{15}The two-fluid version of this flow, which combines compressibility and variable density effects, has been considered in Ref. 22. The initial setup of the flow is shown in Fig. 1. The background density ( $ \rho 0 *$) is constant on either side of the interface and follows an error function profile across the interface,

^{15,24}

*Y*is the parameter controlling the steepness of the error function profile and is set to 34.

_{v}^{15}The background pressure $ p 0 *$ can be obtained from the hydrostatic condition as follows:

*p*is the pressure at the interface. Finally, the background temperature ( $ T 0 *$) is computed from the ideal gas equation of state (1d). The dimensional variables $ \rho * , \u2009 p *$, and $ T *$ are normalized by their respective interface ( $ x 1 * = 0$) values: $ \rho I = ( \rho l * + \rho h * ) / 2$,

_{I}*p*, and $ T I = p I / ( \rho I R * )$. Coordinates $ x i *$ are normalized with the wavelength of the initial perturbation of the interface,

_{I}*λ*. The velocity $ u i *$ is normalized by $ A g \lambda $, and time is non-dimensionalized by $ \lambda A g$ ( $ \tau = t * A g \lambda $). This time normalization is consistent with the linear analysis

^{21}and has been shown to collapse the evolution of various quantities in variable-density buoyancy driven turbulence.

^{31–33}The non-dimensional form of the governing equation (1)–(3) can then be written as

*γ*is the ratio of specific heats. The Reynolds, Froude, and Prandtl numbers are defined as

*A*is the Atwood number and

*M*is the interface Mach number based on the free fall velocity over the distance

*λ*and the isothermal speed of sound. The non-dimensional parameters

*A*and

*M*are defined as

^{21}wherein the linearized Navier–Stokes equations are examined for the case of compressible ideal fluid. The growth rate exponent is found to depend on the Mach number as defined here. The corresponding incompressible limit is achieved in the limit of speed of sound tending to infinity. For the Rayleigh–Taylor instability, the limits of infinitesimal

*A*and

*M*are not equivalent. Instability still exists at

*M*= 0, but it does not at

*A*= 0. The static Mach number is different from the acoustic Mach number (which is zero initially for RTI). The isothermal sound speed is used in the definition of

*M*instead of the isentropic sound speed (a factor of

*γ*difference), as the ratio of specific heats has a different role in the evolution of the instability. Additional discussion on the role of different non-dimensional parameters for the compressible RTI can be found in Ref. 23. A perturbation Reynolds number (

*Re*) can be defined based on the potential terminal velocity (

_{p}*U*),

_{p}*M*quantifies compressibility effects. This is in contrast to the thermal equilibrium case, wherein

*M*characterizes effects of both compressibility and the background stratification level. For the thermal equilibrium case, increasing Mach number results in increasing density stratification, since

*M*also appears in the exponent of the background density variation. Changing the Mach number results in a change in pressure at the interface, which corresponds to pure compressibility effect, and also to a variation in stratification.

^{14}The background density and pressure profiles for

*A*= 0.04 at different

*M*values are shown in Fig. 2. It is evident from the figure that the background density on either side of the interface is constant and independent of the Mach number. The background pressure profile is linear away from the interface, and the pressure gradient ( $ d p o / d x 1$) increases quadratically with

*M*. It is worth noting that, in contrast to the isochoric background state, both density and pressure exhibit exponential variation away from the interface in the isothermal case.

^{14}Moreover, unlike the current case, the density stratification level increases with increasing

*M*for the isothermal background case.

## III. COMPUTATIONAL APPROACH

In ICF applications and astrophysical flows, RTI occurs over a wide range of Mach, Knudsen and Reynolds numbers. The long-term objective of the present work is to examine kinetic effects on RTI that are beyond the scope of the Navier–Stokes equations. Toward this end, we perform direct numerical simulations (DNS) of single mode Rayleigh–Taylor instability using the gas kinetic method (GKM).^{34–37} The GKM approach has a fundamental advantage over traditional NS based methods, as the fluxes are calculated directly from the particle distribution function allowing for description of certain features those are not easily captured by traditional continuum approaches. In future works, we will employ the non-equilibrium version of GKM to examine RTI under more complex kinetic effects.

*t*is time. Here, $ E \u2192$ is the particle acceleration due to external forces. The function

*f*lies in a six-dimensional phase space consisting of the three spatial coordinates and three velocity coordinates. The term on the RHS [ $ ( \u2202 f \u2202 t ) collisions$] is modeled using the Bhatnagar–Gross–Krook model,

^{38}

^{,}

*g*is the Maxwellian distribution. The Boltzmann–BGK model is solved in the current method to obtain the numerical fluxes. The details of GKM can be found in the following references.

^{34–36}Here, we provide the main features.

*U*contains the conservation variables $ [ \rho , \rho u i , \rho e ]$ and

*F*is the flux of the variables

*U*. These macroscopic quantities are stored at the cell center; however, for calculation of fluxes we need them at cell interface. Thus, a reconstruction step is required to interpolate cell center values at the cell interface. A fifth order WENO scheme is used for the purpose of this reconstruction.

^{39}The fluxes at the cell interface are then evaluated from

*ξ*is an internal energy variable, and $ d \Xi = d c i d \xi $ is a volume element in phase space. These fluxes are calculated as moments of the particle distribution function

*f*. The integration here is done over velocity coordinates in the phase space and over all other internal degrees of freedom,

*ψ*being the vector of the moments corresponding to fluxes of each conserved quantity,

*f*at the cell interface. After solving for

*f*, we update the macroscopic values at the cell center based on the fluxes as follows:

*B*is assigned the value of $ 0.25 \lambda $, so that an initial linear growth of the instability can be expected.

^{24}Periodic boundary conditions are applied in the streamwise direction. The top and bottom walls are treated as isothermal no-slip boundaries. Sixth order compact filtering is used to perform the simulations.

^{15,40}The simulation parameters, domain sizes, and grid resolution are listed in Table I.

Cases . | Re
. _{p} | M . | A . | Grid points . | $ L x 1$ . | $ L x 2$ . | $ L x 3$ . |
---|---|---|---|---|---|---|---|

Two-dimensional simulation | |||||||

Re1M1 | 1000 | 0.0855 | 0.04 | 256 × 3072 × 4 | $ 8 \lambda $ | λ | 1 |

Re1M2 | 1000 | 0.2 | 0.04 | 256 × 3072 × 4 | $ 8 \lambda $ | λ | 1 |

Re1M3 | 1000 | 0.3 | 0.04 | 256 × 3072 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M1 | 5000 | 0.0855 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M2 | 5000 | 0.2 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M3 | 5000 | 0.3 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M4 | 5000 | 0.4 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M1 | 8000 | 0.0855 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M2 | 8000 | 0.2 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M3 | 8000 | 0.3 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M4 | 8000 | 0.4 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Three-dimensional simulation | |||||||

Re1M1-3D | 1000 | 0.0855 | 0.04 | 128 × 128 × 1024 | $ 8 \lambda $ | λ | λ |

Re1M3-3D | 1000 | 0.3 | 0.04 | 128 × 128 × 1024 | $ 8 \lambda $ | λ | λ |

Cases . | Re
. _{p} | M . | A . | Grid points . | $ L x 1$ . | $ L x 2$ . | $ L x 3$ . |
---|---|---|---|---|---|---|---|

Two-dimensional simulation | |||||||

Re1M1 | 1000 | 0.0855 | 0.04 | 256 × 3072 × 4 | $ 8 \lambda $ | λ | 1 |

Re1M2 | 1000 | 0.2 | 0.04 | 256 × 3072 × 4 | $ 8 \lambda $ | λ | 1 |

Re1M3 | 1000 | 0.3 | 0.04 | 256 × 3072 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M1 | 5000 | 0.0855 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M2 | 5000 | 0.2 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M3 | 5000 | 0.3 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re2M4 | 5000 | 0.4 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M1 | 8000 | 0.0855 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M2 | 8000 | 0.2 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M3 | 8000 | 0.3 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Re3M4 | 8000 | 0.4 | 0.04 | 512 × 4096 × 4 | $ 8 \lambda $ | λ | 1 |

Three-dimensional simulation | |||||||

Re1M1-3D | 1000 | 0.0855 | 0.04 | 128 × 128 × 1024 | $ 8 \lambda $ | λ | λ |

Re1M3-3D | 1000 | 0.3 | 0.04 | 128 × 128 × 1024 | $ 8 \lambda $ | λ | λ |

*et al.*

^{15}The spike front ( $ h s = h s * / \lambda $) is located at the position of maximum density gradient along the line $ x 2 = L x 2 / 2$ (note that for these simulations $ L x 2 = \lambda $). Similarly, the bubble front ( $ h b = h b * / \lambda $) is the corresponding location along $ x 2 = 0$. The bubble-spike speeds defined as in Bian

*et al.,*

^{15}

^{,}

*et al.*

^{15}In addition, Fig. 4 shows that the results are also mesh converged at the highest Reynolds and Mach numbers considered here (case Re3M4).

## IV. RAYLEIGH–TAYLOR INSTABILITY-DRIVEN COMPRESSIBLE MIXING

In order to address several open questions in the literature regarding the compressible RTI, the effects of Reynolds and Mach numbers on the evolution of the instability are considered. First, the variations of bubble-spike heights and velocities are shown with varying Mach number for the two-dimensional case. Three-dimensional effects are discussed next, followed by analyses of the vorticity and kinetic energy transport equations in order to develop a physics-based explanation of the findings. In particular, the vorticity transport equation budget is examined both during early (linear) and later stages.

### A. Effect of Mach number on the instability growth for the two-dimensional case

Figures 5(a) and 5(b) present the evolution of bubble and spike heights ( $ h b / s = h b / s * / \lambda $). Unlike the isothermal background case,^{14,30} it is seen that increasing Mach number results in faster instability growth. This is evident in the evolution of both bubble and spike heights, where at a given time, both the bubble and spike fronts have displaced more from the initial interface with increasing Mach number. The destabilization effect of compressibility for the uniform density background seen here is in direct contrast to the stabilizing action witnessed for an initially isothermally stratified RTI. With an isothermal stratification, growth of RTI is inhibited as the magnitude of the stratification parameter is increased, with complete growth suppression at strong enough stratification.

Turning our attention to evolution of bubble and spike velocities shown in Figs. 5(c) and 5(d), we observe that at early times, in the linear regime, the speeds are similar for all Mach numbers. Near the incompressible limit (*M* = 0.0855), the bubble speed then reaches values consistent with the potential flow theory, before a re-acceleration, as established by previous studies.^{24,28} However, as the Mach number increases, the bubble side grows faster, by-passing a quasi-constant velocity stage around the potential flow result altogether and transitions directly to the reacceleration stage. At high enough Reynolds number and for low to moderate Atwood numbers, the bubble side experiences late-time quadratic growth, contrary to the potential flow theory.^{15,24} The Reynolds number here is not large enough for the effect to manifest near the incompressible limit. However, the results suggest that quadratic growth might be achievable on the bubble side at lower Reynolds numbers, at large enough Mach numbers. For instance, at *M* = 0.4, both bubble and spike sides exhibit quasi-linear growth of the tip velocities, which indicate overall quadratic growth. Interestingly, unlike the high Reynolds number incompressible case,^{24} here, the quadratic growth is not accompanied by strong oscillations in the bubble velocity, suggesting in the high Mach number case a growth mechanism other than through discrete vortex pairs being transported to the bubble tip and adding bursts of acceleration. At larger Atwood numbers, for the incompressible case, there is an increase in asymmetry between the bubble and spike sides. The spike side does not reach a plateau around the potential flow theory result. For the compressible case, the destabilization is seen on both the bubble and spike sides.

Figure 6 shows density contours at three different times for different Mach numbers. At *τ* = 2, we see that bubble and spike heights are comparable, but with progression in time, the spike front is displaced more with increasing Mach number. Increasing Mach number also affects the shape of the mixing layer. It can be observed that the spike front is thinner at higher Mach numbers, with larger density gradients developing at the spike tip. In addition, low density spots form near the bubble front and high density spots form near the spike front.

Another feature that comes to light from the evolution plots in Fig. 5 is that at a given time, the spike height is larger than the bubble height. This is what is referred to as bubble-spike asymmetry. The asymmetry is exhibited in Fig. 8(a), where the difference of bubble and spike heights are plotted at different Mach numbers. The height difference remains around zero for $ \tau < 3$, with the asymmetry developing at later times and becoming stronger at higher M values. For the lowest Mach number of *M* = 0.0855, $ h s \u2212 h b$ remains near zero at all times, consistent with the low Atwood incompressible limit.

*x*

_{3}component is non-zero so that $ | \omega \u2192 | = | u 1 , 2 \u2212 u 2 , 1 |$. Figure 8(b) shows the evolution of the ratio $ \Gamma s / \Gamma b$. At low Mach number, this ratio remains close to unity, but it increases with increasing Mach numbers. Evidently, more vorticity is produced near the spike side than the bubble side, resulting in an increase the vortical asymmetry, similar to asymmetry seen in Fig. 8(a).

### B. Three-dimensional effects

This section examines the three-dimensional effects on the flow. The three-dimensional simulations are performed at Re = 1000 for two Mach numbers, M = 0.085 and M = 0.3. Figure 9 shows the evolution of bubble and spike characteristics at the two Mach numbers. The bubble-spike velocities reach the potential theory velocity for 3D near the incompressible limit and exhibit a temporary plateau around this velocity. The bubble-spike fronts re-accelerate from this velocity, owing to the vorticity propagation near the fronts.^{15} However, at higher *M* values, the destabilization effect of compressibility causes the bubble and spike velocities to quickly pass through the potential flow theory value, without reaching a temporary plateau.

Figure 10 compares the bubble-spike heights and velocities evolution in two and three dimensions. It is seen that the destabilization effect with Mach number is more prominent in three dimensions. Three-dimensional effects also result in earlier departure from the potential flow velocity. The density iso-surface for the case with *M* = 0.3 is shown in Fig. 11 at *τ* = 4, as the non-linear effects become strong. Intricate structures in three-dimensional case can be seen compared to two-dimensional simulations, consistent with the results in Ref. 15. This finding supports the inference that the flow would transition to chaotic growth mean behavior at lower Reynolds numbers in the three-dimensional case.

### C. Vorticity evolution

Vorticity has been a key parameter of study for RTI, as the misalignment between pressure and density gradients is responsible for the main instability growth mechanism, via the baroclinic production of vorticity.^{31} For the single mode case, the production of vorticity on the sides of the RTI bubbles and spikes through the Kelvin–Helmholtz instability and subsequent self-propagation of vortex pairs (2D) and rings (3D) toward the layer edges are directly connected to the layer growth re-acceleration^{15,24,28} and transition to mean quadratic growth at late times.^{15,24} Baroclinic production of vorticity appears in many compressible flows and has been intensely studied, e.g., in shock-turbulence interaction,^{41–43} Richtmyer–Meshkov instability,^{44–46} and flows involving chemical reaction.^{47,48} These flows also highlight an additional mechanism of vorticity production in compressible flows, involving vorticity-dilatation interaction. The role of vorticity on compressible RTI growth has been studied by Gauthier,^{18} who highlighted the importance of baroclinic production and the interaction with helicity. The effect of isothermal background stratification on the growth of vorticity is explained in Ref. 14, focusing on the baroclinic production due to the initial background stratification. Wieland *et al.*^{22} were the first to discuss the role of different types of background stratification for the two-fluid case. Bian *et al.*^{15} addressed the single-fluid case and focused on the transport of vorticity toward the bubble tip, as the growth mechanism related to re-acceleration and late time quadratic growth for the single-mode case. Here, we extend these studies by explaining the increased instability growth and asymmetry with Mach number for the single-fluid case with constant initial background density. As explained above, this setup also isolates the intrinsic compressibility effects from those due to stratification.

In Fig. 12, we compare vorticity contours at different Mach numbers. The contours look similar at *τ* = 2, when the evolution is still close to the linear regime, for the three different Mach numbers considered here. There is not much difference in vortical content at this time. At *τ* = 3, the secondary vortices begin to develop, indicating the onset of Kelvin–Helmholtz instability at the deforming interface.^{28} The effect of Mach number manifests at this time. The vortex cores have descended more for higher Mach numbers, stretching the deforming interface. The pair of vortices generated near the centerline split into two pairs those are quickly transported toward the spike and bubble tips. These individual pairs can be seen close to the two fronts at *τ* = 4. The evolution plot in Fig. 13 shows the overall increase in vorticity ( $ \omega \u2192 = \omega \u2192 * A g / \lambda $) as the instability grows, with additional vorticity being generated as secondary Kelvin–Helmholtz instability develops at the deforming interface. With increasing Mach number, there is an increase in total vortical content consistent with more rapid destabilization of the flow.

^{14}

*x*

_{2}in 2D and

*x*

_{2}and

*x*

_{3}in 3D), shown below for an arbitrary quantity $T$,

Baroclinic torque is seen to be the major source of vorticity production. We see different trends of the torque profiles for different Mach numbers. For *M* = 0.0855, we see that the peak values occur near the interface, slightly tilting toward the spike side. However, for the higher Mach number cases, the peak starts shifting more toward the spike side, resulting in more vorticity production on this side. Dilatation effects are negligible at earlier times. As the instability evolves, sharper gradients develop, resulting in prominent dilatational effects near the spike and bubble front. The contribution from viscous dissipation also grows with time. Since baroclinic torque is the primary contributor to vortex growth, we will examine this term in greater detail.

*β*

_{1}represents the baroclinic torque generated due to non-linear interactions between the fluctuating density and pressure gradient.

*β*

_{2}and

*β*

_{3}are the contributions to baroclinic torque resulting from the mean field interacting with the fluctuating field. During the initial stages of RTI evolution, when non-linear effects are not important, the contribution of

*β*

_{1}toward the baroclinic torque is negligible. For an isochoric background state, the

*β*

_{3}component is also negligible during this stage away from the interface. Consequently, the initial vorticity dynamics near the spike and bubble tips is primarily affected by the mean pressure gradient term.

#### 1. Linear analysis for the early time behavior

*β*

_{1}is quadratic in perturbation;

*β*

_{2}depends on the mean pressure gradient; and

*β*

_{3}depends on the background density gradient.

*ρ*and

_{h}*ρ*. To verify this using the numerical data, we focus on

_{l}*β*

_{2}, which gives the non-zero contribution to $ T 3$ during the early times evolution. We define averages around the bubble and spike tip regions as follows (as depicted in Fig. 7):

Figure 15(b) shows the evolutions of $ \beta 2 s pike$ and $ \beta 2 b ubble$. It is seen that for $ \tau < 0.8$, the results are very close for the bubble and spike sides. At slightly later times, differences start to appear, corresponding to the asymmetric vortical growth for the bubble and spike sides. The onset of asymmetry occurs earlier with increasing Mach number, which is consistent with the *M*^{2} weighing of the density in formula (29). Hence, at early times, the baroclinic torque due to the background pressure gradient generates vorticity in an asymmetrical manner between the bubble and spike sides, which is associated with the inherent background pressure variation on the heavy and light fluid sides of the interface.

#### 2. Non-linear analysis

We next investigate late time behavior using the non-linear decomposition presented in Sec. IV C. Figure 16 shows the profiles of the individual portions of the baroclinic torque. *β*_{2} is found to be the dominant term of the three contributions. *β*_{2} has a bubble-spike asymmetry that is also indicative of the differential growth at early times. The non-linear term, *β*_{1}, shows a bubble-spike asymmetry but has an overall positive effect on the growth of vorticity. *β*_{3} also shows a bubble-spike asymmetry and aids in vorticity production near the spike side but destroys vorticity near the bubble side. Hence, it contributes to the overall enhancement of the asymmetry.

Figure 17 shows the growth of the three terms with time. At early times, *β*_{2} is the dominant production mechanism for vorticity, with the other two contributions becoming relevant at later times. Upon increasing Mach number, the enhancement in total baroclinic torque is seen in all three terms. This is different than the isothermal stratification case, where the baroclinic torque due to perturbation gradients was found to be the dominant contributor^{14} at weak background stratifications. On the other hand, terms based on background pressure gradient are found to become prominent at higher values of stratification and to inhibit the growth of vorticity.

### D. Kinetic energy evolution

^{49}

*K*, and turbulent kinetic,

_{m}*k*, energies,

The kinetic energy level clearly increases with Mach number, as the instability grows faster. For the lowest Mach number, *k* asymptotes to a limiting value, at late times. This is consistent with the results for the growth of the layer width for the relatively low *Re* values considered here and has been discussed in Sec. IV A. Increasing Mach number results in increased growth rate in the linear regime and kinetic energy growth at later times, indicative of a faster than linear growth of the layer width. Figure 18(b) shows the profiles of kinetic energy. There is a stark difference in the profiles with increasing Mach number. The increase in global kinetic energy can be seen in the profiles as well. However, the increase in kinetic energy is more prominent near the spike front than the bubble front, once again highlighting the bubble-spike asymmetry.

^{31,49}

^{32}The normalized governing equation for fluctuating kinetic energy can be expressed in the following form:

^{31}

^{,}

^{31}$ R i j = \u27e8 \rho u i \u2033 u j \u2033 \u27e9$ is the Favre averaged Reynolds stresses and $ d = \u2202 u i \u2032 \u2202 x i$ is the fluctuating dilatation. The mass flux–pressure gradient term ( $ K 3$) is the primary source of kinetic energy production for RTI flows. $ K 3$ represents the indirect effect of the gravitational potential and transfers energy from the mean velocity field to the fluctuating field. In contrast to shear flows, the production due to mean velocity gradient ( $ K 1$) is negligible for RTI. The pressure-dilatation term ( $ K 2$) facilitates a reversible exchange between the fluctuating kinetic and internal fields. The energy transfer from the kinetic to the internal field via dissipation ( $ K 4$) is irreversible. Finally, the viscous action ( $ K 5$) can also lead to a two-way exchange between the mean and fluctuating field.

Next we examine the various terms in the budget of turbulent kinetic energy transport equation (36) in Fig. 19. Individual terms in the budget are plotted for three different Mach numbers at *τ* = 4. It is seen that mass flux–pressure gradient acts as the dominant source in transferring energy from the mean flow. For higher Mach numbers, this term develops a bubble-spike asymmetry; hence, there is more transfer of energy at the spike front. Figure 20 shows the turbulent kinetic energy budget for the three-dimensional cases. Again, it can be seen that the mass flux–pressure gradient ( $ K 3$) term acts as the major contributor for evolution of kinetic energy. The profiles of $ K 3$ also exhibit bubble-spike asymmetry.

Figures 21 and 22 show the profiles of mass flux in vertical direction (*a*_{1}) and those of the mean pressure gradient. Both mass flux and pressure gradient are negative at all points, hence implying a transfer of energy from mean flow to the fluctuations. With an increase in Mach number, stronger negative pressure gradient develops near the spike front, indicating stronger extraction of kinetic energy from the mean flow, resulting in asymmetrical growth of turbulent kinetic energy. Figures 23 and 24 show the correspondent profiles for the three-dimensional case. Similar to the 2D case, the asymmetry in $ K 3$ term arises primarily from the asymmetry in mean pressure gradient.

## V. CONCLUSIONS

The effect of compressibility on the Rayleigh Taylor Instability (RTI) is examined using direct numerical simulations (DNS) performed with the gas-kinetic method (GKM). Previous studies of compressibility effects on RTI employ isothermal base profiles, and the resulting flow combines the effects of compressibility and stratification.^{14,30} Subject to these joint effects, it was found that the RTI was stabilized at higher static Mach numbers. These studies do not clearly establish the individual actions of compressibility and stratification on the suppression of the instability. Toward that end, this work isolates the effect of compressibility by employing an initial isochoric (as opposed to isothermal) background profile. Two- and three-dimensional simulations are performed over a range of Reynolds numbers (1000, 8000) at several static Mach numbers (0.085, 0.3, 0.4). Contrary to the conclusions of the previous isothermal background studies, it is found that compressibility, in isolation, has a destabilizing effect. Furthermore, compressibility enhances the asymmetry between bubble and spike growths. The destabilization and asymmetry effects are more pronounced in three-dimensional simulations. Thus, faster destabilization is seen for the three-dimensional cases. Turbulent kinetic energy is also seen to grow faster, consistent with the destabilizing nature of compressibility.

To develop further insight into the two-pronged effect of compressibility, the vorticity dynamics is examined using DNS data at various Mach and Reynolds numbers. Baroclinic torque is found to be the primary source of vorticity production, with dilatational effect becoming important at later times. Three distinct contributions of baroclinic torque are identified—non-linear effects (*β*_{1}), mean pressure-gradient effect (*β*_{2}), and mean density-gradient effect (*β*_{3}). It is found that the contribution of *β*_{2} component is most significant toward vorticity growth, especially at early times. At early times, the *β*_{2} contribution can be directly expressed using the background pressure variation and exhibits different magnitudes on the spike and bubble sides, proportional to *ρ _{h}* and

*ρ*, respectively. This difference inherently initiates the bubble-spike asymmetry seen at later stages. Furthermore, since

_{l}*β*

_{2}depends on the Mach number, the asymmetry tendency is stronger at higher compressibility levels. Once the bubble-spike asymmetry is initiated, the rest of the vorticity budget terms responds and further amplifies this asymmetry.

The turbulent kinetic energy budget is also examined. The mass flux–pressure gradient term ( $ K 3$) is the primary process through which energy is transferred from the mean flow to the fluctuating field. Both *β*_{2} and $ K 3$ depend on the mean pressure gradient, which becomes more negative near the spike and more positive near the bubble side with increasing Mach number. Thus, the effects of compressibility observed in this study can be attributed to the characteristics of the evolution of the mean pressure gradient.

## ACKNOWLEDGMENTS

This work was performed under the auspices of DOE. S.M., B.S., and S.S.G. have been supported through Triad Award No. 628660. LANL, an affirmative action/equal opportunity employer, is managed by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy under Contract No. 89233218CNA000001.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Swapnil Majumder:** Data curation (equal); Formal analysis (equal); Writing – original draft (equal). **Bajrang Sharma:** Formal analysis (equal). **Daniel Livescu:** Conceptualization (equal); Formal analysis (equal); Writing – review & editing (equal). **Sharath S. Girimaji:** Conceptualization (equal); Formal analysis (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings are available from the corresponding author upon reasonable request.

## REFERENCES

*Presented at the 29th CFD, Lecture Series*1998-03, von Karman Insititute for Fluids Dynamics, February 23–27, 1998; available at https://api.semanticscholar.org/CorpusID:116261881