Investigation of strong isothermal stratification effects on multi-mode compressible Rayleigh–Taylor instability

Aslangil, Denis; Wong, Man Long

doi:10.1063/5.0164504

Rayleigh–Taylor instability, RTI, occurs at the interface separating two fluids subjected to acceleration when the density gradient and the acceleration are in opposite directions. Previous scientific research primarily considered RTI under the incompressible assumption, which may not be valid in many high-energy-density engineering applications and astrophysical phenomena. In this study, the compressibility effects of the background isothermal stratification strength on multi-mode two-dimensional RTI are explored using fully compressible multi-species direct numerical simulations. Cases under three different isothermal Mach numbers – $M a = 0.15, 0.3, and 0.45$ – are investigated to explore weakly, moderately, and strongly stratified compressible RTI, respectively, at an Atwood number of 0.04. Unlike incompressible RTI, an increase in the flow compressibility through the strength of the background stratification can suppress the RTI growth and can lead to a termination of the RTI mixing layer growth with a highly molecularly mixed state. Our findings suggest that even at the chosen relatively low Atwood number, the variable-density effects can be significantly enhanced due to an increase in the background stratification for the compressible RTI as different spatial profiles become noticeably asymmetric across the mixing layer for the strongly stratified case. In addition, this study compares the chaotic behavior of the cases by studying the transport of the turbulent kinetic energy as well as the vortex dynamics. The Reynolds number dependence of the results is also examined with three different Reynolds numbers, and the findings for the large-scale mixing and flow quantities of interest are shown to be universal in the range of the Reynolds numbers studied.

I. INTRODUCTION

The Rayleigh–Taylor (RT) instability, or RTI, is an instability at the interface between two fluids subjected to acceleration when the directions of the density gradient and the acceleration have opposite signs.^1,2 The RTI is observed in nature from atmospheric, oceanic, and geologic flows to astrophysical phenomena, such as supernova explosions. It also occurs in engineering applications with high-energy-density (HED) processes such as inertial confinement fusion (ICF).³ As a result, the RTI is of great interest in the mixing community, and it has been the subject of many experimental and numerical research studies.^4–8 Although the flow instability is observed between the highly compressible flows in many natural phenomena and engineering systems, most of the numerical studies benefited from using the incompressible assumption to reduce computational expenses. In this study, we present the flow compressibility effects due to the background stratification on the multi-species two-dimensional (2D) multi-mode RTI with results from direct numerical simulations (DNS).

Using the linear stability analysis, RTI has exponential growth at early times when the initial perturbation is small. At later times, the height of the multi-mode incompressible RTI mixing layer, h, under constant acceleration has self-similar growth that is defined as⁹^,

h \approx α A g t^{2},

(1)

where t is the time, g is the acceleration, α is the growth rate parameter, and A is the Atwood number representing the density contrast between the two mixing fluids. By assuming that both fluids are ideal gases, A is defined as

A = \frac{W_{h} - W_{l}}{W_{h} + W_{l}} = \frac{ρ_{h} - ρ_{l}}{ρ_{h} + ρ_{l}},

(2)

where W_h and W_l are the molecular weights of the heavy and light fluids, respectively, while ρ_h and ρ_l are the initial species densities of the heavy and light fluids. However, Eq. may not be applicable for approximating the RTI growth, which is observed in HED applications and astrophysics phenomena. For example, when the acceleration is not constant and varies as a function of time, the validity of Eq. is still an open question.^10–13 The equation is also derived based on a similarity analysis that ignores any compressibility effects, so the validity of the equation is lost in the compressible regime.¹⁴

There have been several investigations on single-mode and multi-mode compressible RTI in previous literature. Wieland et al.¹⁵ carried out 2D single-mode compressible RTI DNS using adaptive wavelet-based mesh refinement. It was shown that an increase in the background stratification strength, which determines the flow compressibility of the mixing fluids, leads to the ceasing of the mixing layer growth for the single-mode RTI. Simulations of 2D single-mode compressible RTI with compact finite difference scheme and hyper-viscosity model were also performed by Luo et al.¹⁶ to examine the stabilizing and destabilizing effects of the instability under different Mach and Atwood numbers. The work is further extended to study the effects of the Atwood number and stratification parameter on 2D compressible RTI with multi-mode perturbations.¹⁷ Mellado et al.¹⁸ conducted three-dimensional (3D) large-eddy simulations of compressible RTI and found that intrinsic compressibility effects were limited as the turbulent Mach number is bounded in the range of 0.25–0.6. A DNS of 3D multi-mode compressible RT turbulence at a moderate Atwood number, A = 0.25, was conducted using the auto-adaptive multidomain Chebyshev–Fourier–Fourier numerical method by Gauthier,¹⁹ and the simulation results were compared to cases with either a lower Reynolds number or a lower Atwood number within the Boussinesq approximation. In a recent work by Luo and Wang,²⁰ the mixing and energy transfer of 3D compressible RT turbulence indicates that the direct and reverse subgrid-scale (SGS) fluxes can be strengthened by compression and expansion motions, respectively.

Multi-species mixing is also of great interest in the scientific community as it is observed in broader engineering applications, and it has been shown that the flow dynamics could be significantly different than single fluid flows.²¹ In addition to RTI, there is also a wide study on multi-species mixing induced by the Richtmyer–Meshkov instability (RMI),^22–24 which is an interfacial instability similar to RTI but induced by impulsive accelerations, such as shock waves. While the mixing of fluids induced by RMI is commonly studied with compressible simulations due to the compressibility effects of shock waves, high-resolution RMI simulations²⁵ show that under a normal configuration, the flows within the RMI mixing layers are quasi-incompressible except around times of shock–interface interactions. The incompressible multi-species variable-density mixing dynamics of the turbulence aiming to mimic the core regions of the RTI and RMI mixing layers have been studied using several different homogeneous flow configurations as well.^21,26–30 These studies lack information about compressibility effects, such as those from the background stratification, which can potentially play a significant role in the multi-species mixing processes.

For the compressible RTI, instability and multi-species mixing can be largely affected by many factors, such as background stratification, fluid properties, and acoustic effects. The mixing process can also generate strong shock waves through the piston-like motions of the bubbles and spikes.³¹ The growth of compressible RTI can be affected either through fluid compressibility or flow compressibility.³² Fluid compressibility is governed by the fluids' properties, such as the ratios of specific heat capacities, while flow compressibility is associated with the thermodynamic states independent of the fluid properties, such as the interface pressure. The former can be considered as dynamic or intrinsic, which is an effect of the speed of sound; and the latter can be viewed as static, as this compressibility effect results in a stable stratification.^19,33

The effects of compressibility originating from the flow compressibility are investigated in this work using 2D DNS of multi-mode RTI with different isothermal background stratification strengths, where some of the cases are strongly stratified. The flow compressibility can be controlled by the initial background stratification, which is characterized by the isothermal Mach number, Ma,^15,34 defined as

M a = \frac{\sqrt{g λ_{0}}}{c_{0}},

(3)

where

c_{0} = \sqrt{p^{I} / ρ^{I}}

is the initial isothermal speed of sound defined by the interface pressure, p^I, and interface density, ρ^I, and λ₀ is the characteristic length scale of the problem. This isothermal Mach number definition decouples the fluid compressibility characterized by the values of the ratios of the specific heats and the flow compressibility characterized by the thermodynamic state of the system as explained in the work by Reckinger et al.³⁴ It should be pointed out that when the interface pressure tends to infinity, i.e.,

p^{I} \to \infty

⁠, the RTI approaches the incompressible limit as

M a \to 0

⁠.^3,32,35 The background stratification strength is also commonly represented by the stratification parameter, Sr, in previous literature,^16,17,19,20 where

S r = M a^{2}

if the same characteristic length scale is used.

In many previous papers on multi-mode RTI, the domain width is commonly used as the characteristic length scale in the definition of the isothermal Mach number [see Eq. (3) and/or stratification parameter]. However, the characteristic wavelength of the initial perturbation is a better length scale to characterize the RTI problem, at least during the early times of the flow evolution, which is similar to the incompressible RTI study by Livescu et al.³⁶ Hence, in the definition of Ma, we choose to use the characteristic wavelength of the perturbations instead of the domain width, which is different than many previous multi-mode compressible RTI studies. For example, the largest background stratification 3D RTI multi-mode cases have Sr = 6 reported based on the domain width in the work by Gauthier.¹⁹ This, in fact, corresponds to around $S r \approx 0.04$ or $M a \approx 0.2$ based on the characteristic wavelength (assuming the problem domain width is $2 π$ as they used a pseudo-spectral code). Similarly, Sr = 0.1 − Sr = 1.0 based on the shown domain width in the work by Luo and Wang¹⁷ is equivalent to $S r \approx 0.013$ − $S r \approx 0.13$ ⁠, or $M a \approx 0.11$ − $M a \approx 0.35$ ⁠, if the characteristic wavelength of the perturbations is used as the reference length. In addition, Luo and Wang²⁰ extended their 2D work to 3D for their moderate Atwood number A = 0.5 by performing numerical simulations up to Sr = 3.0 based on a domain width that corresponds to around eight times smaller stratification parameter when the characteristic wavelength is used in the definition of Sr. Our study on the 2D multi-mode RTI includes cases with Ma = 0.45 where the background stratification strength is stronger than the problems in many previous multi-mode compressible RTI studies.

Our choice of the characteristic length scale can also provide a more direct and fair comparison between our multi-mode compressible RTI results with the single-mode compressible RTI results in previous works, such as the Ma = 0.3 case in the work by Wieland et al.¹⁵ To confirm that there is no re-growth at the late times of the flow evolution for the moderately and strongly stratified cases, the ensemble-averaged statistical results of the multi-mode cases in this study are obtained with DNS performed up to much longer non-dimensional times compared to the single-mode cases in the aforementioned work. The findings of this study suggest that in the case with a strong background stratification, Eq. (1) loses its validity even for the RTI under constant acceleration.

The paper is structured as follows: in Sec. II, we detail the governing equations and numerical methods, followed by the setup of DNS cases in Sec. III. Then, we present the analysis of the dynamics of the compressible mixing in Sec. IV. Section V shows and discusses the flow energetics and small-scale vortical motions. Finally, in Sec. VI, we summarize the results of this study.

II. PROBLEM DESCRIPTION

A. Governing equations

The governing equations are the compressible multi-species Navier–Stokes equations

\frac{\partial ρ Y_{i}}{\partial t} + \nabla \cdot (ρ u Y_{i}) = - \nabla \cdot J_{i},

(4)

\frac{\partial ρ u}{\partial t} + \nabla \cdot (ρ uu + p δ) = \nabla \cdot τ + ρ g,

(5)

\begin{matrix} \frac{\partial E}{\partial t} + \nabla \cdot [(E + p) u] = \nabla \cdot (τ \cdot u - q_{c} - q_{d}) + ρ g \cdot u, \end{matrix}

(6)

where ρ is the mixture density,

u = {[u, v]}^{T} = {[u_{1}, u_{2}]}^{T}

is the velocity vector, p is the pressure, and E is the total energy per unit volume of the fluid mixture. Y_i is the mass fraction of species

i \in {1, 2, \dots, N}

⁠, with N the total number of species. All Y_i sum up to one by definition, and g is the body force vector per unit mass. In this work,

g = {[g, 0]}^{T}

is assumed, where g is a constant parameter.

J_{i}

is the diffusive mass flux for species i.

τ, q_{c}

⁠, and

q_{d}

⁠, which, respectively, are viscous stress tensor, conductive heat flux, and inter-species diffusional enthalpy flux.

δ

is the identity tensor.

The mixture is assumed to be ideal and calorically perfect, with

E = ρ (e + \frac{1}{2} u \cdot u),

(7)

p = (γ - 1) ρ e, e = c_{v} T,

(8)

where e and T are the mixture specific internal energy and temperature; and γ and c_v are the ratio of specific heats and specific heat at constant volume of the mixture.

In this work, a binary mixture with N = 2 is studied for each case. The multi-component diffusive mass flux of species i is given by Fick's law

J_{i} = - ρ D_{i} \nabla Y_{i}, i = 1, 2,

(9)

where

D_{1} = D_{2} = D

is the binary diffusion coefficient. The baro-diffusion and thermo-diffusion are not considered in this work.

The mixture is also assumed to be Newtonian with viscous stress tensor,

τ

⁠, given by

τ = 2 μ S + (μ_{v} - \frac{2}{3} μ) δ (\nabla \cdot u),

(10)

where μ and μ_v are the shear viscosity and bulk viscosity of the mixture, respectively. The effects of the bulk viscosity are neglected in this work. S is the strain-rate tensor given by

S = \frac{1}{2} [\nabla u + {(\nabla u)}^{T}] .

(11)

The conductive flux and the inter-species diffusional enthalpy flux³⁷ are given by

q_{c} = - κ \nabla T,

(12)

q_{d} = \sum_{i = 1}^{N} h_{i} J_{i},

(13)

where κ is the thermal conductivity of the mixture. h_i is the specific enthalpy of species i

h_{i} = c_{p, i} T,

(14)

where

c_{p, i}

is the specific heat at constant pressure of species i.

Assuming that all species are at pressure and temperature equilibria, the mixing rules for the fluid properties γ, c_v, μ, and κ are given by

γ = \frac{c_{p}}{c_{v}}, c_{p} = \sum_{i = 1}^{N} Y_{i} c_{p, i}, c_{v} = \sum_{i = 1}^{N} Y_{i} c_{v, i},

