Compressible Rayleigh–Taylor instability subject to isochoric initial background state

Majumder, Swapnil; Sharma, Bajrang; Livescu, Daniel; Girimaji, Sharath S.

doi:10.1063/5.0163886

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,³ 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.¹⁸ 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.²¹ 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²¹ or viscosities²³) 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.¹⁴ 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 = \frac{\sqrt{g λ}}{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₀ is the isothermal speed of sound. For the isothermal case, M defines the exponent of the background density variation.²¹ Weiland et al.¹⁴ 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.²² 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.¹⁴ 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 = (ρ_{h} - ρ_{l}) / (ρ_{h} + ρ_{l})$ ⁠, with ρ_l and ρ_h 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 M values for the isothermal background case.¹⁴ 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²⁴ 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²⁹ 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.³⁰ and Bian et al.¹⁵ 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.¹⁵ 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_p. That is, the viscous dissipation does not overcome vortex acceleration due to induced vortical velocity, at any given Atwood number. In addition, Bian et al.¹⁵ 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.¹⁵ 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

We examine RTI that develops within the region separating high and low density zones in a single fluid. To isolate the intrinsic compressibility effects from those due to background stratification, the background density is uniform in each of the two zones, with the background conditions in hydrostatic equilibrium. In the continuum limit, the equations governing the flow are the compressible Navier–Stokes equations, in an ideal compressible fluid,

(1)

\begin{matrix} \frac{\partial ρ^{*}}{\partial t^{*}} + \frac{\partial ρ u_{j}^{*}}{\partial x_{j}^{*}} = 0, \end{matrix}

(1a)

\begin{matrix} \frac{\partial ρ u_{i}^{*}}{\partial t^{*}} + \frac{\partial ρ u_{i}^{*} u_{j}^{*}}{\partial x_{j}^{*}} = - \frac{\partial p^{*}}{\partial x_{i}^{*}} + ρ^{*} g_{i}^{*} + \frac{\partial τ_{i j}^{*}}{\partial x_{j}^{*}}, \end{matrix}

(1b)

\begin{matrix} \frac{\partial ρ e^{*}}{\partial t^{*}} + \frac{\partial ρ^{*} e^{*} u_{j}^{*}}{\partial x_{j}^{*}} = - \frac{\partial p^{*} u_{i}^{*}}{\partial x_{i}^{*}} + ρ u_{i}^{*} g_{i}^{*} + \frac{\partial τ_{i j}^{*} u_{i}^{*}}{\partial x_{j}^{*}} - \frac{\partial q_{j}^{*}}{\partial x_{j}^{*}}, \end{matrix}

(1c)

\begin{matrix} p^{*} = ρ^{*} R^{*} T^{*}, \end{matrix}

(1d)

\begin{matrix} e^{*} = \frac{1}{2} u_{i}^{*} u_{i}^{*} + c_{v}^{*} T^{*} . \end{matrix}

(1e)

Here, the superscript

^{*}

is used to denote the dimensional variables. The density of the fluid is denoted by

ρ^{*}

⁠, velocity in direction i by

u_{i}^{*}

⁠, temperature by

T^{*}

⁠, and pressure by

p^{*}

⁠.

g_{i}^{*} = (- 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 (⁠

τ_{i j}^{*}

⁠) and the heat flux vector (⁠

q_{i}^{*}

⁠) are given by

τ_{i j}^{*} = 2 μ^{*} S_{i j}^{*} - \frac{2}{3} μ^{*} \frac{\partial u_{k}^{*}}{\partial x_{k}^{*}} δ_{i j}, S_{i j}^{*} = \frac{1}{2} (\frac{\partial u_{i}^{*}}{\partial x_{j}^{*}} + \frac{\partial u_{j}^{*}}{\partial x_{i}^{*}}),

(2)

q_{j}^{*} = - κ^{*} \frac{\partial T^{*}}{\partial x_{j}^{*}} .

(3)

Here,

μ^{*}

and

κ^{*}

are the dynamic viscosity and the coefficient of thermal conductivity of the fluid, respectively, and δ_ij is the Kronecker delta function.

A. Background state

Rayleigh–Taylor instability (RTI) evolution is strongly influenced by the background state of the thermodynamic variables.²² 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.¹⁵ 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 (⁠

ρ_{0}^{*}

⁠) is constant on either side of the interface and follows an error function profile across the interface,^15,24

ρ_{0}^{*} = 0.5 [1 + erf (Y_{v} \frac{x_{1}^{*}}{L_{x_{1}^{*}}})] (ρ_{h}^{*} - ρ_{l}^{*}) + ρ_{l}^{*},

(4)

where

ρ_{h}^{*}

and

ρ_{l}^{*}

are the dimensional densities on the top and bottom sides of the interface, respectively. Y_v is the parameter controlling the steepness of the error function profile and is set to 34.¹⁵ The background pressure

p_{0}^{*}

can be obtained from the hydrostatic condition as follows:

p_{0}^{*} = p_{I} - ρ^{*} g^{*} x_{1}^{*},

(5)

where p_I 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

ρ^{*}, p^{*}

⁠, and

T^{*}

are normalized by their respective interface (⁠

x_{1}^{*} = 0

⁠) values:

ρ_{I} = (ρ_{l}^{*} + ρ_{h}^{*}) / 2

⁠, p_I, and

T_{I} = p_{I} / (ρ_{I} R^{*})

⁠. Coordinates

x_{i}^{*}

are normalized with the wavelength of the initial perturbation of the interface, λ. The velocity

u_{i}^{*}

is normalized by

\sqrt{A g λ}

⁠, and time is non-dimensionalized by

\sqrt{\frac{λ}{A g}}

(⁠

τ = t^{*} \sqrt{\frac{A g}{λ}}

⁠). This time normalization is consistent with the linear analysis²¹ 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 can then be written as

(6)

\begin{matrix} \frac{\partial ρ}{\partial t} + \frac{\partial ρ u_{j}}{\partial x_{j}} = 0, \end{matrix}

(6a)

\begin{matrix} \frac{\partial ρ u_{i}}{\partial t} + \frac{\partial ρ u_{i} u_{j}}{\partial x_{j}} = - \frac{1}{A M^{2}} \frac{\partial p}{\partial x_{i}} + \frac{1}{F r^{2}} ρ g_{i} + \frac{1}{R e} \frac{\partial τ_{i j}}{\partial x_{j}}, \end{matrix}

(6b)

\begin{matrix} \frac{\partial ρ e}{\partial t} + \frac{\partial ρ e u_{j}}{\partial x_{j}} = - \frac{1}{A M^{2}} \frac{\partial p u_{i}}{\partial x_{i}} + \frac{1}{F r^{2}} ρ u_{i} g_{i} + \frac{1}{R e} \frac{\partial τ_{i j} u_{i}}{\partial x_{j}} \\ - \frac{γ}{(γ - 1) ARePr M^{2}} \frac{\partial}{\partial x_{j}} (κ \frac{\partial T}{\partial x_{j}}), \end{matrix}

(6c)

\begin{matrix} p = ρ T, \end{matrix}

(6d)

\begin{matrix} e = \frac{1}{2} u_{i} u_{i} + \frac{1}{(γ - 1) A M^{2}} c_{v} T . \end{matrix}

(6e)

FIG. 1.

View large Download slide

Schematic of the initial RTI setup. $ρ_{h}^{*}$ and $ρ_{l}^{*}$ are the densities on the top and bottom sides of the interface, respectively. $L_{x_{1}}, L_{x_{2}}$ ⁠, and $L_{x_{3}}$ are the domain sizes in directions x₁, x₂, and x₃, respectively.

The viscous stress tensor is defined by

τ_{i j} = 2 μ S_{i j} - \frac{2}{3} μ \frac{\partial u_{k}}{\partial x_{k}} δ_{i j}, S_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}) .

(7)

Here,

κ = κ^{*} / κ_{r}

is the heat conduction coefficient,

μ = μ^{*} / μ_{r}

is the dynamic viscosity,

c_{v} = c_{v}^{*} (γ - 1) / R^{*}

is the specific heat at constant volume, and γ is the ratio of specific heats. The Reynolds, Froude, and Prandtl numbers are defined as

R e = \frac{ρ_{I} λ \sqrt{A g λ}}{μ_{r}}; F r^{2} = \frac{A g λ}{g λ} = A; P r = \frac{μ_{r} γ R^{*}}{k_{r} (γ - 1)} .

(8)

The non-dimensional form of the background profiles is then given by

\begin{matrix} ρ_{0} = A [1 + erf (Y_{v} x_{1})] + (1 - A) . \end{matrix}

(9a)

\begin{matrix} p = 1 - M^{2} ρ_{0} x_{1} . \end{matrix}

(9b)

Here, 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

A = \frac{ρ_{h} - ρ_{l}}{ρ_{h} + ρ_{l}}; M = \sqrt{\frac{ρ_{I} g λ}{p_{I}}} .

(10)

The definition of Mach number stems from the analysis of Livescu,²¹ 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_p) can be defined based on the potential terminal velocity (U_p),

R e_{p} = \sqrt{\frac{A}{1 + A}} R e; U_{p} = \sqrt{\frac{A g λ}{1 + A}} .

(11)

For an isochoric background state, the static Mach number 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.¹⁴ 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.¹⁴ Moreover, unlike the current case, the density stratification level increases with increasing M for the isothermal background case.

FIG. 2.

Background (a) density and (b) pressure profiles for A = 0.04 at different Mach numbers. (a) ρ0 and (b) p0.

View large Download slide

Background (a) density and (b) pressure profiles for A = 0.04 at different Mach numbers. (a) ρ₀ and (b) p₀.

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.

The Boltzmann equation describes the transport of a single-particle distribution function

f (\vec{x}, \vec{c}, t)

⁠,

\frac{\partial f}{\partial t} + \vec{c} . \nabla f + \vec{E} . \nabla_{c} f = {(\frac{\partial f}{\partial t})}_{collisions},

(12)

where

\vec{x}

is position,

\vec{c}

is velocity, and t is time. Here,

\vec{E}

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 [

{(\frac{\partial f}{\partial t})}_{collisions}

] is modeled using the Bhatnagar–Gross–Krook model,³⁸^,

{(\frac{\partial f}{\partial t})}_{collisions} = \frac{g - f}{τ},

(13)

where 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.

The governing equation in GKM can be cast as follows:

\frac{\partial}{\partial t} \int_{Ω} U d x + \oint_{A} \vec{F} . d \vec{A} = 0,

(14)

where the vector U contains the conservation variables

[ρ, ρ u_{i}, ρ 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.³⁹ The fluxes at the cell interface are then evaluated from

F_{i} = {[F_{ρ}, F_{ρ u_{i}}, F_{E}]}^{T} = \int_{- \infty}^{\infty} c_{i} ψ f (c_{i}, t, ξ) d Ξ,

(15)

where the fluxes of the three conserved quantities—mass, momentum, and energy—are given by

F_{ρ}, F_{ρ u_{i}}, F_{E}

⁠, respectively.

F_{i} = [F, G, H]

are the fluxes in the three coordinated directions, respectively. ξ is an internal energy variable, and

d Ξ = d c_{i} d ξ

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,

ψ = (1, c_{i}, \frac{1}{2} (c_{1}^{2} + c_{2}^{2} + c_{3}^{2} + ξ^{2})) .

(16)

The B–BGK equation is solved using the method of characteristics to obtain f at the cell interface. After solving for f, we update the macroscopic values at the cell center based on the fluxes as follows:

\begin{matrix} U_{i}^{n + 1} = U_{i}^{n} - \frac{1}{Δ x} \int_{t}^{t + \nabla t} (F_{i + 1 / 2} (t) - F_{i - 1 / 2} (t)) d t \\ - \frac{1}{Δ y} \int_{t}^{t + \nabla t} (G_{i + 1 / 2} (t) - G_{i - 1 / 2} (t)) d t \\ - \frac{1}{Δ z} \int_{t}^{t + \nabla t} (H_{i + 1 / 2} (t) - H_{i - 1 / 2} (t)) d t . \end{matrix}

(17)

Using the gas kinetic method, both two-dimensional and three-dimensional RTI cases are simulated. A density perturbation is superposed on the background field. For the 2D case, this has the form,

ρ^{'} = B cos (\frac{2 π x_{2}}{λ}),

(18)

while for the 3D case, it is

ρ^{'} = B cos (\frac{2 π x_{2}}{λ} + \frac{2 π x_{3}}{λ}) .

(19)

The perturbation amplitude B is assigned the value of

0.25 λ

⁠, so that an initial linear growth of the instability can be expected.²⁴ 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.

TABLE I.

List of simulations and parameters.

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 λ$	λ	1
Re1M2	1000	0.2	0.04	256 × 3072 × 4	$8 λ$	λ	1
Re1M3	1000	0.3	0.04	256 × 3072 × 4	$8 λ$	λ	1
Re2M1	5000	0.0855	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re2M2	5000	0.2	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re2M3	5000	0.3	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re2M4	5000	0.4	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M1	8000	0.0855	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M2	8000	0.2	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M3	8000	0.3	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M4	8000	0.4	0.04	512 × 4096 × 4	$8 λ$	λ	1
Three-dimensional simulation
Re1M1-3D	1000	0.0855	0.04	128 × 128 × 1024	$8 λ$	λ	λ
Re1M3-3D	1000	0.3	0.04	128 × 128 × 1024	$8 λ$	λ	λ

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 λ$	λ	1
Re1M2	1000	0.2	0.04	256 × 3072 × 4	$8 λ$	λ	1
Re1M3	1000	0.3	0.04	256 × 3072 × 4	$8 λ$	λ	1
Re2M1	5000	0.0855	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re2M2	5000	0.2	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re2M3	5000	0.3	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re2M4	5000	0.4	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M1	8000	0.0855	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M2	8000	0.2	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M3	8000	0.3	0.04	512 × 4096 × 4	$8 λ$	λ	1
Re3M4	8000	0.4	0.04	512 × 4096 × 4	$8 λ$	λ	1
Three-dimensional simulation
Re1M1-3D	1000	0.0855	0.04	128 × 128 × 1024	$8 λ$	λ	λ
Re1M3-3D	1000	0.3	0.04	128 × 128 × 1024	$8 λ$	λ	λ

A rigorous comparison study is performed to ensure the accuracy of the GKM code for RTI investigation. The present simulations are compared against the results of Bian et al.¹⁵ The spike front (⁠

h_{s} = h_{s}^{*} / λ

⁠) 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}} = λ

⁠). Similarly, the bubble front (⁠

h_{b} = h_{b}^{*} / λ

⁠) is the corresponding location along