(15)

μ = \frac{\sum_{i = 1}^{N} μ_{i} Y_{i} / \sqrt{W_{i}}}{\sum_{i = 1}^{N} Y_{i} / \sqrt{W_{i}}}, κ = \frac{\sum_{i = 1}^{N} κ_{i} Y_{i} / \sqrt{W_{i}}}{\sum_{i = 1}^{N} Y_{i} / \sqrt{W_{i}}},

(16)

where

c_{v, i}

⁠, μ_i, κ_i, and W_i are the specific heat at constant volume, dynamic viscosity, thermal conductivity, and molecular weight of species i, respectively. The specific gas constant R_i of species i is given by

R_{i} = \frac{R_{g}}{W_{i}},

(17)

where R_g is the universal gas constant. The ratios of specific heats of both species, γ_i, are assumed to be the same at the value of 1.4 in this study. Thus, the mixture has a constant and uniform mixture γ at the value of 1.4. In addition, the heavier fluid with larger molecular weight is assigned as the first species, with i = 1.

B. Non-dimensional problem parameters

In the initial setup of the RTI problem, the heavier and the lighter fluids have initial species densities

ρ_{1, 0}

and

ρ_{2, 0}

⁠, respectively, at the interface, and the pressure is continuous at the interface with the value p^I. The kinematic viscosity, ν_i, of both species are initially assumed to be identical at the interface, i.e.,

ν_{i, 0} = ν_{0}

⁠, while the dynamic viscosity of each species is assumed to be constant in time and uniform in space, which is given by

μ_{i} = ρ_{i, 0} ν_{0} .

(18)

The RTI problem can be described by a few important non-dimensional parameters. First of all, the initial Atwood number of A = 0.04 [defined by Eq. ] is chosen for all cases in this study. The problem Reynolds number, Re₀, is another parameter used for the initialization of the problem and is given by

R e_{0} = \sqrt{\frac{g λ_{0}^{3}}{ν_{0}^{2}}},

(19)

where λ₀ is the characteristic wavelength of the initial perturbation. Three different problem Reynolds numbers –

R e_{0} = 3187.5, R e_{0} = 6375.0

⁠, and

R e_{0} = 12 750.0

– are chosen for the multi-mode cases with different background stratification strengths. In this work, the isothermal Mach number,

M a = \sqrt{g λ_{0}} / c_{0}

[also defined by Eq. but repeated here], that characterizes the background stratification strength, is chosen to be 0.15, 0.30, and 0.45, which represent the weakly stratified, moderately stratified, and strongly stratified cases, respectively.

The fluid properties are described by the Schmidt number, Sc, between the two fluids and the species Prandtl number, Pr_i,

S c = \frac{ν_{0}}{D}, P r_{i} = \frac{c_{p, i} μ_{i}}{κ_{i}},

(20)

where Sc = 1,

P r_{1} = 0.92

⁠, and

P r_{2} = 1.08

for all cases in this work.

C. Initial and boundary conditions

The computational domain has a size of

[0, L_{x}] \times [0, L_{y})

⁠, where the y direction is periodic. The aspect ratio of the domain is

L_{x} / L_{y} = 4

for the Ma = 0.15 and Ma = 0.30 cases and

L_{x} / L_{y} = 2

for the Ma = 0.45 case. The interface between the light and heavy fluid is located at

x / λ_{0} = 0

⁠. The flow is initialized with a hydrostatic isothermal background state at the initially uniform background temperature, T₀, that satisfies

\frac{\partial p^{H}}{\partial x} = - ρ^{H} g,

(21)

where p^H and ρ^H are the background hydrostatic pressure and density, respectively.

In order to well resolve the initial mixing width between the light and heavy fluids, the interface is initially diffused. The interface is smoothed with the error function by following a previous work³⁴ to set the mass fraction of the lighter fluid,

Y_{2}^{H}

⁠, using a characteristic initial thickness of the interface δ as

Y_{2}^{H} (x) = \frac{1}{2} [1 + erf (\frac{x}{δ})] .

(22)

Therefore, the background hydrostatic pressure, p^H, is given by

\begin{matrix} ln (\frac{p^{H} (x)}{p^{I}}) = \frac{g}{T_{0}} \int_{0}^{x} \frac{2}{(R_{2} - R_{1}) erf (ξ / δ) + (R_{1} + R_{2})} d ξ, \end{matrix}

(23)

where p^I is the pressure at the interface if the problem is initialized with a discontinuous interface. Similarly, the interface density for a discontinuous RTI is defined as

ρ^{I} = \frac{p^{I}}{R^{I} T_{0}},

(24)

where R^I is defined as

R^{I} = \frac{R_{g}}{(W_{1} + W_{2}) / 2} .

(25)

Narrow-band multi-mode perturbation with nine modes is added to the initial location of the interface between the fluids for the mass fraction fields. The perturbed displacement in the x direction, η, from the original interface location is given by

η (y) = \frac{η_{0}}{3} \sum_{m = 12}^{20} [cos (\frac{2 π m y}{L_{y}} + ϕ_{m})],

(26)

where m is the mode number and the middle mode number m_d is 16 in this work. The characteristic wavelength, λ₀, is defined as

λ_{0} = L_{y} / m_{d}

⁠.

ϕ_{m} \in [0, 2 π)

is the random phase shift of mode m to add randomness over different realizations. As a result, the perturbed mass fraction field of the lighter fluid is given by

Y_{2}^{H} (x, y) = \frac{1}{2} [1 + erf (\frac{x - η (y)}{δ})] .

(27)

The characteristic initial thickness of the interface δ is set to be 4% of λ₀ for all cases in this work.

Subsonic outflow boundary conditions are used at the left and right boundaries with the Navier–Stokes characteristic nonreflecting method.^34,38 The characteristic-based boundary method makes use of the local one-dimensional inviscid (LODI) relations to prevent any spurious inward-propagating acoustic waves. The target pressure values for the nonreflecting method are calculated from the hydrostatic equation at the domain boundaries.

III. DETAILS OF DIRECT NUMERICAL SIMULATIONS

A. Numerical methods

All multi-species simulations in this work are performed using Hydrodynamics Adaptive Mesh Refinement Simulator (HAMeRS).³⁹ It is a hydrodynamics code with Cartesian adaptive mesh refinement and local time stepping (or sub-cycling) capabilities. Conservative flux correction at coarse-fine grid interfaces⁴⁰ is implemented for any kind of conservative fluxes in the telescoping form. The parallelization of the code and the management of Cartesian grid cells are facilitated by the Structured Adaptive Mesh Refinement Application Infrastructure (SAMRAI) library^41–45 from Lawrence Livermore National Laboratory (LLNL). Feature-based sensors are implemented to identify regions of interest for adaptive mesh refinement (AMR). These include gradient and wavelet sensors⁴⁶ on different flow fields. The 13-point fourth-order dispersion-relation-preserving (DRP) scheme by Bogey and Bailly⁴⁷ re-cast in the conservative flux difference form^39,48 is used for the discretization of the convective flux. Explicit tenth-order finite difference schemes are used to discretize the diffusive and viscous terms as source terms of the equations. Navier–Stokes nonreflecting characteristic boundary conditions^34,38 with the ghost cell approach⁴⁹ are implemented for outflow domain boundaries. A third-order total variation diminishing Runge–Kutta (RK-TVD) scheme⁵⁰ is adopted for time integration with the local time stepping for each grid level.

B. Grid settings and realizations

Various cases with different combinations of background stratification strength characterized by the isothermal Mach number, Ma, and problem Reynolds number, Re₀, are designed for the DNS runs. Table I lists all of these cases for the 2 D multi-mode RTI problem. The index following the first letter, Ma: 015, 030, and 045, denotes the Mach numbers, 0.15, 0.30, and 0.45, respectively; and the index following the second letter, Re: 1, 2, and 3, represents the Reynolds number, Re₀: 3187.5, 6375.0, or 12 750.0, respectively. We should note that while the range of Re₀ chosen in this work is lower than that in the previous 2D single-mode compressible RTI paper by Wieland et al.¹⁵ (in the range of 25 500–102 000), these multi-mode DNS are computationally more expensive, as they are conducted in a much larger domain with each mode number being repeated around 12–20 times in the periodic dimension. Also note that if the domain width L_y instead of λ₀ is used in Eqs. (3) and (19), similar to some previous studies^17,19,20 on multi-mode RTI, both the problem isothermal Mach numbers and Reynolds numbers would be at higher values. For instance, the case Ma045Re3 would have an isothermal Mach number of 1.8 and a Reynolds number of 816 000. The isothermal Mach number based on L_y is also included in Table I for each case for reference.

TABLE I.

Parameters for the DNS cases. N_x and N_y are the numbers of grid cells in the x- and y-directions, respectively, for the base grid. For reference, the values in the parentheses in the second column represent values for the alternative isothermal Mach number, $M a^{*} = \sqrt{g L_{y}} / c_{0}$ ⁠, defined by the domain width in some previous studies^17,19,20 instead of the perturbation wavelength λ₀.

Cases	Ma	Re₀	$L_{x} / L_{y}$	Base grid resolution, $N_{x} \times N_{y}$	Number of mesh levels
Ma015Re1	0.15 (0.6)	$3187.5$	4	$1024 \times 256$	3
Ma030Re1	0.30 (1.2)	$3187.5$	4	$1024 \times 256$	3
Ma045Re1	0.45 (1.8)	$3187.5$	2	$1024 \times 512$	2
Ma015Re2	0.15 (0.6)	$6375.0$	4	$1024 \times 256$	3
Ma030Re2	0.30 (1.2)	$6375.0$	4	$1024 \times 256$	3
Ma045Re2	0.45 (1.8)	$6375.0$	2	$1024 \times 512$	2
Ma015Re3	0.15 (0.6)	12 750.0	4	$2048 \times 512$	3
Ma030Re3	0.30 (1.2)	12 750.0	4	$2048 \times 512$	3
Ma045Re3	0.45 (1.8)	12 750.0	2	2048 × 1024	2

Cases	Ma	Re₀	$L_{x} / L_{y}$	Base grid resolution, $N_{x} \times N_{y}$	Number of mesh levels
Ma015Re1	0.15 (0.6)	$3187.5$	4	$1024 \times 256$	3
Ma030Re1	0.30 (1.2)	$3187.5$	4	$1024 \times 256$	3
Ma045Re1	0.45 (1.8)	$3187.5$	2	$1024 \times 512$	2
Ma015Re2	0.15 (0.6)	$6375.0$	4	$1024 \times 256$	3
Ma030Re2	0.30 (1.2)	$6375.0$	4	$1024 \times 256$	3
Ma045Re2	0.45 (1.8)	$6375.0$	2	$1024 \times 512$	2
Ma015Re3	0.15 (0.6)	12 750.0	4	$2048 \times 512$	3
Ma030Re3	0.30 (1.2)	12 750.0	4	$2048 \times 512$	3
Ma045Re3	0.45 (1.8)	12 750.0	2	2048 × 1024	2

In addition, 16 realizations of simulations are performed for each combination case of Re₀ and Ma for better statistical convergence. Three levels of meshes, with two levels of adaptive mesh refinement, are used for the Ma = 0.15 and Ma = 0.30 cases, while the Ma = 0.45 case only has two levels of meshes. For the realization simulations with $R e_{0} = 3187.5$ and $R e_{0} = 6375.0$ ⁠, base grids of the Ma = 0.15 and Ma = 0.30 cases have 256 grid points in the y direction, while that of the Ma = 0.45 case has 512 grid points. Finer base grids are used for the highest Reynolds number simulations with $R e_{0} = 12 750.0$ to resolve smaller finest scales in the simulations. At this highest Re₀, the base grids of the Ma = 0.15 and Ma = 0.30 cases have 512 grid points, and that of the Ma = 0.45 case has 1024 grid points in the y direction. A constant refinement ratio of 1:2 is used across different grid levels in each direction. A finer base grid is used in Ma = 0.45 cases as the flow experiences larger pressure and density gradients even far away from the mixing layer width due to the stronger background stratification. A value sensor based on the mass fraction is used to identify mixing regions for AMR. Acoustic Courant–Friedrichs–Lewy (CFL) number of 0.5 is used for the Ma = 0.15 case, and 0.2 is used for the stronger background stratification cases. The diffusive CFL number is half of that of the acoustic CFL number for each case. The grid convergence analysis is also performed on one of the realizations for each combination case of Re₀ and Ma, and the results are compared with those obtained from simulations with base grids refined by a factor of 2 in each direction. Through grid sensitivity analysis, it is verified that all quantities of interest studied in the later sections are well grid converged for the corresponding realization of each case. While only one realization is used in the grid convergence analysis, all statistical results presented in this work are ensemble-averaged results over all realizations for each case.

IV. MIXING DYNAMICS

A. Flow visualization of mass fraction

In Fig. 1, the effects of stratification strength on the RTI are illustrated by comparing the visualizations of the mass fraction fields at three different non-dimensional time instants at

R e_{0} = 12 750.0

⁠, which is the largest problem Reynolds number. The non-dimensional or normalized time,

t^{*} = t / t_{r}

⁠, is introduced using the characteristic timescale defined as

t_{r} = \sqrt{λ_{0} / g},

(28)

to facilitate direct comparison of the cases with different background stratification strengths, or Ma. As A is kept constant in this study, we do not include A in the definition of t_r to account for the Atwood number effect, and the initial perturbation wavelength is used as the characteristic length scale in the definition of t_r. The first-row panel of Fig. 1 presents contours of the mass fraction of the heavy fluid, Y_h, at an early time,