x_{2} = 0

⁠. The bubble-spike speeds defined as in Bian et al.,¹⁵^,

F r_{b / s} = u_{b / s}^{*} / \sqrt{A / (1 + A) g λ} .

(20)

These values are obtained from the vertical velocity of the front locations. Figure 3 shows the comparisons for the evolutions of bubble/spike heights and speeds against those in Ref. 15 for case Re1M1. The results highlight convergence with increasing the mesh resolution and are very close to those reported in Bian et al.¹⁵ In addition, Fig. 4 shows that the results are also mesh converged at the highest Reynolds and Mach numbers considered here (case Re3M4).

FIG. 3.

Comparison of (a) bubble height, (b) spike height, (c) bubble speed, and (d) spike speed for Rep = 1000 with Bian et al.15

View large Download slide

Comparison of (a) bubble height, (b) spike height, (c) bubble speed, and (d) spike speed for Re_p = 1000 with Bian et al.¹⁵

FIG. 4.

View large Download slide

Grid convergence of: (a) bubble height, (b) spike height, (c) bubble speed, and (d) spike speed, for case with Re_p = 8000, M = 0.4.

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}^{*} / λ$ ⁠). 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.

FIG. 5.

View large Download slide

Effect of Mach number on evolution of (a) bubble height, (b) spike height, (c) bubble velocity, and (d) spike velocity, at Re_p = 8000, two-dimensional case.

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,²⁴ 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.

FIG. 6.

View large Download slide

Density contours are shown for cases with Re_p = 8000 and M = 0.0855 in (a)–(c), M = 0.2 in (d)–(f), and M = 0.3 in (g)–(i). (a) τ = 2, (b) τ = 3, (c) τ = 4, (d) τ = 2, (e) τ = 3, (f) τ = 4, (g) τ = 2, (h) τ = 3, and (i) τ = 4.

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 $τ < 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} - h_{b}$ remains near zero at all times, consistent with the low Atwood incompressible limit.

We next look at evolution of vorticity field (⁠

\vec{ω} = \nabla \times \vec{u}

⁠) near the bubble and spike sides, in the anticipation of the in-depth vorticity analysis in Sec. IV C. To quantify the vorticity accumulating near the spike and bubble fronts, we define circulation in a region around the fronts as

Γ_{s} = \int_{- L_{x_{1} / 2}}^{L_{x_{1} / 2}} \int_{0}^{λ / 4} | \vec{ω |} d x_{2} d x_{1},

(21)

Γ_{b} = \int_{- L_{x_{1} / 2}}^{L_{x_{1} / 2}} \int_{λ / 4}^{λ / 2} | \vec{ω |} d x_{2} d x_{1} .

(22)

The regions of integration around the bubble and spike fronts are shown in Fig. 7. The region from

0 < x_{2} < λ / 4

denotes the bubble side and from

λ / 4 < x_{2} < λ / 2

denotes the spike side of the instability. For the 2D case, only the x₃ component is non-zero so that

| \vec{ω} | = | u_{1, 2} - u_{2, 1} |

⁠. Figure 8(b) shows the evolution of the ratio

Γ_{s} / Γ_{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).

FIG. 7.

View large Download slide

Domain of integration for taking averages around bubble and spike fronts.

FIG. 8.

View large Download slide

Effect of Mach number on asymmetrical growth between bubble and spike based on (a) spike bubble/height difference and (b) vortical content ratio. All cases are two-dimensional, with Re_p = 8000.

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.¹⁵ 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.

FIG. 9.

View large Download slide

Effect of Mach number on the evolution of (a) bubble height, (b) spike height, (c) bubble speed, and (d) spike speed for 3D RTI.

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.

FIG. 10.

View large Download slide

Comparison of the evolution of (a) bubble height, (b) spike height, (c) bubble speed, and (d) spike speed between two- and three-dimensional cases. The Reynolds number for both the cases is 1000.

FIG. 11.

View large Download slide

Density iso-surfaces for 3D case: Re_p = 1000 and M = 0.3 at τ = 4.

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.³¹ 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,¹⁸ 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.²² were the first to discuss the role of different types of background stratification for the two-fluid case. Bian et al.¹⁵ 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.²⁸ 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 (⁠ $\vec{ω} = {\vec{ω}}^{*} \sqrt{A g / λ}$ ⁠) 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.

FIG. 12.

View large Download slide

Vorticity contours are shown for cases with Re_p = 8000 and M = 0.0855 in (a)–(c), M = 0.2 in (d)–(f), and M = 0.3 in (g)–(i). (a) τ = 2, (b) τ = 3, (c) τ = 4, (d) τ = 2, (e) τ = 3, (f) τ = 4, (g) τ = 2, (h) τ = 3, and (i) τ = 4.

FIG. 13.

Effect of Mach number on growth of total vorticity content, ω ¯ = 1 L x 1 L x 2 ∫ − L x 1 / 2 L x 1 / 2 ∫ − L x 2 / 2 L x 2 / 2 | ω → | d x 1 d x 2.

View large Download slide

Effect of Mach number on growth of total vorticity content, $\bar{ω} = \frac{1}{L_{x_{1}} L_{x_{2}}} \int_{- L_{x_{1} / 2}}^{L_{x_{1} / 2}} \int_{- L_{x_{2} / 2}}^{L_{x_{2} / 2}} | \vec{ω} | d x_{1} d x_{2}$ ⁠.