t^{*} = 10

⁠. The mass fraction contours are essentially the same for all cases at this time. This indicates that the effects of the background stratification strength are not dominant at this early stage, at least for the mass fraction fields. The second-row panel of Fig. 1 presents Y_h at

t^{*} = 20

during the flow evolution. At this relatively late time, the differences between the contours of the Y_h start becoming noticeable among the three cases where there is a slow-down in the growth of the mixing layer for the moderate and high Ma cases (Ma030Re3 and Ma045Re3). At a later time,

t^{*} = 60

⁠, the low Ma number case (Ma015Re3) has the largest mixing width compared to the other two cases. The bubbles (the penetration of lighter fluid into heavier fluid) and spikes (the penetration of heavier fluid into lighter fluid) are still penetrating into the unmixed regions without getting considerably mixed at

t^{*} = 60

(third row of the Fig. 1). The figures also suggest an enhanced mixing within the mixing layer for the high Ma case (Ma045Re3) where the mass fraction variations, and also the density gradients, within the mixing layer are less sharp. Such reduction in the density gradients due to molecular mixing is expected to further prevent any re-growth for the high Ma case. An increase in the background stratification strength also leads to shorter and more homogeneously mixed mixing layers for the moderately and strongly stratified cases compared to the weakly stratified case. Such effects are most noticeable for the strongly stratified case. In addition, as the investigated Atwood number in this study is small, there is no noticeable difference between the mass, volume, and mole fraction visualizations in this study. A comparison of the spatial profiles based on mass, volume, and mole fractions can be found in Appendix A.

FIG. 1.

$Effects of the isothermal background stratification strength on the heavy fluid mass fraction, Yh, of the Ma = 0.15 (Ma015Re3), Ma = 0.30 (Ma030Re3), and Ma = 0.45 (Ma045Re3) cases at R e 0 = 12 750.0 from left to right at different times. Top row: t * = 10; middle row: t * = 20; and bottom row: t * = 60.$

View large Download slide

Effects of the isothermal background stratification strength on the heavy fluid mass fraction, Y_h, of the Ma = 0.15 (Ma015Re3), Ma = 0.30 (Ma030Re3), and Ma = 0.45 (Ma045Re3) cases at $R e_{0} = 12 750.0$ from left to right at different times. Top row: $t^{*} = 10$ ⁠; middle row: $t^{*} = 20$ ⁠; and bottom row: $t^{*} = 60$ ⁠.

We also present the Reynolds number effects on the heavy fluid mass fraction at three different problem Re₀ for the three different Ma cases at $t^{*} = 20$ and $t^{*} = 50$ in Figs. 2 and 3, respectively. It should be reminded that the problem Re₀ is mainly controlled by the initial interface kinematic shear viscosity, ν₀, of the flow. From the figures, it can be seen that an increase in Re₀ has a limited effect on the heights of the bubbles and spikes for all Ma cases, while a higher Re₀ leads to finer structures in the mass fraction field. The small Re₀ dependence of the mixing layer thickness and other large scales suggests that the slow-down in the growth of the mixing layer under strong background stratification, i.e., for the higher Ma cases, is at least valid in the range of Re₀ investigated in this study. This small Reynolds number dependence behavior of the large scales is also consistent with the incompressible and compressible single-mode RTI DNS results by Wei and Livescu⁵¹ and Wieland et al.,¹⁵ respectively, where the former and latter show insignificant Reynolds number dependence behavior starting from Re = 7500.0 and Re = 25 500.0, respectively.

FIG. 2.

$Reynolds number effects on the heavy fluid mass fraction, Yh, of the Ma = 0.15 (top row), Ma = 0.30 (middle row), and Ma = 0.45 (bottom row) cases at R e 0 = 3187.5 , R e 0 = 6375.0, and R e 0 = 12 750.0 from left to right at t * = 20.$

View large Download slide

Reynolds number effects on the heavy fluid mass fraction, Y_h, of the Ma = 0.15 (top row), Ma = 0.30 (middle row), and Ma = 0.45 (bottom row) cases at $R e_{0} = 3187.5, R e_{0} = 6375.0$ ⁠, and $R e_{0} = 12 750.0$ from left to right at $t^{*} = 20$ ⁠.

FIG. 3.

$Reynolds number effects on the heavy fluid mass fraction, Yh, of the Ma = 0.15 (top row), Ma = 0.30 (middle row), and Ma = 0.45 (bottom row) cases at R e 0 = 3187.5 , R e 0 = 6375.0, and R e 0 = 12 750.0 from left to right at t * = 50.$

View large Download slide

Reynolds number effects on the heavy fluid mass fraction, Y_h, of the Ma = 0.15 (top row), Ma = 0.30 (middle row), and Ma = 0.45 (bottom row) cases at $R e_{0} = 3187.5, R e_{0} = 6375.0$ ⁠, and $R e_{0} = 12 750.0$ from left to right at $t^{*} = 50$ ⁠.

B. Time evolution of mixing quantities

The time evolution of different mixing quantities is discussed in this subsection. To study the effects of the compressibility through background stratification on the large-scale mixing in the RTI flows, the growth of the mixing layer width is presented in Fig. 4(a) for different Ma cases at the largest problem Re₀. The mixing layer width, W, is defined using the mass fractions of the heavy and light fluids, Y_h and Y_l (⁠

Y_{l} = 1 - Y_{h}

⁠), as⁵²^,

W = 4 \int {\bar{Y}}_{h} {\bar{Y}}_{l} d x,

(29)

where

\bar{\cdot}

represents the averaging in all homogeneous directions, i.e., y direction in this work, over all realizations. When

t^{*} < 10

⁠, i.e., when the compressibility effect on W is not observable yet, the evolution of the mixing width collapses well with the chosen t_r for all three Ma cases. As shown in Fig. 4(a), the non-dimensional mixing layer width,

W / λ_{0}

⁠, continuously grows for the lowest Ma case (Ma015Re3) at least until the end of the simulation time, which is close to the incompressible RTI growth.⁹ However, we observe a slower growth than a self-similar quadratic growth described in Eq. . Such slower growth is an additional evidence of the suppression effects of the isothermal background stratification on compressible RTI, even at relatively low Ma numbers. RTI.⁹ For the medium Ma case (Ma030Re3), the stronger background stratification starts to show effects at

t^{*} \approx 15

as the growth of W slows down compared to the case Ma015Re3. For the largest Ma case (Ma045Re3), the background stratification kicks in at an even earlier time,

t^{*} \approx 10

⁠, to suppress the growth of W compared to the low and medium Ma cases at the same Re₀. For this case, the growth of the mixing layer saturates at

t^{*} \approx 30

⁠, where W remains stationary after that. A similar slow-down of the mixing layer was also observed for the 2D single-mode compressible RTI cases under strong background stratifications in previous works,¹⁵ which was performed at least until

t^{*} \approx 20

⁠. However, they did not observe a complete cessation of the RTI growth for their Ma = 0.3 and Ma = 0.6 cases. This study not only extends the slow-down observation to the multi-mode compressible RTI but also shows much less possibility of re-growth for the moderately and strongly stratified cases, as the DNS cases are performed up to a very late time,

t^{*} = 100

⁠, when the values of W well converge to constant values, and the mixing layers become molecularly well mixed, as discussed below.

FIG. 4.

View large Download slide

Time evolution of (a) the normalized mixing width, $W / λ_{0}$ ⁠, and (b) the normalized bubble and spike heights, $h_{b, s} / λ_{0}$ ⁠, at $R e_{0} = 12 750.0$ ⁠. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3). The lines for the spikes are marked with crosses, and those for the bubbles are marked with circles.

In addition to the integral mixing width, Fig. 4(b) shows the heights of the bubbles, h_b, and spikes, h_s, that are defined with the 99% and 1% thresholds of $\bar{Y_{h}}$ ⁠, respectively. Similar to W, both quantities are normalized by λ₀ in the figure. As expected, the behavior of $h_{b, s}$ is similar to W, although $h_{b, s}$ is more sensitive to the dynamics of the mixing layer edges, especially at late times (⁠ $t^{*} > 60)$ ⁠. Similar to W, $h_{b, s}$ also shows cessation in their growth rates for the moderately and strongly stratified cases (Ma030Re3 and Ma045Re3) after $t^{*} \approx 40$ and $t^{*} \approx 30$ ⁠, respectively. Both h_b and h_s have almost the same growth in time but in opposite directions for all cases, and this quasi-symmetry in growth is consistent with the incompressible single-mode RTI DNS results of Wei and Livescu⁵¹ at the same Atwood number before the reacceleration regime. Interestingly, Wieland et al.¹⁵ found that there is an asymmetry in the growth of the bubble and spike heights for the single-mode compressible RTI, where the spikes are consistently at larger heights than the bubbles for cases at different Ma.

The mixedness parameter, Θ, is broadly used to estimate the molecular mixing state of the mixing layer of the hydrodynamic instabilities.^25,52–54 Θ is defined as

Θ = \frac{\int \bar{Y_{h} Y_{l}} d x}{\int {\bar{Y}}_{h} {\bar{Y}}_{l} d x},

(30)

which has a value between zero and unity for the cases from the fully unmixed to fully mixed mixing layers, respectively. Θ can also be viewed as the ratio of the amount of fluids molecularly mixed to the amount of fluids entrained through convection. Its asymptotic value for the 2D incompressible RTI was reported⁵⁴ to be around 0.55, which indicates a balance between molecular mixing and fluid entrainment across the mixing layer at late times. For the low Ma case (Ma015Re3), which has the weakest background stratification, Θ tends asymptotically to this previously reported incompressible value, as shown in Fig. 5(a). However, for the larger Ma cases (Ma030Re3 and Ma045Re3), Θ approaches the value of unity, which indicates a state with well-mixed fluids in the mixing layer. The state of a well molecularly mixed layer caused by the large scalar dissipation represents a significant reduction of the hydrodynamically unstable density variations. This leads to a very low possibility of re-growth for the mixing layer, which is consistent with the observations on the saturation of the growth rates of W and

h_{b, s}

at late times for the moderately and strongly stratified cases, as presented in Fig. 4.

FIG. 5.

View large Download slide

Time evolution of (a) the mixedness, Θ, and (b) the normalized domain-integrated scalar dissipation rate, $χ^{*}$ ⁠, at $R e_{0} = 12 750.0$ ⁠. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

The normalized domain-integrated scalar dissipation rate,

χ^{*}

⁠, is defined as the scalar dissipation rate, χ, integrated across the mixing layer and normalized with the mixing width, W, and the characteristic timescale, t_r,

χ = D \nabla Y_{h} \cdot \nabla Y_{h} = - D \nabla Y_{h} \cdot \nabla Y_{l},

(31)

χ^{*} = \frac{\int \bar{χ} t_{r} d x}{W} .

(32)

χ is the dissipation rate of the scalar variance based on mass fractions and is dominated by the small scales of the flows. Figure 5(b) presents the time evolution of χ for all Ma cases under investigation at the largest Re₀. Initially, the normalized scalar dissipation rate is high for each case due to large mass fraction gradients within the mixing layer. The scalar dissipation rate decreases as the molecular diffusion process smooths out the sharp mass fraction variations at the instability interface. At this early time, the flow experiences a large-scale stirring from the hydrodynamic instability. Such stirring leads to a large-scale entrainment of pure fluids at the mixing layer. As the flow becomes very chaotic across the mixing layer, the scalar dissipation rate again reaches a peak value at around

t^{*} = 18

for each case. Although the normalized scalar dissipation rate has a similar evolution for all three different Ma cases being studied, it should be noted that the normalization factor, W, continuously grows for the lowest Ma case (Ma015Re3), and

χ^{*}

stays slightly high. This suggests that the domain-integrated scalar dissipation rate is always higher than others for this weakly stratified case, as there is always pure fluid penetration into the continuously growing mixing layer. In addition,

W / λ_{0}

⁠, Θ, and

χ^{*}

have smooth evolution that indicates a well statistical convergence for these parameters. We note that, as discussed above, h_b and h_s become highly sensitive to the mixing layer edge dynamics and may require more than 16 realizations to have smoother evolution at

t^{*} > 60

⁠.

It is also important to explore the generality of the above results by investigating the effects of the Reynolds number on the mixing statistical quantities. Figure 6 presents the Reynolds number effects on W and $h_{b, s}$ for the cases with the weakest (Ma015Re1, Ma015Re2, and Ma015Re3) and strongest (Ma045Re1, Ma045Re2, and Ma045Re3) background stratification strengths using different Re₀. As seen from the figure, Re₀ has a very limited effect on the growth of the mixing layer as the values of W approach similar heights for all three Re₀ no matter whether Ma = 0.15 or Ma = 0.45. This is consistent with the time evolution of the bubble and spikes as shown earlier in the visualizations of the mass fraction field. This lack of dependence of the large-scale characteristics, such as the mixing width and the bubble/spike heights (not shown in the figure), on the Reynolds number was also shown by Wieland et al.¹⁵ and Wei and Livescu⁵¹ for the single-mode compressible RTI and incompressible RTI, respectively. These findings indicate that the cease of the mixing layer growth due to strong background stratification is highly universal for the multi-mode compressible RTI and independent of the Re₀ at least for the investigated range of Re₀. While the other two mixing metrics, Θ, and $χ^{*}$ ⁠, show some sensitivities to Re₀, in particular at early times when the flow is not chaotic and still developing. For example, $χ^{*}$ gets a slightly larger value for the case Ma015Re1 compared to the case Ma015Re3 at $t^{*} \approx 18$ as the diffusion transport coefficient is larger for the lower Reynolds number cases with the same Ma. At later times of the flow evolution, both Θ and $χ^{*}$ collapse as the flow becomes fully developed. This further assures that the ceasing of the mixing and the reduction of the mass fraction variations across the mixing layer caused by strong background stratification strength are quite universal behavior as they do not show much Reynolds number dependence, at least for the range of Re₀ chosen in this study. It should be noted that $χ^{*}$ does not have significant Reynolds number dependence at late times. Although the gradients of the mass fraction are smaller for the lower Reynolds number cases, $χ^{*}$ is defined as a quantity that consists of the product of the squared mass fraction gradient magnitude and the mass diffusivity, D, which scales with ν₀ as Sc is kept constant and is larger for the low Reynolds number cases. As a result, at early times, such as $t^{*} \approx 18$ ⁠, while the flow is yet to be developed, $χ^{*}$ is slightly larger for the low Re₀ cases as these cases have larger mass diffusivity. At later times, $χ^{*}$ becomes independent of the Reynolds number as the flow has developed, and the larger mass fraction gradients compensate for the smaller mass diffusivity for the high Reynolds number cases. This is also consistent with the flow visualization of the heavy fluid mass fraction field, where a higher value of Re₀ leads to finer features and larger gradients in the mass fraction field.

FIG. 6.

View large Download slide

Time evolution of (a) the normalized mixing width, $W / λ_{0}$ ⁠, (b) the mixedness, Θ, and (c) the normalized domain-integrated scalar dissipation rate, $χ^{*}$ ⁠, of the Ma = 0.15 cases (top row) and the Ma = 0.45 cases (bottom row) with different Reynolds numbers. Red solid line: $R e_{0} = 3187.5$ ⁠; green dashed line: $R e_{0} = 6375.0$ ⁠; and blue dash-dotted line: $R e_{0} = 12 750.0$ ⁠.

C. Spatial profiles of mixing quantities

To further investigate the mixing behavior of the compressible RTI, we present the mean spatial profiles of the heavy fluid mass fraction, ${\bar{Y}}_{h}$ ⁠, in Fig. 7 for all Ma cases at the largest Re₀. At the very early time, $t^{*} = 5$ ⁠, the profiles of all cases have a monotonic transition from heavy fluid region to light fluid region. However, at $t^{*} = 10$ ⁠, the profiles of ${\bar{Y}}_{h}$ show non-monotonic shapes as the pure fluids exchange locations within the mixing layer due to fluid entrainment while the molecular mixing is limited. When $t^{*} \geq 20$ ⁠, the enhanced molecular mixing smooths out the ${\bar{Y}}_{h}$ profiles toward more monotonic shapes again for all cases. It is important to note that at these late times, increased strength of the background stratification leads to larger gradients of the ${\bar{Y}}_{h}$ profiles near the edges of the mixing layer. This is also noticeable in the ${\bar{Y}}_{h}$ visualizations in Fig. 1 for the higher Ma cases (Ma030Re3 and Ma045Re3), where there is a more distinct separation between the pure fluid regions and the mixing layers compared to the low Ma case (Ma015Re3).

FIG. 7.

$Effects of the isothermal background stratification strength on the mean spatial profiles of the heavy fluid mass fraction, Y ¯ h, at R e 0 = 12 750.0 at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).$

View large Download slide

Effects of the isothermal background stratification strength on the mean spatial profiles of the heavy fluid mass fraction, ${\bar{Y}}_{h}$ ⁠, at $R e_{0} = 12 750.0$ at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

Next, we present the spatial profiles of the heavy fluid mass fraction variance, $\bar{Y_{h}^{'} Y_{h}^{'}}$ ⁠, at different time instants in Fig. 8, where $(\cdot)'$ represents the Reynolds fluctuation of a given quantity. At $t^{*} = 5$ ⁠, there is a symmetric distribution with a single peak for each case, and at $t^{*} = 10$ ⁠, there are double peaks which are consistent with the results of the non-monotonic ${\bar{Y}}_{h}$ profiles. A more interesting behavior is observed at later times, i.e., when $t^{*} \geq 20$ ⁠, as the $\bar{Y_{h}^{'} Y_{h}^{'}}$ profiles become asymmetric even with our chosen very small Atwood number. Results of the incompressible numerical simulations^{9,21,29,30,36} and quasi-incompressible experiments^55–57 of the buoyancy-driven flows suggest highly symmetric flow evolution when A < 0.2. Recently, the DNS results by Prine et al.⁵⁸ showed that for single-mode compressible RTI at A = 0.25, an increase in the background stratification strength leads to a significant difference between the evolution of the locations of the spikes and bubbles, as the former grows much faster than the latter. In addition, they showed that an increase in compressibility through the background stratification strength (higher Ma) prevents the formation of vortical structures at the bubble side of the flows, while the vortical structures are more resistant to the effects of background stratification at the spike side. In this study, the Atwood number is kept very small (A = 0.04), so one would not expect any asymmetric behavior on the spatial profiles of the mass fraction variance. Our results of the low Ma case, Ma015Re3, are consistent with the earlier observations on incompressible or quasi-incompressible RTI flows, as there is no significant difference between the evolution of the h_b and h_s. However, the mixing state of the flow surprisingly exhibits differences between the heavier and lighter fluid regions when Ma is increased, where the mass fraction variance is larger in the heavier fluid regions compared to the lighter fluid regions. Figure 9 shows the normalized mean spatial profiles of the scalar dissipation rate, $\bar{χ} t_{r}$ ⁠, between the cases Ma015Re3, Ma030Re3, and Ma045Re3. The scalar dissipation rate profiles stay quasi-symmetric until the late times of the flow evolution as the large-scale asymmetry in $\bar{Y_{h}^{'} Y_{h}^{'}}$ has not been reflected in the small-scale dissipation quantity. At $t^{*} = 50$ ⁠, the spatial profile becomes asymmetric for the case Ma045Re3, which can be attributed to the asymmetric $\bar{Y_{h}^{'} Y_{h}^{'}}$ profile as the scalar dissipation rate would be enhanced in regions where the variance is larger.

FIG. 8.

$Effects of the isothermal background stratification strength on the spatial profiles of the heavy fluid mass fraction variance, Y h ′ Y h ′ ¯, at R e 0 = 12 750.0 at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).$

View large Download slide

Effects of the isothermal background stratification strength on the spatial profiles of the heavy fluid mass fraction variance, $\bar{Y_{h}^{'} Y_{h}^{'}}$ ⁠, at $R e_{0} = 12 750.0$ at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

FIG. 9.

Effects of the isothermal background stratification strength on the normalized mean spatial profiles of the scalar dissipation rate, χ ¯ t r, at R e 0 = 12 750.0 at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

View large Download slide

Effects of the isothermal background stratification strength on the normalized mean spatial profiles of the scalar dissipation rate, $\bar{χ} t_{r}$ ⁠, at $R e_{0} = 12 750.0$ at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

The normalized mean density profiles, $\bar{ρ} / ρ^{I}$ ⁠, are shown in Fig. 10 for all Ma cases at the highest Re₀. At $t^{*} = 0$ ⁠, the mean density profiles near the interface at $x / λ_{0} = 0$ are non-monotonic with inverse gradients for all cases because of the initial conditions of unstable stratifications. In addition, due to the exponential-like relation between the density field and the height of the domain for these isothermally stratified cases, the initial density stratification is not symmetric, and such asymmetry becomes dominant by an increase in the isothermal background stratification strength. As a result, the initial background stratification would also contribute to the observed asymmetry in $\bar{Y_{h}^{'} Y_{h}^{'}}$ (see Fig. 8) for the high Ma cases. For incompressible RTI under constant acceleration, the unstable stratification persists indefinitely during the simulation or experiment as the pure fluids have constant densities and the pure fluids continuously fill the growing mixing layer. However, for the isothermally stratified compressible RTI, the densities of the pure fluids are not uniform and are dependent on location in the domain, as derived in Eq. (21) and shown in the initial conditions of Fig. 10. As a result, as seen from Fig. 10(c), the density profile rapidly becomes monotonic without any inverse gradient for the strongest Ma case (Ma045Re3) after $t^{*} = 10$ because of the coupled effects of mixing and background stratification. From the time evolution of the mixing width, we can see that the mixing layer continues to grow beyond $t^{*} = 10$ until $t^{*} \approx 30$ for this case. There is a similar observation for the moderately stratified case Ma030Re3, where the mixing width temporarily still increases after the mean density profile becomes monotonic at around $t^{*} = 20$ ⁠. This may be due to the inertia of bubbles/spikes, the existence of locally unstable density gradients that are not represented by the mean density profile because of the 2D chaotic motions, and also possibly the strong expansion and compression motions.²⁰ It is also notable that even in our weakest background stratification case, Ma015Re3, the mean density profile eventually becomes monotonic without inverse gradient at late time $t^{*} = 50$ ⁠. This suggests that the mixing layer growth may also cease for this low Ma case as well if the simulation is run for a longer time, which is significantly different than incompressible RTI evolution.

FIG. 10.

Normalized mean spatial profiles of the density, ρ ¯ / ρ I, of cases with different Ma at R e 0 = 12 750.0 at different times. Cyan solid line: t * = 5; red dashed line: t * = 10; green dash-dotted line: t * = 20; and blue dotted line: t * = 50; thin gray dashed line: initial conditions at t * = 0.

View large Download slide

Normalized mean spatial profiles of the density, $\bar{ρ} / ρ^{I}$ ⁠, of cases with different Ma at $R e_{0} = 12 750.0$ at different times. Cyan solid line: $t^{*} = 5$ ⁠; red dashed line: $t^{*} = 10$ ⁠; green dash-dotted line: $t^{*} = 20$ ⁠; and blue dotted line: $t^{*} = 50$ ⁠; thin gray dashed line: initial conditions at $t^{*} = 0$ ⁠.

V. FLOW ENERGETICS

The RTI at the interface leads to a rapid and efficient conversion of potential energy into kinetic energy.^21,36 For the incompressible case, the continuous growth of the mixing layer of RTI is caused by this energy conversion. This section investigates the effects of background stratification strength on the multi-mode compressible RTI evolution further through flow energetics.

A. Time evolution and spatial profiles of turbulent kinetic energy

Figure 11(a) presents the normalized domain-integrated turbulent kinetic energy,

{TKE}^{*}

⁠, which is the turbulent kinetic energy,

TKE

⁠, integrated across the mixing layer and normalized with the mixing width, W, and the reference turbulent kinetic energy,

{TKE}_{r}

⁠.

{TKE}_{r}

is defined with the interface density, ρ^I, and the gravity wave speed,

\sqrt{g λ_{0}}

⁠. The normalized integrated turbulent kinetic energy and the related quantities are given as

TKE = \frac{1}{2} ρ u_{i}^{″} u_{i}^{″},

(33)

{TKE}_{r} = \frac{1}{2} ρ^{I} \sqrt{g λ_{0}} \sqrt{g λ_{0}},

(34)

{TKE}^{*} = \frac{\int \bar{TKE} d x}{W \cdot {TKE}_{r}},

(35)

where

u^{″} = u - \tilde{u}

represents the Favre fluctuations of the velocity field using the Favre mean (density-weighted) of the velocity,

\tilde{u} = \bar{ρ u} / \bar{ρ}

⁠. As seen in Fig. 11(a),

{TKE}^{*}

approaches a constant value for the low Ma case (Ma015Re3) that is similar to the incompressible RTI flows.⁹ However, for the cases with stronger background stratification (Ma030Re3 and Ma045Re3),

{TKE}^{*}

decays when

t^{*} > 15

⁠. This decay is consistent with the cessation of the growth of the RTI mixing width, which is also an indication of the end of conversion from potential energy to kinetic energy.

FIG. 11.

Time evolution of (a) the normalized domain-integrated turbulent kinetic energy, TKE *, and (b) the Reynolds normal stress anisotropy, ⟨ R ̃ 11 ⟩ / ⟨ R ̃ 22 ⟩, at R e 0 = 12750.0. Red solid line: Ma = 0.15 at the interface (Ma015Re3); and green dashed line: Ma = 0.30 (Ma030Re3); blue dash-dotted line: Ma = 0.45 (Ma045Re3).

View large Download slide

Time evolution of (a) the normalized domain-integrated turbulent kinetic energy, ${TKE}^{*}$ ⁠, and (b) the Reynolds normal stress anisotropy, $⟨ {\tilde{R}}_{11} ⟩ / ⟨ {\tilde{R}}_{22} ⟩$ ⁠, at $R e_{0} = 12750.0$ ⁠. Red solid line: Ma = 0.15 at the interface (Ma015Re3); and green dashed line: Ma = 0.30 (Ma030Re3); blue dash-dotted line: Ma = 0.45 (Ma045Re3).

In a buoyancy-driven flow, the anisotropy of the Reynolds normal stress is another important quantity to examine, as the flow has directionality due to the constant acceleration. Both incompressible numerical^{9,12,13,21,27} and quasi-incompressible experimental measurements^56,57 show that buoyancy-driven flows are highly anisotropic in general. The anisotropy of the 2D flows can be quantified with the use of the Reynolds stress tensor, which is defined as