The vorticity (⁠

\vec{ω} = \nabla \times \vec{u}

⁠) transport equation in a compressible flow can be derived by taking the curl of the momentum equation . The governing equation for vorticity magnitude (⁠

Ω = | \vec{ω} |

⁠) can be expressed in the following form:¹⁴

\begin{matrix} \frac{D Ω}{D t} = \underset{T_{1}}{\underset{︸}{\hat{ω_{i}} ω_{j} S_{i j}}} - \underset{T_{2}}{\underset{︸}{Ω S_{k k}}} + \underset{T_{3}}{\underset{︸}{\hat{ω_{i}} ϵ_{ijk} \frac{1}{A M^{2} ρ^{2}} \frac{\partial ρ}{\partial x_{j}} \frac{\partial p}{\partial x_{k}}}} \\ + \underset{T_{4}}{\underset{︸}{\hat{ω_{i}} / R e \frac{\partial^{2} ω_{i}}{\partial x_{j} \partial x_{j}}}} - \underset{T_{5}}{\underset{︸}{\frac{2 \hat{ω_{i}}}{R e} ϵ_{ijk} \frac{\partial}{\partial x_{j}} [S_{k l} \frac{\partial l n (1 / ρ)}{\partial x_{l}}]}} . \end{matrix}

(23)

Here,

\hat{ω_{i}} = ω_{i} / Ω

is the vorticity unit vector. Vortex stretching (⁠

T_{1}

⁠) is the primary vorticity production mechanism in constant density incompressible flows. The vortex stretching term vanishes in two-dimensional flows. Vorticity-dilatation (⁠

T_{2}

⁠) quantifies the dilatational effect on vorticity growth and is absent in constant density incompressible flows. The baroclinic torque (⁠

T_{3}

⁠) generated due to misaligned pressure and density gradients is the dominant mechanism of vorticity production in incompressible RTI flows. Terms

T_{4}

and

T_{5}

represent the viscous diffusion and variable density effects on vorticity transport.

Profiles of the individual budget terms are considered and comparisons are made for different Mach numbers, at different times. Figure 14 presents the relative magnitude of the different terms in the vorticity transport equation budget for the 2D cases. The notation

⟨ \cdot ⟩

denotes the averaging operator in the homogeneous directions (x₂ in 2D and x₂ and x₃ in 3D), shown below for an arbitrary quantity

T

⁠,

⟨ T ⟩ (x_{1}, t) = \frac{1}{L_{x_{2}} L_{x_{3}}} \int_{- L_{x_{3} / 2}}^{L_{x_{3} / 2}} \int_{- L_{x_{2} / 2}}^{L_{x_{2} / 2}} T (x_{1}, x_{2}, x_{3}, t) d x_{2} d x_{3} .

(24)

FIG. 14.

View large Download slide

Vorticity budget contributions for cases with Re_p = 8000 and M = 0.0855 in (a)–(c), M = 0.2 in (d)–(f), and M = 0.3 in (g)–(i).

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.

To explore the effect of baroclinic torque on vorticity dynamics, we decompose the total pressure and density fields into mean and fluctuation contributions

p = ⟨ p ⟩ + p'; ρ = ⟨ ρ ⟩ + ρ' .

(25)

Note that the spatial derivatives of the mean quantities in these directions vanish. Then, the baroclinic torque can be decomposed as

T_{3} = \hat{ω_{i}} ϵ_{ijk} \frac{1}{A M^{2} ρ^{2}} (\underset{β_{1}}{\underset{︸}{\frac{\partial ρ'}{\partial x_{j}} \frac{\partial p'}{\partial x_{k}}}} + \underset{β_{2}}{\underset{︸}{\frac{\partial ⟨ p ⟩}{\partial x_{j}} \frac{\partial ρ'}{\partial x_{k}}}} + \underset{β_{3}}{\underset{︸}{\frac{\partial p'}{\partial x_{j}} \frac{\partial ⟨ ρ ⟩}{\partial x_{k}}}}) .

(26)

Here, β₁ represents the baroclinic torque generated due to non-linear interactions between the fluctuating density and pressure gradient. β₂ and β₃ 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 β₁ toward the baroclinic torque is negligible. For an isochoric background state, the β₃ 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

For a two dimensional flow, the three contributions arising from the baroclinic torque are

\begin{matrix} β_{1} = \frac{\hat{ω_{3}}}{A M^{2} ρ^{2}} (\frac{\partial ρ'}{\partial x_{2}} \frac{\partial p'}{\partial x_{1}} - \frac{\partial ρ'}{\partial x_{1}} \frac{\partial p'}{\partial x_{2}}), \\ β_{2} = \frac{\hat{ω_{3}}}{A M^{2} ρ^{2}} (\frac{\partial ρ'}{\partial x_{2}} \frac{\partial ⟨ p ⟩}{\partial x_{1}}), \\ β_{3} = - \frac{\hat{ω_{3}}}{A M^{2} ρ^{2}} (\frac{\partial p'}{\partial x_{2}} \frac{\partial ⟨ ρ ⟩}{\partial x_{1}}) . \end{matrix}

(27)

Here, β₁ is quadratic in perturbation; β₂ depends on the mean pressure gradient; and β₃ depends on the background density gradient.

In the linear limit, the mean pressure and density can be taken at the initial background state, which away from the interface is given by

⟨ p ⟩ = p_{0} (x_{1}) = 1 - M^{2} ρ_{0} x_{1}, ⟨ ρ ⟩ = ρ_{0} .

(28)

Under the linear approximation, the baroclinic term contributions can be simplified away from the interface as

\begin{matrix} T_{3} \approx \frac{\hat{ω_{3}}}{A M^{2} ρ_{0}^{2}} \frac{\partial p_{0}}{\partial x_{1}} \frac{\partial ρ'}{\partial x_{2}} \\ = {\begin{matrix} \frac{\hat{ω_{3}}}{A M^{2} ρ_{0}^{2}} (1 - M^{2} ρ_{h}) \frac{\partial ρ'}{\partial x_{2}} bubble side, \\ \frac{\hat{ω_{3}}}{A M^{2} ρ_{0}^{2}} (1 - M^{2} ρ_{l}) \frac{\partial ρ'}{\partial x_{2}} spike side . \end{matrix} \end{matrix}

(29)

The notation

\bar{.}

is used for averaging in the three spatial directions (⁠

x_{1}, x_{2}, x_{3}

⁠), shown below for an arbitrary quantity

T

⁠,

\begin{matrix} \bar{T} (t) = \frac{1}{L_{x_{1}} L_{x_{2}} L_{x_{3}}} \int_{- L_{x_{1} / 2}}^{L_{x_{1} / 2}} \int_{- L_{x_{2} / 2}}^{L_{x_{2} / 2}} \int_{- L_{x_{3} / 2}}^{L_{x_{3} / 2}} T (x_{1}, x_{2}, x_{3}, t) d x_{3} d x_{2} d x_{1} . \end{matrix}

(30)

We first compare the evolution of

\bar{β_{2}}

with that of

\bar{T_{3}}

⁠. Figure 15(a) shows that non-linear effects start to affect the evolution of

⟨ T_{3} ⟩

at around

τ = 1.4

⁠. Thus, at early stages, the linear approximation holds, which shows that the baroclinic term is likely larger on the spike side, due to the difference between ρ_h and ρ_l. To verify this using the numerical data, we focus on β₂, 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):