{\tilde{R}}_{i j} = \frac{\bar{ρ u_{i}^{″} u_{j}^{″}}}{\bar{ρ}} .

(36)

The anisotropy of the Reynolds normal stress components can be measured from

\frac{{\tilde{R}}_{11} + ϵ}{{\tilde{R}}_{22} + ϵ},

(37)

where ϵ that is close to zero is added to the denominator to prevent division by zero. The ratio increases when a larger portion of

\bar{TKE}

is contributed from

\bar{ρ} {\tilde{R}}_{11}

⁠, while the ratio is zero when there is no contribution to

\bar{TKE}

from that component. ϵ is also added to the numerator such that the ratio is one not only when the Reynolds normal stress is isotropic but also when

\bar{TKE}

is zero, i.e., when the flows are laminar. For simplicity, ϵ is not shown in any expression of anisotropy in the following text. To measure the Reynolds normal stress anisotropy over time, the anisotropy averaged in the central part of the mixing layer,

⟨ {\tilde{R}}_{11} ⟩ / ⟨ {\tilde{R}}_{22} ⟩

⁠, can be used, where

⟨ \cdot ⟩

donates a spatial averaging within the interior chaotic region of the mixing layer that satisfies

4 {\bar{Y}}_{h} {\bar{Y}}_{l} \geq 0.9 .

(38)

At the early stage of the instability, it is shown in Fig. 11(b) that

⟨ {\tilde{R}}_{11} ⟩ / ⟨ {\tilde{R}}_{22} ⟩

attains its largest value where more than 90% of

\bar{TKE}

is contributed from

\bar{ρ} {\tilde{R}}_{11}

for all three Ma cases at the highest Re₀. This is expected as the potential energy is first converted to this Reynolds normal stress component before the energy is transferred to other components through the pressure–strain redistribution. Later, the Reynolds normal stress anisotropy in the central part of the mixing layer rapidly decays when the flow develops for all three cases. However, it stays larger than one for the low Ma case (Ma015Re3) as the growth of the mixing layer and TKE continues. On the other hand, our DNS results suggest that the flows become quasi-isotropic for the higher Ma cases with stronger background stratification strengths (Ma030Re3 and Ma045Re3), as the anisotropy metrics fluctuate around one. It is worth mentioning that the isotropization of the Reynolds normal stress is also rare in RMI,^25,59,60 where the post-shocked flows stay anisotropic in late times, even in the absence of any accelerations. We believe the rapid isotropization of the Reynolds normal stress components for the higher Ma cases is an outcome of the cease of the RTI growth caused by the compressible effects, where the pressure-strain redistribution term is more effective in equalizing TKE among the Reynolds normal stress components when the production term of TKE becomes small or is even no longer present.

Figure 12 presents the normalized mean spatial profiles of turbulent kinetic energy, $\bar{TKE} / {TKE}_{r}$ ⁠, at four distinct times for all three Ma cases at the highest Re₀. Although the compressibility effects through background stratification on the mixing layer growth are small at $t^{*} = 5$ ⁠, we can still observe that a stronger background stratification strength leads to a smaller peak of the spatial profile. Such reduction of the magnitude of $\bar{TKE} / {TKE}_{r}$ due to stronger background stratification strength becomes more and more observable at later times. In addition, Fig. 13 presents the spatial profiles of the Reynolds normal stress anisotropy at different time instants for all cases at the highest problem Re₀. Initially, the peaks of ${\tilde{R}}_{11} / {\tilde{R}}_{22}$ around $x / λ_{0} = 0$ are greater than 10 for all cases, which indicates that the Reynolds normal stress components are highly anisotropic around the middle portion of the mixing layer. Nevertheless, the Reynolds normal stress anisotropy of the flow slightly decreases to the value of roughly 1.5 at the core of the mixing layer at a later time, $t^{*} = 20$ ⁠. At a very late time, $t^{*} = 50$ ⁠, the spatial profiles of the Reynolds normal stress anisotropy fluctuates around unity for the higher Ma cases (Ma030Re3 and Ma045Re3), whereas the flow stays larger than one within a large portion of the mixing layer for the low Ma case (Ma015Re3). This shows that the flows are more isotropic for the higher Ma cases compared to that of the case Ma015Re3, which is also consistent with the time evolution of the anisotropy in Fig. 11(b).

FIG. 12.

Effects of the isothermal background stratification strength on the normalized mean spatial profiles of turbulent kinetic energy, TKE ¯ / TKE r, at R e 0 = 12 750.0 at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

View large Download slide

Effects of the isothermal background stratification strength on the normalized mean spatial profiles of turbulent kinetic energy, $\bar{TKE} / {TKE}_{r}$ ⁠, at $R e_{0} = 12 750.0$ at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

FIG. 13.

Effects of the isothermal background stratification strength on the spatial profiles of Reynolds normal stress anisotropy, R ̃ 11 / R ̃ 22, at R e 0 = 12 750.0 at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

View large Download slide

Effects of the isothermal background stratification strength on the spatial profiles of Reynolds normal stress anisotropy, ${\tilde{R}}_{11} / {\tilde{R}}_{22}$ ⁠, at $R e_{0} = 12 750.0$ at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

B. Transport of turbulent kinetic energy

Researchers have extensively studied the kinetic energy transport for turbulence induced by incompressible RTI and buoyancy-driven turbulence.^{21,26,36,53,61,62} On the other hand, the analysis of kinetic energy budgets for compressible RTI is still limited except some very recent works.^20,63,64 In this subsection, the budgets of the turbulent kinetic energy are analyzed for the multi-mode compressible RTI under different background stratification strengths.

The transport equation of the turbulent kinetic energy per unit mass,

k = {\tilde{R}}_{i i} / 2

⁠, of the studied problem has the following one-dimensional (1D) mean form after taking the average in the homogeneous y direction over all realizations:

\begin{matrix} \underset{term (I)}{\underset{︸}{\frac{\partial \bar{ρ} k}{\partial t}}} \underset{term (II)}{\underset{︸}{+ \frac{\partial (\bar{ρ} \tilde{u} k)}{\partial x}}} = \underset{term (III)}{\underset{︸}{a_{1} (\frac{\partial \bar{p}}{\partial x} - \frac{\partial {\bar{τ}}_{11}}{\partial x}) - \bar{ρ} {\tilde{R}}_{11} \frac{\partial \tilde{u}}{\partial x}}} \\ \underset{term (IV)}{\underset{︸}{- \frac{1}{2} \frac{\partial (\bar{ρ u_{i}^{″} u_{i}^{″} u^{″}})}{\partial x} - \frac{\partial (\bar{u^{'} p^{'}})}{\partial x} + \frac{\partial (\bar{u_{i}^{'} τ_{i 1}^{'}})}{\partial x}}} \\ \underset{term (V)}{\underset{︸}{+ \bar{p^{'} \frac{\partial u_{i}^{'}}{\partial x_{i}}}}} \underset{term (VI)}{\underset{︸}{- \bar{τ_{i j}^{'} \frac{\partial u_{i}^{'}}{\partial x_{j}}}}}, \end{matrix}

(39)

where the LHS consists of rate of change [term (I)] and convection [term (II)]. The RHS consists of production [term (III)], turbulent transport [term (IV)], pressure-dilatation [term (V)], and dissipation [term (VI)]. a is the turbulent mass flux velocity given by

a = \bar{ρ' u^{'}} / \bar{ρ}

⁠, where

u^{'}

and

ρ'

are the Reynolds fluctuations of the velocity and density fields, respectively. Also, note that

\bar{ρ} k = \bar{TKE}

⁠.

The transport equation of k is closely related to that of the Favre mean kinetic energy per unit mass, which is defined as

\tilde{K} = \tilde{u} \cdot \tilde{u} / 2

⁠. The averaged transport equation of

\tilde{K}

in the 1D mean form is given by

\begin{array}{l} \frac{\partial \bar{ρ} \tilde{K}}{\partial t} + \frac{\partial (\bar{ρ} \tilde{u} \tilde{K})}{\partial x} = \bar{ρ} \tilde{u} g - \frac{\partial \bar{ρ} \tilde{u} {\tilde{R}}_{11}}{\partial x} - \tilde{u} \frac{\partial \bar{p}}{\partial x} + \bar{ρ} {\tilde{R}}_{11} \frac{\partial \tilde{u}}{\partial x} + \tilde{u} \frac{{\bar{τ}}_{11}}{\partial x} . \end{array}

(40)

Unlike k and its transport equation,

\tilde{K}

and its transport equation are dependent on the chosen reference frame. The term

\bar{ρ} \tilde{u} g

denotes the production to convert the potential energy into the mean Favre mean kinetic energy. As

\tilde{u} = \bar{u} + a

⁠, the term can be split into two components

\bar{ρ} \tilde{u} g = \underset{term (I)}{\underset{︸}{\bar{ρ} \bar{u} g}} + \underset{term (II)}{\underset{︸}{\bar{ρ} a_{1} g}} .

(41)

While term (I) changes under the Galilean transformation, term (II) does not. Similarly,

\tilde{u} {\bar{p}}_{, 1}

can also be decomposed as

\tilde{u} \frac{\partial \bar{p}}{\partial x} = \bar{u} \frac{\partial \bar{p}}{\partial x} + a_{1} \frac{\partial \bar{p}}{\partial x} .

(42)

By inspecting the transport equations of

\tilde{K}

and k, it can be seen that

\bar{ρ} a_{1} g

is an agent that is invariant under the Galilean transformation of the reference frame to convert potential energy into Favre mean kinetic energy if both a₁ and g have the same signs (or a₁ is positive for the problem being studied). Moreover, the transfer of energy between

\tilde{K}

and k is only contributed by

a_{1} {\bar{p}}_{, 1}

and

\bar{ρ} {\tilde{R}}_{11} {\tilde{u}}_{, 1}

⁠. The turbulent mass flux velocity component in the streamwise direction, a₁, undoubtedly plays a critical role in the conversion of potential energy to mean and turbulent kinetic energies in the buoyancy-driven flows. Due to its critical role in the flow energetics of acceleration-driven multi-species flow instabilities, the turbulent mass flux velocity is studied in many previous works on incompressible RTI,^36,65 buoyancy-driven variable-density turbulence,^21,26 and RMI.^60,66,67

Figure 14 shows the time evolution of the normalized streamwise turbulent mass flux velocity component,

{⟨ a_{1} ⟩}^{*}

⁠, which is a₁ averaged in the middle part of the mixing layer and normalized by the gravity wave speed for all Ma cases at the highest Re₀

{⟨ a_{1} ⟩}^{*} = \frac{⟨ a_{1} ⟩}{\sqrt{g λ_{0}}} .

(43)

FIG. 14.

View large Download slide

Time evolution of the normalized mean of turbulent mass flux velocity in the middle part of the mixing layer, ${⟨ a_{1} ⟩}^{*}$ ⁠, at $R e_{0} = 12 750.0$ ⁠. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

As seen in the figure, ${⟨ a_{1} ⟩}^{*}$ is the largest (and is always positive) for the lowest Ma case (Ma015Re3) as the mixing layer continuously grows, which means the pure fluid motions always convert the potential energy to the kinetic energy. However, ${⟨ a_{1} ⟩}^{*}$ rapidly decays at $t^{*} \approx 10$ for the cases Ma030Re3 and Ma045Re3. More interestingly, ${⟨ a_{1} ⟩}^{*}$ becomes negative for these higher Ma cases [see Fig. 14(b)]. A negative value of the mean turbulent mass flux velocity indicates that there are fluid motions in the reserve directions compared to the earlier times in the middle portion of the domain where the heavier fluid moves up in the negative x direction, whereas lighter fluid migrates down in the opposite direction. This motion is consistent among different realizations and is evident after the ensemble averaging. This behavior is attributed to the result of the inertia of the flow when the stronger background stratification terminates the growth of the mixing layer, which leads to a rapid decay in ${⟨ a_{1} ⟩}^{*}$ ⁠. Such rapid inertial motion causes undershoot in the time evolution of ${⟨ a_{1} ⟩}^{*}$ to obtain a negative value.

The time evolution of the normalized domain-integrated production term, ${TKE}_{prod}^{*}$ ⁠, and dissipation term, ${TKE}_{diss}^{*}$ ⁠, of k for different Ma cases at the highest Re₀ is shown in Fig. 15. The production and dissipation terms are given by terms (III) and (VI), respectively, in Eq. (39). These two transport terms are normalized with ${TKE}_{r} \cdot W / t_{r}$ ⁠. From the figure, it can be seen that the time evolution of ${TKE}_{prod}^{*}$ is very similar to that of ${⟨ a_{1} ⟩}^{*}$ ⁠, as the production is essentially only contributed by the component $a_{1} {\bar{p}}_{, 1}$ ⁠, which is marked with circles in the plots. As a result, ${TKE}_{prod}^{*}$ is also negative and acts as an anti-production term, except at very early times for the higher Ma cases (Ma030Re3, Ma045Re3). On the other hand, ${TKE}_{diss}^{*}$ is always negative to dissipate $\bar{ρ} k$ ⁠.

FIG. 15.

View large Download slide

Time evolution of the normalized total production term ${TKE}_{prod}^{*}$ ⁠, the normalized production component associated with the turbulent mass flux velocity, and the normalized dissipation term, ${TKE}_{diss}^{*}$ of the turbulent kinetic energy $\bar{ρ} k$ at $R e_{0} = 12 750.0$ ⁠. The terms are given by Eq. (39). Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3). The lines for the normalized total production are not marked with any symbols, while the lines for the normalized dissipation are marked with crosses. The normalized production component associated with the turbulent mass flux velocity is represented by the circles only, and no lines.

The spatial profiles of the normalized budgets of the turbulent kinetic energy, $\bar{ρ} k$ (or $\bar{TKE}$ ⁠), are shown in Fig. 16 for the three Ma cases at the highest problem Re₀ at $t^{*} = 15$ and $t^{*} = 25$ ⁠. The budget terms given by Eq. (39) are normalized by ${TKE}_{r} / t_{r}$ ⁠. From the figure, it can be seen that the production, turbulent transport, and dissipation terms are the only relevant terms in the budgets for different cases at the chosen times, while other terms are negligible. At $t^{*} = 15$ ⁠, the production and the turbulent transport terms are dominant in the middle part of the mixing layer for the two lower Ma cases (Ma015Re3 and Ma030Re3). In this core region of the mixing layer, the turbulent transport terms are very negative. While the production term is still positive for the case Ma045Re3 at $t^{*} = 15$ in the middle part of the layer, it is small compared to the dissipation and the turbulent transport terms, and the overall rate of change of $\bar{ρ} k$ is very negative. At the edge of the mixing layer, the turbulent transport term is the only dominant term and is positive for each case. As a result, the spatial profiles of $\bar{ρ} k$ continue to spread at this time for all cases. At a later time, $t^{*} = 25$ ⁠, the budgets of the Ma015Re3 and Ma030Re3 cases are still quite similar compared to those at $t^{*} = 15$ ⁠, except that there are more fluctuations in the spatial profiles of the turbulent transport term. As the flows are more chaotic at this time, more realizations may be needed for smoother, or statistically better converged, profiles for this term. However, the general shape of this term at $t^{*} = 25$ can still be observed, which is quite similar to that at $t^{*} = 15$ for each case. Its effect is still to transfer $\bar{ρ} k$ from the core to the edges of the mixing layers for all cases, including case Ma045Re3. The production term of the cases Ma015Re3 and Ma030Re3 remains largely positive to continuously convert potential energy into kinetic energy, similar to the observation that this term is large over a long period of time in other studies of incompressible RTI and buoyancy-driven mixing.^21,26,36 Nevertheless, the production term has become negative across a large portion of the mixing layer for case Ma045Re3 due to negative turbulent mass flux velocity at this time. All relevant terms are negative at the core of the mixing layer, and hence the peak of $\bar{ρ} k$ (or $\bar{TKE}$ ⁠) decreases rapidly for this case around this time instant.

FIG. 16.

View large Download slide

Comparison of the normalized budget terms of turbulent kinetic energy, ρk, of Ma = 0.15 (Ma015Re3) (top row), Ma = 0.30 (Ma030Re3) (middle row), and Ma = 0.45 (Ma045Re3) (bottom row) cases at $R e_{0} = 12 750.0$ at two different times. The terms are given by Eq. (39). Cyan solid line: production [term (III)]; red dashed line: pressure-dilatation [term (V)]; green dash-dotted line: turbulent transport [term (IV)]; blue dash-dot-dotted line: dissipation [term (VI)]; orange dash-triple-dotted line: negative of convection due to streamwise velocity component associated with turbulent mass flux; thin black dotted line: summation of all terms (rate of change in the simulation frame).

C. Vortex dynamics

The study of vortex dynamics, particularly the evolution of the vorticity, can provide an additional dimension to examine the effects of background stratification on the multi-mode compressible RTI. Gauthier¹⁹ studied the vortex dynamics through helicity for the multi-mode compressible RTI, but only one background stratification strength was considered. While the effects of background stratification on the vorticity dynamics of compressible RTI were addressed by Wieland et al.,¹⁵ it was only limited to single-mode RTI. In this section, the effects of background stratification strength on the vorticity evolution for the multi-mode compressible RTI are studied. Figure 17 shows a qualitative comparison of the contours of the only nonzero component (z-component) of vorticity for the 2D flows, $ω = \nabla \times u$ ⁠, normalized by t_r for all three Ma cases at the largest Re₀, i.e., Ma015Re3, Ma030Re3, and Ma045Re3, at three different time instants: $t^{*} = 10, t^{*} = 20$ ⁠, and $t^{*} = 60$ ⁠. The first-row panel of the figure compares the contours between the three cases at $t^{*} = 10$ ⁠. At this time, the flows are still not very chaotic, as the roll-up of the bubbles and spikes is still very coherent for all cases. In addition, the shapes of the vortices are essentially the same for all cases, similar to the mass fraction contour plots in Fig. 1. However, it is also observed that an increase in the background stratification strength through Ma decreases the magnitude of the vorticity, as indicated by the minimum and maximum values of the colorbars. The reduced strength of vorticity with Ma is an early indication of background stratification effects on compressible RTI. At the later time, $t^{*} = 20$ ⁠, not only the strength but also the shapes of the vortices have become different between the three cases. The slow-down in the growth of the mixing layer due to the background stratification effect is reflected in the vortical structures where the structures are confined in a narrower region with higher Ma. At the much later time, $t^{*} = 60$ (shown by the third row of Fig. 17), the low Ma case (Ma015Re3) continues to show the largest vorticity magnitude within the mixing layer. However, the strength of the normalized vorticity is reduced compared to its strength at the earlier time, $t^{*} = 20$ ⁠. Although the flow is more chaotic at this time instant, the large-scale bubbles and spikes are still distinguishable and continue to penetrate into the unmixed regions. On the other hand, the flows have less distinguishable large-scale structures for the higher Ma cases (Ma030Re3 and Ma045Re3) due to larger levels of molecular mixing, and this is consistent with the observation of Y_h in Fig. 1. The effects of the background stratification strength on the evolution of the vorticity for these multi-mode compressible RTI cases are similar to those found by Wieland et al.¹⁵ for single-mode compressible RTI where the vorticity magnitude decreases and the propagation of the “fronts” of nonzero vorticity is suppressed with a larger Ma.

FIG. 17.

View large Download slide

Effects of the isothermal background stratification strength on the normalized vorticity, $ω / t_{r}$ ⁠, of the Ma = 0.15 (Ma015Re3), Ma = 0.30 (Ma030Re3), and Ma = 0.45 (Ma045Re3) cases from left to right at $R e_{0} = 12 750.0$ at different times. Top row: $t^{*} = 10$ ⁠; middle row: $t^{*} = 20$ ⁠; and bottom row: $t^{*} = 60$ ⁠.

The enstrophy, Ω, can be used to quantify the background stratification effects through Ma on the small-scale vortical structures. Ω is a measure of the vorticity strength within the mixing layer and is largely contributed by the small scales of the flows. It is defined as

Ω = ρ ω \cdot ω .

(44)

The normalized enstropy,

Ω^{*}

⁠, is defined as the normalization of the domain-integrated enstrophy by the mixing width, W, the interface density, ρ^I, and the characteristic timescale, t_r

Ω^{*} = \frac{t_{r}^{2} \int \bar{Ω} d x}{W ρ^{I}} .

(45)

Figure 19(a) presents the time evolution of

Ω^{*}

for the three different Ma cases at the highest Re₀. At early times,

t^{*} < 10

⁠, the time evolution of

Ω^{*}

collapses well for the three different cases. At later times, the strength of the vortical motions is the largest for the lowest Ma case (Ma015Re3), which is consistent with the vorticity contour plots. At very late times,

t^{*} > 40

⁠, while the

Ω^{*}

is still the largest for the case Ma015Re3, it decreases in this case, as well as other cases. Note that at this time period, the normalization factor W is still increasing for case Ma015Re3, while the growth of the mixing layer widths has already slowed down for case Ma030Re3 and ceased for case Ma045Re3. Figure 18 presents the spatial profiles of the enstrophy for the cases Ma015Re3, Ma030Re3, and Ma045Re3. We observe the strongest enstrophy at the center of the mixing layer for all three cases with quasi-symmetric shapes. The decrease in the enstrophy strength with stronger background stratification strength is consistent with Fig. 19(a).

FIG. 18.

Effects of the isothermal background stratification strength on the normalized mean spatial profiles of enstrophy, Ω ¯ / Ω r, at R e 0 = 12 750.0 at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

View large Download slide

Effects of the isothermal background stratification strength on the normalized mean spatial profiles of enstrophy, $\bar{Ω} / Ω_{r}$ ⁠, at $R e_{0} = 12 750.0$ at different times. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

FIG. 19.

View large Download slide

Time evolution of (a) the normalized domain-integrated enstrophy, $Ω^{*}$ ⁠, and (b) the normalized domain-integrated root mean square of baroclinic torque, ${\sqrt{\bar{B_{3}^{2}}}}^{*}$ ⁠, at $R e_{0} = 12 750.0$ ⁠. Red solid line: Ma = 0.15 (Ma015Re3); green dashed line: Ma = 0.30 (Ma030Re3); and blue dash-dotted line: Ma = 0.45 (Ma045Re3).

In addition to the enstrophy, the baroclinic torque is also studied as it is the production term of the vorticity,

ω

⁠. For 2D flows, this vector only has one non-zero component, which is given by

B_{3} = B_{z} = \frac{1}{ρ^{2}} \nabla ρ \times \nabla p .

(46)

The normalized domain-integrated root mean square (rms) of the baroclinic torque,

{\sqrt{\bar{B_{3}^{2}}}}^{*}

⁠, is given by

{\sqrt{\bar{B_{3}^{2}}}}^{*} = \frac{t_{r}^{2} \int \sqrt{\bar{B_{3}^{2}}} d x}{W} .

(47)

From Fig. 19(b), it can be seen that the time evolution of

{\sqrt{\bar{B_{3}^{2}}}}^{*}

follows closely that of the enstrophy, which is expected as this quantity is the source of vorticity. The strength of baroclinic torque decreases with stronger background stratification strength, which is similar to the single-mode compressible RTI findings.¹⁵ For the higher Ma cases Ma030Re3 and Ma045Re3,

{\sqrt{\bar{B_{3}^{2}}}}^{*}

asymptotically converges to a value around 0.6–0.8, indicating that there is still a consistent enstrophy production within the mixing layer until the end of the simulations. This asymptotic behavior is mainly attributed to the remaining highly unmixed fluids near the edges of the mixing layers, which are indicated by the double peaks in the spatial profiles of

\bar{Y_{h}^{'} Y_{h}^{'}}

in Fig. 8 for the two cases at late times.

Furthermore, the small-scale vortical motions of the flows are sensitive to the Reynolds number, and more details are discussed in Appendix B.

VI. CONCLUSION

In this study, we explored the background isothermal stratification effects on the 2D multi-mode compressible RTI mixing and energetics with fully compressible multi-species direct numerical simulations. The RTI is driven by the initial inverse density gradient at the perturbed two-fluid interface under constant acceleration. Although the background stratification strength has a small effect on the magnitude of the initial density jump at the interface compared with the Atwood number effect, it determines the background density gradient coupled with the initial jump. The background density gradient was found to largely affect the late-time mixing layer evolution using three cases under different strengths of background stratification. In this study, the strength of the stratification is controlled by the background isothermal Mach number, Ma, and they are chosen as 0.15, 0.30, and 0.45 for weakly, moderately, and strongly stratified cases, respectively. Similar to the 2D single-mode compressible isothermally stratified RTI,¹⁵ an increase in the background isothermal Mach number leads to a suppression of the growth of the RTI mixing layer. Simulations with durations that are substantially longer than those of the previous single-mode compressible RTI study show that the multi-mode RTI growth ceases even under relatively moderate stratification strength (Ma = 0.30). When the mixing layer grows, the non-monotonic behavior of the mean density gradient disappears rapidly for the strongly stratified cases, and the mixing layer growth slows down quickly. Eventually, after reaching a certain mixing layer height, the flow becomes well-mixed within the mixing layer as the pure fluid entrainment is ceased, and the high degree of molecular dissipation dominates the mixing process. In addition, the suppression of the mixing layer growth due to strong background stratification rapidly redistributes the energy within the two components of the Reynolds normal stress and leads to a quasi-isotropic flow.

The flow energetics of the multi-mode compressible RTI were also found to be significantly different than the incompressible RTI scenario. The TKE reaches a self-similar growing behavior for the incompressible RTI, whereas the TKE follows a non-linear evolution with rapid increase and decay stages due to the background stratification effects for the moderately and strongly stratified multi-mode compressible RTI cases. In contrast to the self-similar behavior of TKE in incompressible RTI,⁹ where the turbulent mass flux velocity is always in the same direction of acceleration, the turbulent mass flux velocity is in the opposite direction of the acceleration (i.e., leads to a negative value of the turbulent mass flux velocity if the acceleration is positive) at late times for these two cases, so the production term behaves as a TKE sink, rather than a source term as in incompressible RTI. The spatial profiles of turbulent mass flux velocity also suggest that the mass flux velocity plays a role in removing TKE at the core region of the mixing layer for the moderately and strongly stratified cases at late times when the mixing layer growth ceases. Finally, the Reynolds number dependence was investigated using three different problem Reynolds numbers for each isothermal Mach number case. Our results suggest that the findings on the effects of isothermal background stratification through Ma on most of the large-scale statistical quantities, such as the mixing width, bubble/spike heights, and TKE, appear to be universal as they do not show strong dependency on the Reynolds number. In addition, even the scalar dissipation rate was also found to be insensitive to the Reynolds number under proper scaling.