β_{2 s pike} = \frac{4}{λ L_{x_{1}}} \int_{λ / 4}^{λ / 2} \int_{- L_{x_{1}} / 2}^{L_{x_{1}} / 2} β_{2} d x_{1} d x_{2},

(31)

β_{2 b ubble} = \frac{4}{λ L_{x_{1}}} \int_{0}^{λ / 4} \int_{- L_{x_{1}} / 2}^{L_{x_{1}} / 2} β_{2} d x_{1} d x_{2} .

(32)

FIG. 15.

View large Download slide

Contrast between the evolutions of (a) $⟨ β_{2} ⟩$ and $⟨ T_{3} ⟩$ and (b) $β_{2 s pike}$ and $β_{2 b ubble}$ ⁠.

Figure 15(b) shows the evolutions of $β_{2 s pike}$ and $β_{2 b ubble}$ ⁠. It is seen that for $τ < 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² 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. β₂ is found to be the dominant term of the three contributions. β₂ has a bubble-spike asymmetry that is also indicative of the differential growth at early times. The non-linear term, β₁, shows a bubble-spike asymmetry but has an overall positive effect on the growth of vorticity. β₃ 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.

FIG. 16.

View large Download slide

Spatial variation of $⟨ β_{1} ⟩$ in (a)–(c), $⟨ β_{2} ⟩$ in (d)–(f) and $⟨ β_{3} ⟩$ in (g)–(i) at different times.

Figure 17 shows the growth of the three terms with time. At early times, β₂ 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¹⁴ 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.

FIG. 17.

Evolution of (a) β 1 ¯, (b) β 2 ¯, and (c) β 3 ¯ (the three contributions of the baroclinic torque) for three different Mach numbers for the 2D cases: Rep = 8000.

View large Download slide

Evolution of (a) $\bar{β_{1}}$ ⁠, (b) $\bar{β_{2}}$ ⁠, and (c) $\bar{β_{3}}$ (the three contributions of the baroclinic torque) for three different Mach numbers for the 2D cases: Re_p = 8000.

D. Kinetic energy evolution

The instability dynamics is considered in greater detail by examining the mean and turbulent kinetic energies. The velocity field is decomposed into a Favre average and fluctuation:⁴⁹

u_{i} = \tilde{u_{i}} + u_{i}^{″} .

(33)

The total kinetic energy

K = \frac{1}{2} u_{i} u_{i}

can then be partitioned into the mean, K_m, and turbulent kinetic, k, energies,

K_{m} = \frac{1}{2} \bar{ρ} \tilde{u_{i}} \tilde{u_{i}}; k = \frac{1}{2} ⟨ ρ u_{i}^{″} u_{i}^{″} ⟩ .

(34)

Figure 18(a) shows the evolution of turbulent kinetic energy,

\bar{k}

⁠, at different Mach numbers,

\bar{k} = \frac{1}{L_{x_{1}}} \int_{- L_{x_{1} / 2}}^{L_{x_{1} / 2}} k d x_{1} .

(35)

FIG. 18.

View large Download slide

Behavior of turbulent kinetic energy at different Mach numbers: (a) and (c) time evolution of $\bar{k}$ ⁠; (b) and (d) k variation across the layer at τ = 4, for 2D and 3D, respectively. (a) $\bar{k}$ evolution, (b) k profile, (c) $\bar{k}$ evolution, and (d) k profile.

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.

The normalized evolution equation for the mean kinetic energy is given as follows:^31,49

\begin{matrix} \frac{\partial K_{m}}{\partial t} + \frac{\partial K_{m} \tilde{u_{j}}}{\partial x_{j}} + \frac{\partial}{\partial x_{j}} [⟨ ρ u_{i}^{″} u_{j}^{″} \tilde{u_{i}} ⟩ + \frac{1}{A M^{2}} \bar{p} \tilde{u_{j}} - \frac{1}{R e} ⟨ τ_{i j} ⟩ \tilde{u_{i}}] \\ = R_{i j} \frac{\partial \tilde{u_{j}}}{\partial x_{j}} + \frac{1}{A M^{2}} \bar{p} \frac{\partial \tilde{u_{k}}}{\partial x_{k}} - \frac{1}{R e} ⟨ τ_{i j} ⟩ \frac{\partial \tilde{u_{i}}}{\partial x_{j}} - a_{i} \frac{1}{A M^{2}} \frac{\partial \bar{p}}{\partial x_{i}} \\ + \frac{1}{R e} a_{k} \frac{\partial ⟨ τ_{k j} ⟩}{\partial x_{j}} + \frac{1}{F r^{2}} \bar{ρ} \tilde{u_{i}} g δ_{1 i} . \end{matrix}

In RTI, the kinetic energy in the mean flow is generated due to the gravity work term (⁠

⟨ ρ u_{k} g δ_{1 k} ⟩

⁠).³² The normalized governing equation for fluctuating kinetic energy can be expressed in the following form:³¹^,

\begin{matrix} \frac{\partial k}{\partial t} + \frac{\partial (k \tilde{u_{j}})}{\partial x_{j}} + \frac{\partial}{\partial x_{j}} [\frac{1}{2} ⟨ ρ u_{i}^{″} u_{i}^{″} u_{j}^{″} ⟩ + ⟨ u_{i}^{'} (\frac{1}{A M^{2}} p' δ_{i j} - \frac{1}{R e} τ_{i j}^{'}) ⟩] = \underset{K 1}{\underset{︸}{- R_{i j} \frac{\partial \tilde{u_{j}}}{\partial x_{j}}}} + \underset{K 2}{\underset{︸}{\frac{1}{A M^{2}} ⟨ p' d ⟩}} + \underset{K 3}{\underset{︸}{\frac{1}{A M^{2}} a_{i} \frac{\partial \bar{p}}{\partial x_{i}}}} - \underset{K 4}{\underset{︸}{\frac{1}{R e} ⟨ τ_{i k}^{'} \frac{\partial u_{i}^{'}}{\partial x_{k}} ⟩}} - \underset{K 5}{\underset{︸}{\frac{1}{R e} a_{i} \frac{\partial ⟨ τ_{k j} ⟩}{\partial x_{j}}}} . \end{matrix}

(36)

Here,