The findings of this study suggest a more complex flow evolution for the compressible RTI than the widely studied incompressible RTI. Future studies may investigate the effects of background stratification strength on the more realistic 3D compressible RTI and buoyancy-driven turbulence.

ACKNOWLEDGMENTS

D.A. would like to thank the U.S. National Science Foundation CBET Fluid Dynamics Program (Award No. 2234415) for the financial support that made this work possible. For the computations, this work used the Extreme Science and Engineering Discovery Environment (XSEDE),⁶⁸ which is supported by National Science Foundation Grant No. ACI-1548562. This work also used Stampede2 at the Texas Advanced Computing Center (TACC) through Allocation No. PHY200090 from the Advanced Cyber infrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation Grant Nos. 2138259, 2138286, 2138307, 2137603, and 2138296.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Denis Aslangil: Funding acquisition (lead); Investigation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Man Long Wong: Investigation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

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

APPENDIX A: COMPARISON OF THE SPATIAL PROFILES BASED ON MASS FRACTION AND VOLUME FRACTION

For flows with large density variations, the flow metrics based on mass, volume, or mole fractions could lead to different conclusions. In this study, the Atwood number is kept low, and we do not anticipate any sensitivity for a selection of these three options. The comparison of the mean spatial profiles of the heavy fluid mass fraction, ${\bar{Y}}_{h}$ ⁠, and heavy fluid volume fraction, ${\bar{Z}}_{h}$ ⁠, at different time instants is presented in Fig. 20. The heavy fluid volume fraction is given by $Z_{h} = ρ Y_{h} / ρ_{h}$ ⁠, where ρ_h is the species density of the heavy fluid. From the figure, it can be seen that the ${\bar{Y}}_{h}$ profiles are basically identical to the ${\bar{Z}}_{h}$ profiles. Figures 21 and 22 show that the variance and the normalized mean spatial profiles of the scalar dissipation rate based on mass and volume fractions are also essentially identical, respectively. We believe the spatial profiles based on mass fraction and volume fraction are basically the same because the Atwood number of the flow is low. Note that the volume fraction is equivalent to the mole fraction for species given by the ideal gas equation of state.

FIG. 20.

$Comparison of the mean spatial profiles of the heavy fluid mass fraction, Y ¯ h, and volume fraction, Z ¯ h, at different times at R e 0 = 12 750.0. The mass fraction mean profiles are represented by lines, and the volume fraction mean profiles are represented by symbols. Red solid line (circles): mass fraction (volume fraction) of Ma015Re3; and blue dash-dotted line (crosses): mass fraction (volume fraction) of Ma045Re3.$

View large Download slide

Comparison of the mean spatial profiles of the heavy fluid mass fraction, ${\bar{Y}}_{h}$ ⁠, and volume fraction, ${\bar{Z}}_{h}$ ⁠, at different times at $R e_{0} = 12 750.0$ ⁠. The mass fraction mean profiles are represented by lines, and the volume fraction mean profiles are represented by symbols. Red solid line (circles): mass fraction (volume fraction) of Ma015Re3; and blue dash-dotted line (crosses): mass fraction (volume fraction) of Ma045Re3.

FIG. 21.

$Comparison of the spatial profiles of the heavy fluid mass fraction variance, Y h ′ Y h ′ ¯, and volume fraction variance, Z h ′ Z h ′ ¯, at different times at R e 0 = 12 750.0. The mass fraction variance profiles are represented by lines, and the volume fraction variance profiles are represented by symbols. Red solid line (circles): mass fraction (volume fraction) of Ma015Re3; and blue dash-dotted line (crosses): mass fraction (volume fraction) of Ma045Re3.$

View large Download slide

Comparison of the spatial profiles of the heavy fluid mass fraction variance, $\bar{Y_{h}^{'} Y_{h}^{'}}$ ⁠, and volume fraction variance, $\bar{Z_{h}^{'} Z_{h}^{'}}$ ⁠, at different times at $R e_{0} = 12 750.0$ ⁠. The mass fraction variance profiles are represented by lines, and the volume fraction variance profiles are represented by symbols. Red solid line (circles): mass fraction (volume fraction) of Ma015Re3; and blue dash-dotted line (crosses): mass fraction (volume fraction) of Ma045Re3.

FIG. 22.

$Comparison of the normalized mean spatial profiles of the scalar dissipation rate, χ ¯ t r, based on mass fraction or volume fraction at different times at R e 0 = 12 750.0. The profiles based on mass fraction are represented by lines and the profiles based on volume fraction are represented by symbols. Red solid line (circles): mass fraction (volume fraction) of Ma015Re3; and blue dash-dotted line (crosses): mass fraction (volume fraction) of Ma045Re3.$

View large Download slide

Comparison of the normalized mean spatial profiles of the scalar dissipation rate, $\bar{χ} t_{r}$ ⁠, based on mass fraction or volume fraction at different times at $R e_{0} = 12 750.0$ ⁠. The profiles based on mass fraction are represented by lines and the profiles based on volume fraction are represented by symbols. Red solid line (circles): mass fraction (volume fraction) of Ma015Re3; and blue dash-dotted line (crosses): mass fraction (volume fraction) of Ma045Re3.

APPENDIX B: MORE ANALYSIS ON THE REYNOLDS NUMBER DEPENDENCE

In contrast to large scales such as bubble/spike heights, the vorticity has large sensitivity to the Reynolds number. In Figs. 23 and 24, the contours of the normalized vorticity component at different Re₀ are compared at two different time instants, $t^{*} = 20$ and $t^{*} = 50$ ⁠. The magnitude of the vorticity clearly increases with Re₀. It was shown that baroclinic torque and viscous diffusion are two major terms in the budgets of vorticity.¹⁵ As a result, as Re₀ increases, the viscous diffusion has a smaller effect due to reduced viscosity to offset the production effect of the baroclinic term, and thus, the vorticity is generated at a higher rate.

FIG. 23.

View large Download slide

Reynolds number effects on the normalized vorticity, $ω / t_{r}$ ⁠, of the Ma = 0.15 (top row), Ma = 0.30 (middle row), and Ma = 0.45 (bottom row) cases with $R e_{0} = 3187.5, R e_{0} = 6375.0$ ⁠, and $R e_{0} = 12 750.0$ from left to right at $t^{*} = 20$ ⁠.

FIG. 24.

View large Download slide

Reynolds number effects on the normalized vorticity, $ω / t_{r}$ ⁠, of the Ma = 0.15 (top row), Ma = 0.30 (middle row), and Ma = 0.45 (bottom row) cases with $R e_{0} = 3187.5, R e_{0} = 6375.0$ ⁠, and $R e_{0} = 12 750.0$ from left to right at $t^{*} = 50$ ⁠.

Figure 25 shows the time evolution of the Reynolds number dependence of the normalized turbulent kinetic energy (⁠ ${TKE}^{*}$ ⁠), enstrophy (⁠ $Ω^{*}$ ⁠), and baroclinic torque (⁠ ${\sqrt{\bar{B_{3}^{2}}}}^{*}$ ⁠). Similar to the mixing width or bubble/spike heights, ${TKE}^{*}$ which is also a large-scale quantity, shows limited sensitivity to Re₀. On the other hand, both $Ω^{*}$ and ${\sqrt{\bar{B_{3}^{2}}}}^{*}$ are mainly contributed by the smallest scales of the flow. Consistent with the large Reynolds number dependence observed for the vorticity contours, these vorticity-related small-scale quantities also show a large dependency on Re₀.

FIG. 25.

View large Download slide

Time evolution of (a) the normalized domain-integrated turbulent kinetic energy, ${TKE}^{*}$ ⁠, (b) the normalized domain-integrated enstrophy, $Ω^{*}$ ⁠, and (c) the normalized domain-integrated root mean square of baroclinic torque, ${\sqrt{\bar{B_{3}^{2}}}}^{*}$ ⁠, of the Ma = 0.15 cases (top row) and the Ma = 0.45 cases (bottom row) with different Reynolds numbers. Red solid line: $R e_{0} = 3187.5$ ⁠; green dashed line: $R e_{0} = 6375.0$ ⁠; and blue dash-dotted line: $R e_{0} = 12 750.0$ ⁠.

REFERENCES

1.

L.

Rayleigh

, “

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

Taylor

, “

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

,”

Proc. R. Soc. A

201

,

192

–

196

(

1950

).

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

Google Scholar

Crossref

3.

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

4.

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

5.

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

6.

A.

Banerjee

, “

Rayleigh–Taylor Instability: A status review of experimental designs and measurement diagnostics

,”

J. Fluids Eng.

142

,

120801

(

2020

).

https://doi.org/10.1115/1.4048349

Google Scholar

Crossref

7.

O.

Schilling

, “

Progress on understanding Rayleigh–Taylor flow and mixing using synergy between simulation, modeling, and experiment

,”

J. Fluids Eng.

142

,

120802

(

2020

).

https://doi.org/10.1115/1.4048518

Google Scholar

Crossref

8.

Y.

Zhou

,

R. J.

Williams

,

P.

Ramaprabhu

,

M.

Groom

,

B.

Thornber

,

A.

Hillier

,

W.

Mostert

,

B.

Rollin

,

S.

Balachandar

,

P. D.

Powell

et al, “

Rayleigh–Taylor and Richtmyer–Meshkov instabilities: A journey through scales

,”

Physica D

423

,

132838

(

2021

).

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

Google Scholar

Crossref

9.

J.

Ristorcelli

and

T.

Clark

, “

Rayleigh–Taylor turbulence: Self-similar analysis and direct numerical simulations

,”

J. Fluid Mech.

507

,

213

–

253

(

2004

).

https://doi.org/10.1017/S0022112004008286

Google Scholar

Crossref

10.

D.

Aslangil

,

A.

Banerjee

, and

A. G. W.

Lawrie

, “

Numerical investigation of initial condition effects on Rayleigh–Taylor instability with acceleration reversals

,”

Phys. Rev. E

94

,

053114

(

2016

).

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

Google Scholar

Crossref

PubMed

11.

G.

Boffetta

,

M.

Magnani

, and

S.

Musacchio

, “

Suppression of Rayleigh–Taylor turbulence by time-periodic acceleration

,”

Phys. Rev. E

99

,

033110

(

2019

).

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

Google Scholar

Crossref

PubMed

12.

D.

Aslangil

,

Z.

Farley

,

A. G. W.

Lawrie

, and

A.

Banerjee

, “

Rayleigh–Taylor instability with varying periods of zero acceleration

,”

J. Fluids Eng.

142

,

121103

(

2020

).

https://doi.org/10.1115/1.4048348

Google Scholar

Crossref

13.

D.

Aslangil

,

A. G. W.

Lawrie

, and

A.

Banerjee

, “

Effects of variable deceleration periods on Rayleigh–Taylor instability with acceleration reversals

,”

Phys. Rev. E

105

,

065103

(

2022

).

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

Google Scholar

Crossref

PubMed

14.

S.

Gauthier

and

B.

Le Creurer

, “

Compressibility effects in Rayleigh–Taylor instability-induced flows

,”

Philos. Trans. R. Soc., A

368

,

1681

–

1704

(

2010

).

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

Google Scholar

Crossref

15.

S. A.

Wieland

,

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

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.

T.

Luo

and

J.

Wang

, “

Effects of Atwood number and stratification parameter on compressible multi-mode Rayleigh–Taylor instability

,”

Phys. Fluids

33

,

115111

(

2021

).

https://doi.org/10.1063/5.0071437

Google Scholar

Crossref

18.

J. P.

Mellado

,

S.

Sarkar

, and

Y.

Zhou

, “

Large-eddy simulation of Rayleigh–Taylor turbulence with compressible miscible fluids

,”

Phys. Fluids

17

,

076101

(

2005

).

https://doi.org/10.1063/1.1965130

Google Scholar

Crossref

19.

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

20.

T.

Luo

and

J.

Wang

, “

Mixing and energy transfer in compressible Rayleigh–Taylor turbulence for initial isothermal stratification

,”

Phys. Rev. Fluids

7

,

104608

(

2022

).

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

Google Scholar

Crossref

21.

D.

Aslangil

,

D.

Livescu

, and

A.

Banerjee

, “

Effects of Atwood and Reynolds numbers 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

22.

R. D.

Richtmyer

, “

Taylor instability in shock acceleration of compressible fluids

,”

Commun. Pure Appl. Math.

13

,

297

–

319

(

1960

).

https://doi.org/10.1002/cpa.3160130207

Google Scholar

Crossref

23.

E. E.

Meshkov

, “

Instability of the interface of two gases accelerated by a shock wave

,”

Fluid Dyn.

4

,

101

–

104

(

1972

).

https://doi.org/10.1007/BF01015969

Google Scholar

Crossref

24.

B.

Thornber

,

J.

Griffond

,

O.

Poujade

,

N.

Attal

,

H.

Varshochi

,

P.

Bigdelou

,

P.

Ramaprabhu

,

B.

Olson

,

J.

Greenough

,

Y.

Zhou

et al, “

Late-time growth rate, mixing, and anisotropy in the multimode narrowband Richtmyer–Meshkov instability: The θ-group collaboration

,”

Phys. Fluids

29

,

105107

(

2017

).

https://doi.org/10.1063/1.4993464

Google Scholar

Crossref

25.

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

26.

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

27.

D.

Livescu

and

J. R.

Ristorcelli

, “

Variable-density mixing in buoyancy-driven turbulence

,”

J. Fluid Mech.

605

,

145

–

180

(

2008

).

https://doi.org/10.1017/S0022112008001481

Google Scholar

Crossref

28.

D.

Chung

and

D.

Pullin

, “

Direct numerical simulation and large-eddy simulation of stationary buoyancy-driven turbulence

,”

J. Fluid Mech.

643

,

279

–

308

(

2010

).

https://doi.org/10.1017/S0022112009992801

Google Scholar

Crossref

29.

D.

Aslangil

,

D.

Livescu

, and

A.

Banerjee

, “

Flow regimes in buoyancy-driven homogeneous variable-density turbulence

,” in

Progress in Turbulence VIII

, edited by

R.

Örlü

,

A.

Talamelli

,

J.

Peinke

, and

M.

Oberlack

(

Springer International Publishing

,

Cham

,

2019

), pp.

235

–

240

.

Google Scholar

Crossref

30.

D.

Aslangil

,

D.

Livescu

, and

A.

Banerjee

, “

Variable-density buoyancy-driven turbulence with asymmetric initial density distribution

,”

Physica D

406

,

132444

(

2020

).

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

Google Scholar

Crossref

31.

B. J.

Olson

and

A. W.

Cook

, “

Rayleigh–Taylor shock waves

,”

Phys. Fluids

19

,

128108

(

2007

).

https://doi.org/10.1063/1.2821907

Google Scholar

Crossref

32.

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

33.

M.-A.

Lafay

,

B.

Le Creurer

, and

S.

Gauthier

, “

Compressibility effects on the Rayleigh–Taylor instability between miscible fluids

,”

Europhys. Lett.

79

,

64002

(

2007

).

https://doi.org/10.1209/0295-5075/79/64002

Google Scholar

Crossref

34.

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

35.

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

36.

D.

Livescu

,

J.

Ristorcelli

,

R.

Gore

,

S.

Dean

,

W.

Cabot

, and

A.

Cook

, “

High-Reynolds number Rayleigh–Taylor turbulence

,”

J. Turbul.

10

,

N13

(

2009

).

https://doi.org/10.1080/14685240902870448

Google Scholar

Crossref

37.

F. A.

Williams

,

Combustion Theory

(

CRC Press

,

2018

).

Google Scholar

Crossref

38.

T. J.

Poinsot

and

S. K.

Lele

, “

Boundary conditions for direct simulations of compressible viscous flows

,”

J. Comput. Phys.

101

,

104

–

129

(

1992

).

https://doi.org/10.1016/0021-9991(92)90046-2

Google Scholar

Crossref

39.

M. L.

Wong

, “

High-order shock-capturing methods for study of shock-induced turbulent mixing with adaptive mesh refinement simulations

,” Ph.D. thesis (

Stanford University

,

2019

).

Google Scholar

40.

M. J.

Berger

and

P.

Colella

, “

Local adaptive mesh refinement for shock hydrodynamics

,”

J. Comput. Phys.

82

,

64

–

84

(

1989

).

https://doi.org/10.1016/0021-9991(89)90035-1

Google Scholar

Crossref

41.

B. T. N.

Gunney

and

R. W.

Anderson

, “

Advances in patch-based adaptive mesh refinement scalability

,”

J. Parallel Distrib. Comput.

89

,

65

–

84

(

2016

).

https://doi.org/10.1016/j.jpdc.2015.11.005

Google Scholar

Crossref

42.

B. T. N.

Gunney

,

A. M.

Wissink

, and

D. A.

Hysom

, “

Parallel clustering algorithms for structured AMR

,”

J. Parallel Distrib. Comput.

66

,

1419

–

1430

(

2006

).

https://doi.org/10.1016/j.jpdc.2006.03.011

Google Scholar

Crossref

43.

R. D.

Hornung

,

A. M.

Wissink

, and

S. R.

Kohn

, “

Managing complex data and geometry in parallel structured AMR applications

,”

Eng. Comput.

22

,

181

–

195

(

2006

).

https://doi.org/10.1007/s00366-006-0038-6

Google Scholar

Crossref

44.

R. D.

Hornung

and

S. R.

Kohn

, “

Managing application complexity in the SAMRAI object-oriented framework

,”

Concurrency Comput.

14

,

347

–

368

(

2002

).

https://doi.org/10.1002/cpe.652

Google Scholar

Crossref

45.

A. M.

Wissink

,

R. D.

Hornung

,

S. R.

Kohn

,

S. S.

Smith

, and

N.

Elliott

, “

Large scale parallel structured AMR calculations using the SAMRAI framework

,” in

Supercomputing, ACM/IEEE 2001 Conference

(

IEEE

,

2001

), pp.

22

–

22

.

Google Scholar

Crossref

46.

M. L.

Wong

and

S. K.

Lele

, “

Multiresolution feature detection in adaptive mesh refinement with high-order shock-and interface-capturing scheme

,” AIAA Paper No. AIAA

2016-3810

,

2016

.

Google Scholar

47.

C.

Bogey

and

C.

Bailly

, “

A family of low dispersive and low dissipative explicit schemes for flow and noise computations

,”

J. Comput. Phys.

194

,

194

–

214

(

2004

).

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

Google Scholar

Crossref

48.

A.

Subramaniam

,

M. L.

Wong

, and

S. K.

Lele

, “

A high-order weighted compact high resolution scheme with boundary closures for compressible turbulent flows with shocks

,”

J. Comput. Phys.

397

,

108822

(

2019

).

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

Google Scholar

Crossref

49.

E.

Motheau

,

A.

Almgren

, and

J. B.

Bell

, “

Navier–Stokes characteristic boundary conditions using ghost cells

,”

AIAA J.

55

,

3399

–

3408

(

2017

).

https://doi.org/10.2514/1.J055885

Google Scholar

Crossref

50.

C.-W.

Shu

and

S.

Osher

, “

Efficient implementation of essentially non-oscillatory shock-capturing schemes, II

,”

J. Comput. Phys.

83

,

32

–

78

(

1989

).

https://doi.org/10.1016/0021-9991(89)90222-2

Google Scholar

Crossref

51.

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

52.

D. L.

Youngs

, “

Numerical simulation of turbulent mixing by Rayleigh–Taylor instability

,”

Physica D

12

,

32

–

44

(

1984

).

https://doi.org/10.1016/0167-2789(84)90512-8

Google Scholar

Crossref

53.

W. H.

Cabot

and

A. W.

Cook

, “

Reynolds number effects on Rayleigh–Taylor instability with possible implications for type Ia supernovae

,”

Nat. Phys.

2

,

562

–

568

(

2006

).

https://doi.org/10.1038/nphys361

Google Scholar

Crossref

54.

W.

Cabot

, “

Comparison of two- and three-dimensional simulations of miscible Rayleigh–Taylor instability

,”

Phys. Fluids

18

,

045101

(

2006

).

https://doi.org/10.1063/1.2191856

Google Scholar

Crossref

55.

P.

Ramaprabhu

and

M. J.

Andrews

, “

Experimental investigation of Rayleigh–Taylor mixing at small Atwood numbers

,”

J. Fluid Mech.

502

,

233

–

271

(

2004

).

https://doi.org/10.1017/S0022112003007419

Google Scholar

Crossref

56.

A.

Banerjee

,

W. N.

Kraft

, and

M. J.

Andrews

, “

Detailed measurements of a Rayleigh–Taylor mixing layer from small to intermediate Atwood numbers

,”

J. Fluid Mech.

659

,

127

–

190

(

2010

).

https://doi.org/10.1017/S0022112010002351

Google Scholar

Crossref

57.

B.

Akula

and

D.

Ranjan

, “

Dynamics of buoyancy-driven flows at moderately high Atwood numbers

,”

J. Fluid Mech.

795

,

313

–

355

(

2016

).

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

Google Scholar

Crossref

58.

T.

Prine

,

D.

Aslangil

, and

M. L.

Wong

, “

Coupled effects of iso-thermal stratification strength and Atwood number on 2D single-mode compressible Rayleigh–Taylor instability

,” AIAA Paper No. AIAA 2023-1044,

2023

.

Google Scholar

59.

V. K.

Tritschler

,

B. J.

Olson

,

S. K.

Lele

,

S.

Hickel

,

X. Y.

Hu

, and

N. A.

Adams

, “

On the Richtmyer–Meshkov instability evolving from a deterministic multimode planar interface

,”

J. Fluid Mech.

755

,

429

–

462

(

2014

).

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

Google Scholar

Crossref

60.

M. L.

Wong

,

J. R.

Baltzer

,

D.

Livescu

, and

S. K.

Lele

, “

Analysis of second moments and their budgets for Richtmyer–Meshkov instability and variable-density turbulence induced by reshock

,”

Phys. Rev. Fluids

7

,

044602

(

2022

).

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

Google Scholar

Crossref

61.

A. W.

Cook

and

Y.

Zhou

, “

Energy transfer in Rayleigh–Taylor instability

,”

Phys. Rev. E

66

,

026312

(

2002

).

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

Google Scholar

Crossref

62.

Q.

Zhou

,

Y.-X.

Huang

,

Z.-M.

Lu

,

Y.-L.

Liu

, and

R.

Ni

, “

Scale-to-scale energy and enstrophy transport in two-dimensional Rayleigh–Taylor turbulence

,”

J. Fluid Mech.

786

,

294

–

308

(

2016

).

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

Google Scholar

Crossref

63.

Z.

Zhao

,

N.-S.

Liu

, and

X.-Y.

Lu

, “

Kinetic energy and enstrophy transfer in compressible Rayleigh–Taylor turbulence

,”

J. Fluid Mech.

904

,

A37

(

2020

).

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

Google Scholar

Crossref

64.

D.

Zhao

,

R.

Betti

, and

H.

Aluie

, “

Scale interactions and anisotropy in Rayleigh–Taylor turbulence

,”

J. Fluid Mech.

930

,

A29

(

2022

).

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

Google Scholar

Crossref

65.

D.

Livescu

,

J. R.

Ristorcelli

,

M. R.

Petersen

, and

R. A.

Gore

, “

New phenomena in variable-density Rayleigh–Taylor turbulence

,”

Phys. Scr.

2010

,

014015

.

https://doi.org/10.1088/0031-8949/2010/T142/014015

Crossref

66.

B. J.

Balakumar

,

G. C.

Orlicz

,

J. R.

Ristorcelli

,

S.

Balasubramanian

,

K. P.

Prestridge

, and

C. D.

Tomkins

, “

Turbulent mixing in a Richtmyer–Meshkov fluid layer after reshock: Velocity and density statistics

,”

J. Fluid Mech.

696

,

67

–

93

(

2012

).

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

Google Scholar

Crossref

67.

M.

Mohaghar

,

J.

Carter

,

B.

Musci

,

D.

Reilly

,

J.

McFarland

, and

D.

Ranjan

, “

Evaluation of turbulent mixing transition in a shock-driven variable-density flow

,”

J. Fluid Mech.

831

,

779

–

825

(

2017

).

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

Google Scholar

Crossref

68.

J.

Towns

,

T.

Cockerill

,

M.

Dahan

,

I.

Foster

,

K.

Gaither

,

A.

Grimshaw

,

V.

Hazlewood

,

S.

Lathrop

,

D.

Lifka

,

G. D.

Peterson

,

R.

Roskies

,

J.

Scott

, and

N.

Wilkins-Diehr

, “

XSEDE: Accelerating scientific discovery

,”

Comput. Sci. Eng.

16

,

62

–

74

(

2014

).

https://doi.org/10.1109/MCSE.2014.80

Google Scholar

Crossref

2023

Author(s)

Investigation of strong isothermal stratification effects on multi-mode compressible Rayleigh–Taylor instability

I. INTRODUCTION

II. PROBLEM DESCRIPTION

A. Governing equations

B. Non-dimensional problem parameters

C. Initial and boundary conditions

III. DETAILS OF DIRECT NUMERICAL SIMULATIONS

A. Numerical methods

B. Grid settings and realizations

IV. MIXING DYNAMICS

A. Flow visualization of mass fraction

B. Time evolution of mixing quantities

C. Spatial profiles of mixing quantities

V. FLOW ENERGETICS

A. Time evolution and spatial profiles of turbulent kinetic energy

B. Transport of turbulent kinetic energy

C. Vortex dynamics

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: COMPARISON OF THE SPATIAL PROFILES BASED ON MASS FRACTION AND VOLUME FRACTION

APPENDIX B: MORE ANALYSIS ON THE REYNOLDS NUMBER DEPENDENCE

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Investigation of strong isothermal stratification effects on multi-mode compressible Rayleigh–Taylor instability Free

I. INTRODUCTION

II. PROBLEM DESCRIPTION

A. Governing equations

B. Non-dimensional problem parameters

C. Initial and boundary conditions

III. DETAILS OF DIRECT NUMERICAL SIMULATIONS

A. Numerical methods

B. Grid settings and realizations

IV. MIXING DYNAMICS

A. Flow visualization of mass fraction

B. Time evolution of mixing quantities

C. Spatial profiles of mixing quantities

V. FLOW ENERGETICS

A. Time evolution and spatial profiles of turbulent kinetic energy

B. Transport of turbulent kinetic energy

C. Vortex dynamics

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: COMPARISON OF THE SPATIAL PROFILES BASED ON MASS FRACTION AND VOLUME FRACTION

APPENDIX B: MORE ANALYSIS ON THE REYNOLDS NUMBER DEPENDENCE

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Investigation of strong isothermal stratification effects on multi-mode compressible Rayleigh–Taylor instability