a_{i} = - ⟨ u_{i}^{″} ⟩ = ⟨ ρ' u_{i}^{'} ⟩ / ⟨ ρ ⟩

is the normalized mass-flux,³¹

R_{i j} = ⟨ ρ u_{i}^{″} u_{j}^{″} ⟩

is the Favre averaged Reynolds stresses and

d = \frac{\partial u_{i}^{'}}{\partial 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.

FIG. 19.

View large Download slide

Turbulent kinetic energy budget across the layer for 2D cases, with Mach numbers: (a) M = 0.0855, (b) M = 0.2, and (c) M = 0.3, and Re_p = 8000 at τ = 4.

FIG. 20.

View large Download slide

Turbulent kinetic energy budget across the layer for 3D cases, with Mach numbers: (a) M = 0.0855 and (b) M = 0.3, and Re_p = 1000 at τ = 4.

Figures 21 and 22 show the profiles of mass flux in vertical direction (a₁) 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.

FIG. 21.

Vertical mass flux (a1) profiles for 2D, with Mach numbers: (a) M = 0.0855, (b) M = 0.2, and (c) M = 0.3, and Rep = 8000 at τ = 4.

View large Download slide

Vertical mass flux (a₁) profiles for 2D, with Mach numbers: (a) M = 0.0855, (b) M = 0.2, and (c) M = 0.3, and Re_p = 8000 at τ = 4.

FIG. 22.

Mean pressure gradient ( ∂ ⟨ p ⟩ ∂ x 1) profiles for 2D cases, with Mach numbers: (a) M = 0.0855, (b) M = 0.2, and (c) M = 0.3, and Rep = 8000 at τ = 4.

View large Download slide

Mean pressure gradient (⁠ $\frac{\partial ⟨ p ⟩}{\partial x_{1}}$ ⁠) profiles for 2D cases, with Mach numbers: (a) M = 0.0855, (b) M = 0.2, and (c) M = 0.3, and Re_p = 8000 at τ = 4.

FIG. 23.

Vertical mass flux (a1) profiles for 3D cases, with Mach numbers: (a) M = 0.0855 and (b) M = 0.3, and Rep = 1000 at τ = 4.

View large Download slide

Vertical mass flux (a₁) profiles for 3D cases, with Mach numbers: (a) M = 0.0855 and (b) M = 0.3, and Re_p = 1000 at τ = 4.

FIG. 24.

Mean pressure gradient ( ∂ ⟨ p ⟩ ∂ x 1) profiles for 3D cases, with Mach numbers: (a) M = 0.0855 and (b) M = 0.3, and Rep = 1000 at τ = 4.

View large Download slide

Mean pressure gradient (⁠ $\frac{\partial ⟨ p ⟩}{\partial x_{1}}$ ⁠) profiles for 3D cases, with Mach numbers: (a) M = 0.0855 and (b) M = 0.3, and Re_p = 1000 at τ = 4.

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 (β₁), mean pressure-gradient effect (β₂), and mean density-gradient effect (β₃). It is found that the contribution of β₂ component is most significant toward vorticity growth, especially at early times. At early times, the β₂ contribution can be directly expressed using the background pressure variation and exhibits different magnitudes on the spike and bubble sides, proportional to ρ_h and ρ_l, respectively. This difference inherently initiates the bubble-spike asymmetry seen at later stages. Furthermore, since β₂ 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 β₂ 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

1.

L.

Rayleigh

, “

Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density

,”

Proc. London Math. Soc.

s1–14

,

170

–

177

(

1882

).

https://doi.org/10.1112/plms/s1-14.1.170

Google Scholar

Crossref

2.

G. I.

Taylor

, “

The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I

,”

Proc. R. Soc. London, Ser. A

201

,

192

–

196

(

1950

).

https://doi.org/10.1098/rspa.1950.0052

Google Scholar

Crossref

3.

G.

Boffetta

and

A.

Mazzino

, “

Incompressible Rayleigh–Taylor turbulence

,”

Annu. Rev. Fluid Mech.

49

,

119

–

143

(

2017

).

https://doi.org/10.1146/annurev-fluid-010816-060111

Google Scholar

Crossref

4.

D.

Livescu

, “

Turbulence with large thermal and compositional density variations

,”

Annu. Rev. Fluid Mech.

52

,

309

–

341

(

2020

).

https://doi.org/10.1146/annurev-fluid-010719-060114

Google Scholar

Crossref

5.

R. S.

Craxton

,

K. S.

Anderson

,

T. R.

Boehly

,

V. N.

Goncharov

,

D. R.

Harding

,

J. P.

Knauer

,

R. L.

McCrory

,

P. W.

McKenty

,

D. D.

Meyerhofer

,

J. F.

Myatt

,

A. J.

Schmitt

,

J. D.

Sethian

,

R. W.

Short

,

S.

Skupsky

,

W.

Theobald

,

W. L.

Kruer

,

K.

Tanaka

,

R.

Betti

,

T. J. B.

Collins

,

J. A.

Delettrez

,

S. X.

Hu

,

J. A.

Marozas

,

A. V.

Maximov

,

D. T.

Michel

,

P. B.

Radha

,

S. P.

Regan

,

T. C.

Sangster

,

W.

Seka

,

A. A.

Solodov

,

J. M.

Soures

,

C.

Stoeckl

, and

J. D.

Zuegel

, “

Direct-drive inertial confinement fusion: A review

,”

Phys. Plasmas

22

,

110501

(

2015

).

https://doi.org/10.1063/1.4934714

Google Scholar

Crossref

6.

Y.

Zhou

, “

Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. I

,”

Phys. Rep.

720–722

,

1

–

136

(

2017

).

https://doi.org/10.1016/j.physrep.2017.07.005

Google Scholar

Crossref

7.

Y.

Zhou

, “

Rayleigh-Taylor and Richtmyer-Meshkov instability induced flow, turbulence, and mixing. II

,”

Phys. Rep.

723–725

,

1

–

160

(

2017

).

https://doi.org/10.1016/j.physrep.2017.07.008

Google Scholar

Crossref

8.

S. F.

Gull

, “

The x-ray, optical and radio properties of young supernova remnants

,”

Mon. Not. R. Astron. Soc.

171

,

263

–

278

(

1975

).

https://doi.org/10.1093/mnras/171.2.263

Google Scholar

Crossref

9.

A. S.

Colgate

, “

Prompt gamma rays and x rays from supernovae

,”

Can. J. Phys.

46

,

S476

–

S480

(

1968

).

https://doi.org/10.1139/p68-274

Google Scholar

Crossref

10.

M.

Zingale

,

S. E.

Woosley

,

C. A.

Rendleman

,

M. S.

Day

, and

J. B.

Bell

, “

Three-dimensional numerical simulations of Rayleigh-Taylor unstable flames in Type Ia supernovae

,”

Astrophys. J.

632

,

1021

(

2005

).

https://doi.org/10.1086/433164

Google Scholar

Crossref

11.

A.

Nouri

,

P.

Givi

, and

D.

Livescu

, “

Modeling and simulation of turbulent nuclear flames in type Ia supernovae

,”

Prog. Aerosp. Sci.

108

,

156

–

179

(

2019

).

https://doi.org/10.1016/j.paerosci.2019.04.004

Google Scholar

Crossref

12.

C.-S.

Huang

,

M. C.

Kelley

, and

D. L.

Hysell

, “

Nonlinear Rayleigh-Taylor instabilities, atmospheric gravity waves and equatorial spread F

,”

J. Geophys. Res.: Space Phys.

98

,

15631

–

15642

, https://doi.org/10.1029/93JA00762 (

1993

).

Google Scholar

Crossref

13.

D.

Livescu

, “

Numerical simulations of two-fluid turbulent mixing at large density ratios and applications to the Rayleigh-Taylor instability

,”

Philos. Trans. R. Soc. A

371

,

20120185

(

2013

).

https://doi.org/10.1098/rsta.2012.0185

Google Scholar

Crossref

14.

S.

Weiland

,

P. E.

Hamlington

,

S. J.

Reckinger

, and

D.

Livescu

, “

Effects of isothermal stratification strength on vorticity dynamics for single-mode compressible Rayleigh-Taylor instability

,”

Phys. Rev. Fluids

4

,

093905

(

2019

).

https://doi.org/10.1103/PhysRevFluids.4.093905

Google Scholar

Crossref

15.

X.

Bian

,

H.

Aluie

,

D.

Zhao

,

H.

Zhang

, and

D.

Livescu

, “

Revisiting the late-time growth of single-mode Rayleigh–Taylor instability and the role of vorticity

,”

Physica D

403

,

132250

(

2020

).

https://doi.org/10.1016/j.physd.2019.132250

Google Scholar

Crossref

16.

T.

Luo

,

J.

Wang

,

C.

Xie

,

M.

Wan

, and

S.

Chen

, “

Effects of compressibility and Atwood number on the single-mode Rayleigh-Taylor instability

,”

Phys. Fluids

32

,

012110

(

2020

).

https://doi.org/10.1063/1.5131585

Google Scholar

Crossref

17.

Z.-X.

Hu

,

Y.-S.

Zhang

,

B.

Tian

,

Z.

He

, and

L.

Li

, “

Effect of viscosity on two-dimensional single-mode Rayleigh-Taylor instability during and after the reacceleration stage

,”

Phys. Fluids

31

,

104108

(

2019

).

https://doi.org/10.1063/1.5122247

Google Scholar

Crossref

18.

S.

Gauthier

, “

Compressible Rayleigh–Taylor turbulent mixing layer between Newtonian miscible fluids

,”

J. Fluid Mech.

830

,

211

–

256

(

2017

).

https://doi.org/10.1017/jfm.2017.565

Google Scholar

Crossref

19.

A. G. W.

Lawrie

and

S. B.

Dalziel

, “

Rayleigh–Taylor mixing in an otherwise stable stratification

,”

J. Fluid Mech.

688

,

507

–

527

(

2011

).

https://doi.org/10.1017/jfm.2011.398

Google Scholar

Crossref

20.

S. D. M.

Wykes

and

S. B.

Dalziel

, “

Efficient mixing in stratified flows: Experimental study of a Rayleigh–Taylor unstable interface within an otherwise stable stratification

,”

J. Fluid Mech.

756

,

1027

–

1057

(

2014

).

https://doi.org/10.1017/jfm.2014.308

Google Scholar

Crossref

21.

D.

Livescu

, “

Compressibility effects on the Rayleigh–Taylor instability growth between immiscible fluids

,”

Phys. Fluids

16

,

118

–

127

(

2004

).

https://doi.org/10.1063/1.1630800

Google Scholar

Crossref

22.

S.

Weiland

,

S. J.

Reckinger

,

P. E.

Hamlington

, and

D.

Livescu

, “

Effects of background stratification on the compressible Rayleigh Taylor instability

,” AIAA Paper No. 2017-3974,

2017

.

Google Scholar

23.

S.

Gerashchenko

and

D.

Livescu

, “

Viscous effects on the Rayleigh-Taylor instability with background temperature gradient

,”

Phys. Plasmas

23

,

072121

(

2016

).

https://doi.org/10.1063/1.4959810

Google Scholar

Crossref

24.

T.

Wei

and

D.

Livescu

, “

Late-time quadratic growth in single-mode Rayleigh-Taylor instability

,”

Phys. Rev. E

86

,

046405

(

2012

).

https://doi.org/10.1103/PhysRevE.86.046405

Google Scholar

Crossref

25.

H.

Liang

,

Z.

Xia

, and

H.

Huang

, “

Late-time description of immiscible Rayleigh–Taylor instability: A lattice Boltzmann study

,”

Phys. Fluids

33

,

082103

(

2021

).

https://doi.org/10.1063/5.0057269

Google Scholar

Crossref

26.

A.

Hamzehloo

,

P.

Bartholomew

, and

S.

Laizet

, “

Direct numerical simulations of incompressible Rayleigh–Taylor instabilities at low and medium Atwood numbers

,”

Phys. Fluids

33

,

054114

(

2021

).

https://doi.org/10.1063/5.0049867

Google Scholar

Crossref

27.

H.

Liang

,

X.

Hu

,

X.

Huang

, and

J.

Xu

, “

Direct numerical simulations of multi-mode immiscible Rayleigh-Taylor instability with high Reynolds numbers

,”

Phys. Fluids

31

,

112104

(

2019

).

https://doi.org/10.1063/1.5127888

Google Scholar

Crossref

28.

P.

Ramaprabhu

,

G.

Dimonte

,

Y. N.

Young

,

A. C.

Calder

, and

B.

Fryxell

, “

Limits of the potential flow approach to the single-mode Rayleigh-Taylor problem

,”

Phys. Rev. E

74

,

066308

(

2006

).

https://doi.org/10.1103/PhysRevE.74.066308

Google Scholar

Crossref

29.

G.

Dimonte

et al, “

A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations: The Alpha-Group collaboration

,”

Phys. Fluids

16

,

1668

–

1692

(

2004

).

https://doi.org/10.1063/1.1688328

Google Scholar

Crossref

30.

S. J.

Reckinger

,

D.

Livescu

, and

O. V.

Vasilyev

, “

Comprehensive numerical methodology for direct numerical simulations of compressible Rayleigh–Taylor instability

,”

J. Comput. Phys.

313

,

181

–

208

(

2016

).

https://doi.org/10.1016/j.jcp.2015.11.002

Google Scholar

Crossref

31.

D.

Livescu

,

J. R.

Ristorcelli

,

R. A.

Gore

,

S. H.

Dean

,

W. H.

Cabot

, and

A. W.

Cook

, “

High-Reynolds number Rayleigh–Taylor turbulence

,”

J. Turbul.

10

,

N13

(

2009

).

https://doi.org/10.1080/14685240902870448

Google Scholar

Crossref

32.

D.

Livescu

and

J. R.

Ristorcelli

, “

Buoyancy-driven variable-density turbulence

,”

J. Fluid Mech.

591

,

43

–

71

(

2007

).

https://doi.org/10.1017/S0022112007008270

Google Scholar

Crossref

33.

D.

Aslangil

,

D.

Livescu

, and

A.

Banerjee

, “

Atwood and Reynolds numbers effects on the evolution of buoyancy-driven homogeneous variable-density turbulence

,”

J. Fluid Mech.

895

,

A12

(

2020

).

https://doi.org/10.1017/jfm.2020.268

Google Scholar

Crossref

34.

K.

Xu

, “

A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method

,”

J. Comput. Phys.

171

,

289

–

335

(

2001

).

https://doi.org/10.1006/jcph.2001.6790

Google Scholar

Crossref

35.

K.

Xu

, “

Gas-kinetic schemes for unsteady compressible flow simulations

,” 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

36.

K.

Xu

,

M.

Mao

, and

L.

Tang

, “

A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow

,”

J. Comput. Phys.

203

,

405

–

421

(

2005

).

https://doi.org/10.1016/j.jcp.2004.09.001

Google Scholar

Crossref

37.

K.

Xu

,

L.

Martinelli

, and

A.

Jameson

, “

Gas-kinetic finite volume methods, flux-vector splitting, and artificial diffusion

,”

J. Comput. Phys.

120

,

48

–

65

(

1995

).

https://doi.org/10.1006/jcph.1995.1148

Google Scholar

Crossref

38.

P. L.

Bhatnagar

,

E. P.

Gross

, and

M.

Krook

, “

A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems

,”

Phys. Rev.

94

,

511

–

525

(

1954

).

https://doi.org/10.1103/PhysRev.94.511

Google Scholar

Crossref

39.

G.

Kumar

,

S. S.

Girimaji

, and

J.

Kerimo

, “

WENO-enhanced gas-kinetic scheme for direct simulations of compressible transition and turbulence

,”

J. Comput. Phys.

234

,

499

–

523

(

2013

).

https://doi.org/10.1016/j.jcp.2012.10.005

Google Scholar

Crossref

40.

S. K.

Lele

, “

Compact finite difference schemes with spectral-like resolution

,”

J. Comput. Phys.

103

,

16

–

42

(

1992

).

https://doi.org/10.1016/0021-9991(92)90324-R

Google Scholar

Crossref

41.

J.

Larsson

and

S. K.

Lele

, “

Direct numerical simulation of canonical shock/turbulence interaction

,”

Phys. Fluids

21

,

126101

(

2009

).

https://doi.org/10.1063/1.3275856

Google Scholar

Crossref

42.

J.

Ryu

and

D.

Livescu

, “

Turbulence structure behind the shock in canonical shock–vortical turbulence interaction

,”

J. Fluid Mech.

756

,

R1

(

2014

).

https://doi.org/10.1017/jfm.2014.477

Google Scholar

Crossref

43.

J.

Larsson

,

I.

Bermejo-Moreno

, and

S. K.

Lele

, “

Reynolds- and Mach-number effects in canonical shock–turbulence interaction

,”

J. Fluid Mech.

717

,

293

–

321

(

2013

).

https://doi.org/10.1017/jfm.2012.573

Google Scholar

Crossref

44.

M. L.

Wong

,

D.

Livescu

, and

S. K.

Lele

, “

High-resolution Navier–Stokes simulations of Richtmyer–Meshkov instability with reshock

,”

Phys. Rev. Fluids

4

,

104609

(

2019

).

https://doi.org/10.1103/PhysRevFluids.4.104609

Google Scholar

Crossref

45.

N.

Peng

,

Y.

Yang

,

J.

Wu

, and

Z.

Xiao

, “

Mechanism and modelling of the secondary baroclinic vorticity in the Richtmyer–Meshkov instability

,”

J. Fluid Mech.

911

,

A56

(

2021

).

https://doi.org/10.1017/jfm.2020.1080

Google Scholar

Crossref

46.

G.

Peng

,

N. J.

Zabusky

, and

S. S.

Zhang

, “

Vortex-accelerated secondary baroclinic vorticity deposition and late-intermediate time dynamics of a two-dimensional Richtmyer–Meshkov interface

,”

Phys. Fluids

15

,

3730

–

3744

(

2003

).

https://doi.org/10.1063/1.1621628

Google Scholar

Crossref

47.

A. M.

Steinberg

,

J. F.

Driscoll

, and

S. L.

Ceccio

, “

Three-dimensional temporally resolved measurements of turbulence-flame interactions using orthogonal-plane cinema-stereoscopic PIV

,”

Exp. Fluids

47

,

527

–

547

(

2009

).

https://doi.org/10.1007/s00348-009-0677-7

Google Scholar

Crossref

48.

P. E.

Hamlington

,

A. Y.

Poludnenko

, and

E. S.

Oran

, “

Interactions between turbulence and flames in premixed reacting flows

,”

Phys. Fluids

23

,

125111

(

2011

).

https://doi.org/10.1063/1.3671736

Google Scholar

Crossref

49.

A.

Mittal

and

S. S.

Girmaji

, “

Mathematical framework for analysis of internal energy dynamics and spectral distribution in compressible turbulent flows

,”

Phys. Rev. Fluids

4

,

042601

(

2019

).

https://doi.org/10.1103/PhysRevFluids.4.042601

Google Scholar

Crossref

2023

Author(s)

Compressible Rayleigh–Taylor instability subject to isochoric initial background state

I. INTRODUCTION

II. FLUID GOVERNING EQUATIONS

A. Background state

III. COMPUTATIONAL APPROACH

IV. RAYLEIGH–TAYLOR INSTABILITY-DRIVEN COMPRESSIBLE MIXING

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

B. Three-dimensional effects

C. Vorticity evolution

1. Linear analysis for the early time behavior

2. Non-linear analysis

D. Kinetic energy evolution

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Compressible Rayleigh–Taylor instability subject to isochoric initial background state

I. INTRODUCTION

II. FLUID GOVERNING EQUATIONS

A. Background state

III. COMPUTATIONAL APPROACH

IV. RAYLEIGH–TAYLOR INSTABILITY-DRIVEN COMPRESSIBLE MIXING

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

B. Three-dimensional effects

C. Vorticity evolution

1. Linear analysis for the early time behavior

2. Non-linear analysis

D. Kinetic energy evolution

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only