Effects of Isothermal Stratification Strength on Vorticity Dynamics for Single-Mode Compressible Rayleigh-Taylor Instability

TL;DR

This study investigates how varying isothermal stratification strength influences vorticity dynamics and growth behavior in compressible Rayleigh-Taylor instability through high-resolution simulations, revealing stratification-dependent growth suppression and asymmetry effects.

Contribution

It provides detailed analysis of stratification effects on RTI vorticity dynamics using wavelet-based simulations, highlighting the role of baroclinicity in growth suppression and asymmetry.

Findings

Stronger stratification suppresses RTI growth.

Weak stratification causes bubble-spike asymmetry.

Reynolds number has limited effect on bubble and spike heights.

Abstract

The effects of isothermal stratification strength on vorticity dynamics for single-mode Rayleigh-Taylor instability (RTI) are examined using two dimensional fully compressible wavelet-based direct numerical simulations. The simulations model low Atwood number ($A=0.04$) RTI development for four different stratification strengths, corresponding to Mach numbers from 0.3 (weakly stratified) to 1.2 (strongly stratified), and for three different perturbation Reynolds numbers, from 5,000 to 20,000. All simulations use adaptive wavelet-based mesh refinement to achieve very fine spatial resolutions at relatively low computational cost. For all stratifications, the bubble and spike go through the exponential growth regime, followed by a slowing of the RTI evolution. For the weakest stratification, this slow-down is then followed by a re-acceleration, while for stronger stratifications the suppression of RTI growth continues. Bubble and spike asymmetries are observed for weak stratifications, with bubble and spike growth rates becoming increasingly similar as the stratification strength increases. For the range of cases studied, there is relatively little effect of Reynolds number on bubble and spike heights, although the formation of secondary vortices becomes more pronounced as Reynolds number increases. The underlying dynamics are analyzed in detail through an examination of the vorticity transport equation, revealing that incompressible baroclinicity drives RTI growth for small and moderate stratifications, but increasingly leads to the suppression of vorticity production and RTI growth for stronger stratifications. These variations in baroclinicity are used to explain the suppression of RTI growth for strong stratifications, as well as the anomalous asymmetry in bubble and spike growth rates for weak stratifications.

Figures

Equations

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u_j)}{\partial x_j}$$

$$\frac{\partial (\rho u_i)}{\partial t}+\frac{\partial (\rho u_i{u}_j)}{\partial x_j}$$

$$\frac{\partial (\rho e)}{\partial t}+\frac{\partial (\rho e{u}_j)}{\partial x_j}$$

$$\frac{\partial (\rho Y_i)}{\partial t}+\frac{\partial (\rho Y_i{u}_j)}{\partial x_j}$$

$$p$$

$$e$$

$$R=Y_i{R}_i=\mathcal{R}\frac{Y_i}{W_i}.$$

$$c_{\mathrm{v}}=(c_{\mathrm{v}}{)}_i{Y}_i,$$

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)=2\mu S_{ij}^{\prime },$$

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

$$q_j=-k\frac{\partial T}{\partial x_j},$$

$$s_{ji}=\rho D\frac{\partial Y_i}{\partial x_j},$$

$$\rho (x_1,x_2,0)$$

$$p(x_1,x_2,0)$$

$${\rho }_{0\alpha }(x_1)$$

$$p_{0\alpha }(x_1)$$

$$R{e}_p\equiv \sqrt{\frac{Ag{\lambda }^3}{(1+A){\nu }^2}}\Rightarrow \nu =\sqrt{\frac{Ag{\lambda }^3}{(1+A)R{e}_p^2}},$$

$$A\equiv \frac{W_2-W_1}{W_2+W_1}.$$

$$M=\sqrt{\frac{{\rho }_{I}g\lambda }{p_I}}\Rightarrow M^2=\frac{g\lambda }{R_I{T}_0}.$$

$${\rho }_{0\alpha }^{\ast }(x_1^{\ast })$$

$$p_{0\alpha }^{\ast }(x_1^{\ast })$$

$$\frac{R_I}{R_{\alpha }}=\frac{2W_{\alpha }}{W_1+W_2}=1+(-1{)}^{\alpha }A\text{ for }\alpha =1,2.$$

$${\rho }_0^{\ast }(x_1^{\ast })$$

$$p_0^{\ast }(x_1^{\ast })$$

$$f_{\ge }(x)=\sum _k{c}_k^0{\varphi }_k^0(x)+{\sum _{j=0}^{\infty }\sum _{\alpha =1}^{2^n-1}\sum _l}_{|d_l^{\alpha ,j}|\ge \epsilon ||f||}{d}_l^{\alpha ,j}{\psi }_l^{\alpha ,j}(x),$$

$$\overline{\overline{f}}(t)=\frac{1}{{\lambda }^2}{\int }_0^{\lambda /2}{\int }_{-\lambda }^{\lambda }f(x_1,x_2,t)d{x}_{1}d{x}_2.$$

$$\frac{D{\omega }_i}{Dt}={\omega }_j{S}_{ij}-{\omega }_i{S}_{kk}-{ϵ}_{ijk}\frac{\partial v}{\partial x_j}\frac{\partial p}{\partial x_k}+{ϵ}_{ijk}\frac{\partial }{\partial x_j}\left[v\frac{\partial (2\mu S_{kl}^{\prime })}{\partial x_l}\right],$$

$$\frac{D{\omega }_i}{Dt}={\omega }_j{S}_{ij}+\nu \frac{{\partial }^2{\omega }_i}{\partial x_j\partial x_j}-{\omega }_i{S}_{kk}-{ϵ}_{ijk}\frac{\partial v}{\partial x_j}\frac{\partial p}{\partial x_k}-2\nu {ϵ}_{ijk}\frac{\partial }{\partial x_j}\left(\frac{S_{kl}^{\prime }}{v}\frac{\partial v}{\partial x_l}\right).$$

$$\frac{D{\omega }_i^{\ast }}{D{t}^{\ast }}={\omega }_j^{\ast }{S}_{ij}^{\ast }+\sqrt{\frac{A}{1+A}}\frac{1}{R{e}_p}\frac{{\partial }^2{\omega }_i^{\ast }}{\partial x_j^{\ast }\partial x_j^{\ast }}-{\omega }_i^{\ast }{S}_{kk}^{\ast }-{ϵ}_{ijk}\frac{\partial v^{\ast }}{\partial x_j^{\ast }}\frac{\partial p^{\ast }}{\partial x_k^{\ast }}-\sqrt{\frac{A}{1+A}}\frac{2}{R{e}_p}{ϵ}_{ijk}\frac{\partial }{\partial x_j^{\ast }}\left[S_{kl}^{\prime \ast }\frac{\partial (\ln v^{\ast })}{\partial x_l^{\ast }}\right].$$

$$v^{\prime \ast }({\text{x}}^{\ast },t^{\ast })\equiv v^{\ast }({\text{x}}^{\ast },t^{\ast })-\left[v_M^{\ast }(x_1^{\ast })-1\right],p^{\prime \ast }({\text{x}}^{\ast },t^{\ast })\equiv p^{\ast }({\text{x}}^{\ast },t^{\ast })-\left[p_M^{\ast }(x_1^{\ast })+x_1^{\ast }-\frac{1}{M^2}\right],$$

$$\frac{\partial v^{\prime \ast }}{\partial x_i^{\ast }}=\frac{\partial v^{\ast }}{\partial x_i}-M^2{v}_M^{\ast }{\delta }_{i1}$$

$$\frac{\partial p^{\prime \ast }}{\partial x_i^{\ast }}=\frac{\partial p^{\ast }}{\partial x_i}+(M^2{p}_M^{\ast }-1){\delta }_{i1}$$

$$-{ϵ}_{ijk}\frac{\partial v^{\ast }}{\partial x_j^{\ast }}\frac{\partial p^{\ast }}{\partial x_k^{\ast }}=-{ϵ}_{ijk}\frac{\partial v^{\prime \ast }}{\partial x_j^{\ast }}\frac{\partial p^{\prime \ast }}{\partial x_k^{\ast }}+\left(M^2{p}_M^{\ast }-1\right){ϵ}_{ij1}\frac{\partial v^{\prime \ast }}{\partial x_j^{\ast }}+\left(M^2{v}_M^{\ast }\right){ϵ}_{ij1}\frac{\partial p^{\prime \ast }}{\partial x_j^{\ast }},$$

$$\frac{D{\omega }_i^{\ast }}{D{t}^{\ast }}=$$

$$+\left(M^2{p}_M^{\ast }-1\right){ϵ}_{ij1}\frac{\partial v^{\prime \ast }}{\partial x_j^{\ast }}+\left(M^2{v}_M^{\ast }\right){ϵ}_{ij1}\frac{\partial p^{\prime \ast }}{\partial x_j^{\ast }}-\sqrt{\frac{A}{1+A}}\frac{2}{R{e}_p}{ϵ}_{ijk}\frac{\partial }{\partial x_j^{\ast }}\left[S_{kl}^{\prime \ast }\frac{\partial (\ln v^{\ast })}{\partial x_l^{\ast }}\right],$$

Full text

I Introduction

Rayleigh-Taylor instability (RTI) is formed at the interface of two fluids with different densities when an accelerative force is applied across the interface in the direction of the less dense fluid Rayleigh1884 ; Taylor1950 . Such a scenario arises in a number of engineering and physics problems, including inertial confinement fusion (ICF) Betti1998 ; Craxton2015 , supernova ignition fronts Hachisu1991 ; Zingale2005 ; Cabot2006 ; JordanIV2008 ; Swisher2015 , X-ray bursts Lewin1997 , and various topics in geophysics Sharp1984 ; Houseman1997 ; Lawrie2011 , to name just a few examples. In each of these cases, both flow and fluid compressibility may affect RTI growth; the former is related to the thermodynamic state and the stratification of background density and pressure fields, while the latter is related to the equation of state and differences in the specific heat ratio between the two fluids Livescu2004 (Ref. Gauthier2017 refers to these two effects as “static” and “dynamic” compressiblity, respectively). Moreover, in many of these cases, as well as in many experiments Dalziel1999 ; Ramaprabhu2004 ; Andrews2010 ; Lawrie2011 ; Zhou2017 , the Atwood number (namely, the density ratio of the two fluids normalized to take values between 0 and 1) is small and the Reynolds number is large, resulting in compressible dynamics driven by relatively small density differences over a wide range of length and time scales. In the present study, two dimensional (2D) fully compressible wavelet-based direct numerical simulations (DNS) are used to examine, from a dynamical standpoint, the evolution of low Atwood number RTI for different isothermal stratification strengths (i.e., flow compressibility) and Reynolds numbers.

Until relatively recently, many DNS studies of RTI used incompressible, low-Mach number, Boussinesq, or anelastic approximations to reduce the computational cost, often yielding valuable physical insights (see, e.g., Cook2001 ; Livescu2004 ; Livescu2008 ; Livescu2009 ; Schneider2016 ; Schneider2016b ). However, the present study is one of a growing number of fully compressible DNS analyses of RTI growth and characteristics. Lafay et al. Lafay2007 examined RTI growth in the linear regime for different compressibility strengths (addressing both flow and fluid compressibility), and Gauthier Gauthier2013 examined RTI growth into the nonlinear regime for two different stratification strengths. More recently, Reckinger et al. Reckinger2016 examined single-mode, 2D RTI growth rates for a range of stratification strengths, and Gauthier Gauthier2017 performed a comprehensive study of the dynamics of multi-mode, three-dimensional (3D) RTI for a relatively strongly stratified case. Both of these more recent studies employed variable-resolution numerical methods to achieve high Reynolds numbers within the context of fully compressible DNS; Reckinger et al. Reckinger2016 used the parallel adaptive wavelet collocation method (PAWCM) Reckinger2010 and Gauthier Gauthier2017 used an auto-adaptive multi-domain Chebyshev-Fourier method LeCreurer2008 . Using currently available computational resources, these and other adaptive techniques are unavoidable when performing fully compressible DNS at high Reynolds numbers. The present study correspondingly employs PAWCM to study the effects of flow compressibility and Reynolds number, including high Reynolds numbers, on RTI growth and dynamics.

Based in large part on observations from these prior fully compressible DNS studies, the general effects of compressibility on RTI growth are now relatively well understood. Despite some initial ambiguity regarding the specific impacts of compressibility (dating back, at least, to the studies by Bernstein & Book Bernstein1983 and Baker Baker1983 ), Livescu Livescu2004 used a linear analysis of the Navier-Stokes equations to show that, for isothermal background stratification, flow compressibility is associated with a reduction in the rate of RTI growth, while fluid compressibility is associated with an increase in the growth rate, compared to the corresponding incompressible case. A number of studies have confirmed these results, particularly with respect to the suppression of RTI growth by flow compressibility Lafay2007 ; Yu2008 ; Gauthier2010 ; Wei2012 ; Reckinger2016 ; Gauthier2017 . In particular, as the stratification strength of the background density field increases, an increasing suppression of RTI growth has been observed. Reckinger et al. Reckinger2016 further found that there are asymmetries in the locations and speeds of upward propagating low density fluid (i.e., “bubbles”) and downward propagating high density fluid (i.e., “spikes”), even at relatively small Atwood number, that may be different than in the incompressible limit. It was also shown by Reckinger et al. Reckinger2016 that drag and potential flow models are unable to predict the suppression of RTI growth for strong stratifications.

Compared to the effects of compressibility, Reynolds number effects on RTI growth have received somewhat less attention. Wei & Livescu Wei2012 used the incompressible variable-density form of the Navier-Stokes equations to show that, at early times, RTI growth rates are larger for smaller Reynolds numbers, due to diffusive effects, but that at long times, RTI growth rates are greater for larger Reynolds numbers. For sufficiently high Reynolds numbers, Wei & Livescu Wei2012 found that, at low Atwood numbers, bubble and spike growth rates become quadratic, similar to the multi-mode case and contrary to the ubiquitous “terminal velocity” assumption in previous studies. The single-mode growth rate ceases to depend on the Reynolds number at sufficiently large Reynolds numbers and may constitute an upper limit for the multi-mode growth rate. Using fully compressible DNS, Gauthier Gauthier2017 similarly found that smaller Reynolds numbers are associated with faster early growth rates of the turbulent mixing layer produced by the RTI. At later times, growth rates for higher Reynolds numbers are similar to, or exceed, those of lower Reynolds numbers. These results were, however, obtained for a single stratification strength, and it remains to be seen how these Reynolds number effects depend on stratification strength, if at all.

In order to understand compressibility and Reynolds number effects in more detail, several authors have examined the dynamics of the vorticity during RTI evolution, generally finding that changes in the baroclinic torque are responsible for changes in RTI growth rates. Lafay et al. Lafay2007 examined the linear regime and found that vorticity production decreases as the stratification strength increases. More recently, Schneider and Gauthier Schneider2016a performed a systematic study of vorticity during RTI growth using 3D multimode simulations that employ the Boussinesq approximation. This study showed that there is an increase in the strength of baroclinic torque production with time, although the contribution to the overall dynamics is dwarfed by the effects of nonlinear vortex stretching. Gauthier Gauthier2013 ; Gauthier2017 was the first to examine vorticity using fully compressible DNS, and showed the importance of baroclinic torque in producing vorticity during RTI growth for a single strongly stratified case. However, changes to the relative magnitudes of the various terms in the vorticity transport equation for different stratification strengths are still not completely understood in the fully compressible case.

Despite the improved understanding of compressibility and Reynolds number effects provided by the recent, primarily computational, studies noted above, a number of outstanding questions remain, and the present study is specifically focused on addressing the following: (i) How does the behavior of low Atwood number RTI depend on both stratification strength and Reynolds number?; (ii) What are the dynamical causes of the observed RTI phenomena?; and (iii) How do the dynamics (specifically, the vorticity dynamics) depend on stratification strength? The first question is motivated by the studies of Lafay et al. Lafay2007 and Wei & Livescu Wei2012 ; the former studied the effects of compressibility, but within the linear regime and for only one Reynolds number, while the latter studied a range of Reynolds numbers, but using an incompressible variable-density formulation of the Navier-Stokes equations that precluded the study of compressibility effects. The second and third questions are motivated primarily by the studies of Reckinger et al. Reckinger2016 , Schneider & Gauthier Schneider2016a , and Gauthier Gauthier2017 . The first of these studies observed bubble-spike asymmetries but not their dynamical causes, while the second and third studies both performed extensive analyses of the vorticity dynamics, but using the Boussinesq approximation (i.e., not a fully compressible study) and for only one stratification strength, respectively.

In the present paper, DNS are performed at low Atwood numbers (0.04 here, as compared to 0.1-0.7 in Reckinger2016 ) for different Reynolds numbers and different strengths of initial hydrostatic stratification, corresponding to Mach numbers between 0.3 (weak stratification) and 1.2 (strong stratification) Livescu2013 ; Livescu2004 ; Gerashchenko2016 . The DNS are performed using adaptive mesh refinement based on PAWCM, as described, validated, and implemented for RTI by Reckinger et al. Reckinger2010 ; Reckinger2016 . This method allows high spatial resolution to be used where it is needed (e.g., where density and velocity gradients are large), while reducing the total number of computational collocation points. The present focus on low Atwood numbers is motivated primarily by the observations of quadratic high Reynolds number single-mode RTI growth in regimes with similarly low Atwood number in the study by Wei & Livescu Wei2012 . In order to understand the dynamical causes of the observed results, the various terms in the vorticity transport equation are examined as functions of time and stratification strength.

It should be noted that several simplifications are made here to allow the underlying physics to be more easily understood. In particular, complex interactions of multiple wavelengths are eliminated by applying only single-mode initial perturbations to the unstable interface between the two fluids with differing densities. Moreover, in the classical incompressible case, where the density of both fluids is constant, RTI growth eventually leads to a re-acceleration of the bubble and spike tips, finally resulting in chaotic dynamics and development. In the compressible case, however, matters can become more complicated due to spatial and temporal variations in the background density, pressure, and temperature fields. The effects of changing any of these fields are largely unknown, and thus only isothermal initial stratifications are studied here to eliminate thermal effects, since the initial state is already in thermal equilibrium. Future work will explore the effects of multi-modal perturbations and different stratification types. Finally, the present simulations and analysis are performed in 2D in order to enable the examination of several different stratification strengths and Reynolds numbers. Each such simulation is computationally expensive and performing a similarly expansive study in 3D remains the focus of future research, due primarily to the need for substantially more extensive computational resources. The primary disadvantage of the present 2D approach is the resulting lack of nonlinear vortex stretching in the vorticity dynamics, although the absence of this effect does have the benefit of more clearly revealing the effects of baroclinicity on the dynamics.

The rest of this paper is organized as follows. The next section discusses the problem setup, including the governing equations and initialization of the RTI. Section III provides a brief discussion of how the wavelet-based adaptive method (i.e., PAWCM) was used to complete the simulations. In Section IV, the paper goes in depth into the results of this study, looking at the effects of stratification strength and Reynolds number on RTI growth. In Section V, the dynamics of the vorticity for fully compressible RTI are outlined and examined. Finally, a summary and conclusions are presented in Section VI.

II Description of the Physical Problem

In the present study, RTI occurs through the initial placement of a heavier fluid, denoted by index ‘2’ with molar mass $W_{2}$, above a lighter fluid, denoted by index ‘1’ with molar mass $W_{1}$, in the presence of a gravitational accelerative force. The addition of a perturbation causes the pseudo-stable initial condition to be lost, and the heavier fluid begins to fall into the lighter fluid in a spike-like formation, while the lighter fluid rises into the heavier fluid in a bubble-like formation. In the following, the fully compressible fluid flow equations solved by the DNS are outlined, followed by a description of the initial isothermal hydrostatic stratifications of different strengths (as characterized by a static Mach number).

II.1 Governing Equations

The numerical simulations solve the fully compressible Navier-Stokes equations for two miscible fluids given by Livescu2013

[TABLE]

where $\rho$ is the density, $u_{i}$ is the velocity in the $x_{i}$ direction, $p$ is the pressure, $g_{i}$ is the gravitational acceleration, $\tau_{ij}$ is the viscous stress tensor, $e$ is the specific total energy, $q_{i}$ is the heat flux, $T$ is the temperature, $(c_{\mathrm{p}})_{l}$ is the specific heat capacity at constant pressure for fluid $l$, $s_{ji}$ is the mass flux for fluid $i$ in the $x_{j}$ direction, and $Y_{i}$ is the mass fraction for the $i^{\textrm{th}}$ fluid. Note that, for a two-fluid system, $Y_{2}=1-Y_{1}$, and so Eq. (4) is only solved in the present simulations for $i=2$ (i.e., the heavier fluid). The pressure and caloric ideal gas laws are assumed to hold, so that the pressure and specific total energy can be expressed as

[TABLE]

where $R$ is the mixture gas constant defined in terms of the universal gas constant $\mathcal{R}$ and the molar mass of each fluid, $W_{i}$, as

[TABLE]

In the above expression, the species gas constant is defined as $R_{i}\equiv\mathcal{R}/W_{i}$. The mixture specific heat at constant volume, $c_{\mathrm{v}}$, appearing in Eq. (6) is similarly defined as

[TABLE]

where the specific heats at constant pressure and volume are related by $(c_{\mathrm{p}})_{i}=(c_{\mathrm{v}})_{i}+R_{i}$ and their mixture values by $c_{\mathrm{p}}=c_{\mathrm{v}}+R$. The specific heats at constant volume are assumed constant and the same for the two fluids, so that the mixture specific heat at constant pressure varies with the flow due to the different molar masses of the two fluids.

The viscous stress $\tau_{ij}$ in Eqs. (2) and (3) is assumed to be Newtonian and is given by

[TABLE]

where $S^{\prime}_{ij}=S_{ij}-S_{kk}\delta_{ij}/3$ is the deviatoric strain rate and the dynamic viscosity is given by $\mu=\rho\nu$, with the kinematic viscosity $\nu$ assumed to be constant (i.e., temperature independent and the same for both fluids) such that spatial and temporal variations in $\mu$ are due entirely to variations in $\rho$. The strain rate tensor $S_{ij}$ is given by

[TABLE]

The heat flux in Eq. (3) is written as

[TABLE]

where $k$ is the thermal conductivity, and the species mass flux in Eqs. (3) and (4) is defined as

[TABLE]

where $D$ is the mass diffusivity. For the range of parameters considered here, the baro-diffusion term is small in the mass flux, and Soret and Dufour effects are neglected in the mass and heat fluxes, respectively. Both $k$ and $D$, like the kinematic viscosity $\nu$, are assumed to be constant and temperature-independent, and both Prandtl and Schmidt numbers are unity.

The majority of fluid properties are taken to be the same between the two fluids for simplicity. This includes the kinematic viscosity, $\nu$, the heat conduction coefficient, $k$, and the mass diffusion coefficient, $D$. It should be noted that effects due to bulk viscosity and non-equilibrium thermodynamics are neglected in the simulations. Investigating these effects is beyond the scope of the present study, although Sagert et al. Sagert2015 and Lai et al. Lai2016 have recently made progress in this direction.

The system of equations given by Eqs. (1)-(12) is solved using the PAWCM numerical approach, which is described in Section III, for an RTI with a physical setup as outlined in the following section.

II.2 Initialization of Rayleigh Taylor Instability

The RTI problem is initialized in the DNS by imposing a perturbation on a stratified isothermal background state that is in hydrostatic equilibrium. The gravitational acceleration is assumed to be in the negative $x_{1}$ direction, such that $g_{i}=-g\delta_{i1}$, where $g$ is the magnitude of the gravitational acceleration. The resulting density, $\rho(x_{1},x_{2},t)$, and pressure, $p(x_{1},x_{2},t)$, fields at $t=0$ can be expressed as

[TABLE]

where $\rho_{0}$ and $p_{0}$ are hydrostatic initial background states and $\rho^{\prime}(x_{1},x_{2},0)$ and $p^{\prime}(x_{1},x_{2},0)$ represent the initial perturbations to the background states.

Assuming an isothermal background state at temperature $T_{0}$, the background density and pressure fields for fluid $\alpha$ (where $\alpha=[1,2]$ and summation over Greek indices is not implied) are given by

[TABLE]

where the initial interface between the two fluids lies at $x_{1}=0$, $p_{I}$ is the interfacial pressure, and $R_{\alpha}=\mathcal{R}/W_{\alpha}$ is the gas constant based on the molar mass of fluid $\alpha$. The heavier fluid ($\alpha=2$) is initially located above the interface for $x_{1}>0$ and the lighter fluid ($\alpha=1$) is initially located below the interface for $x_{1}<0$. A corresponding interfacial density is given using the ideal gas law as $\rho_{I}=p_{I}/(R_{I}T_{0})$ where $R_{I}=\mathcal{R}[(W_{1}+W_{2})/2]^{-1}$.

In each of the cases examined here, the kinematic viscosity $\nu=\mu/\rho$, which is constant and the same in both fluids, is set using the perturbation Reynolds number, $Re_{p}$, as

[TABLE]

where $\lambda$ is the wavelength of the applied perturbation and $A$ is the non-dimensional Atwood number defined as

[TABLE]

Note that in the present study, $W_{2}>W_{1}$ in order to generate RTI.

The degree of flow compressibility defined by the thermodynamic conditions enters the RTI problem by affecting both the background stratification and the further development of dilatational (non-zero velocity divergence) effects Livescu2004 ; Livescu2013 . While dilatational effects and their acoustic manifestations are usually characterized by the Mach number denoting the ratio between velocity and sound speed, together with dilatational to solenoidal kinetic energy ratios, stratification strength can also be recast as a Mach number. This can be done by re-expressing $gx_{1}/(R_{\alpha}T_{0})$ in Eqs. (15) and (16), as described below.

In the present study, the relevant incompressible limit is found by simultaneously increasing the background pressure and temperature to cause an increase in the speed of sound such that the density remains unaffected. This incompressible limit is also easily obtained in practice by uniformly heating a fixed volume of fluid. This results in the definition of an isothermal Mach number based on the ratio of the gravity wave speed, $\sqrt{g\lambda}$, and the isothermal speed of sound, $a_{0}=\sqrt{{p_{I}}/{\rho_{I}}}$ Livescu2004 ; Yu2008 . The resulting Mach number, $M$, is then given by

[TABLE]

It should be noted that $M$ is equivalent to stratification strength parameters used in prior studies of flow (or “static”) compressibility Schneider2016 ; Gauthier2017 ; Gauthier2017a , and that larger values of $M$ indicate stronger stratification.

Normalizing $\rho_{0\alpha}$ in Eq. (15) and $p_{0\alpha}$ in Eq. (16) by $\rho_{I}$, $g$, and $\lambda$, the non-dimensional background states can be rewritten as

[TABLE]

where the characteristic pressure is given as $\rho_{I}g\lambda$ and $x^{*}_{1}\equiv x_{1}/\lambda$ is a normalized distance variable. It can be shown that the ratio $R_{I}/R_{\alpha}$ can be written in terms of the Atwood number $A$ as

[TABLE]

Since $\alpha=1$ corresponds to the lighter fluid for which $x^{*}_{1}<0$ initially and $\alpha=2$ corresponds to the heavier fluid for which $x^{*}_{1}>0$, the non-dimensional background states $\rho^{*}_{0}$ and $p^{*}_{0}$ can be written in final form as

[TABLE]

where $\rho^{*}_{0}=\rho_{0}/\rho_{I}$, $p^{*}_{0}=p_{0}/(\rho_{I}g\lambda)$, with $(1-A)$ for $x^{*}_{1}<0$ (i.e., the lighter fluid) and $(1+A)$ for $x^{*}_{1}>0$ (i.e., the heavier fluid). The resulting initial background stratifications are shown for a variety of Mach numbers in Figure 1, where the size of the density difference at $x^{*}_{1}=0$ is determined by the value of $A$ ($A=0.04$ in the present study).

Following the procedure extensively outlined by Reckinger et al. Reckinger2016 , a single-mode velocity perturbation was applied at $t=0$ to initialize the RTI. Although they are not perfect representations of multi-mode engineering problems found in ICF and other practical applications, the present single-mode simulations can nevertheless be used to gain insights into compressibility-driven physics and dynamics. As shown by Reckinger et al. Reckinger2016 , single-mode simulations can expose any numerical directional bias in the code, which is generally hidden in multi-mode simulations. As a result, single-mode simulations allow the opportunity to ensure that the simulations are completely resolved from the initial state through to late times, and also allow simple checks for symmetry and the introduction of extraneous perturbation modes throughout the simulation. In addition, the results of Wei & Livescu Wei2012 show that single-mode RTI may represent the upper bound for the multi-mode growth rate at low Atwood numbers, when the Reynolds number is sufficiently large.

III Details of the Direct Numerical Simulations

Due to the spatial localization of the developing region, the RTI problem lends itself naturally to state-of-the-art adaptive grid numerical methods. In particular, to effectively capture the instability evolution, very long domains are needed to ensure that late-time growth is captured, but very small grid spacing is required to fully resolve the high gradients at the interface of the instability. For a static computational grid with fixed cell size, this results in a very dense grid and incredibly high computational costs. During the majority of the simulation, however, very fine grid resolutions far away from the interface are unnecessary and, as a result, high grid compression ratios can be achieved through the use of adaptive grid approaches. A method that has proven effective at achieving high compression ratios is the Parallel Adaptive Wavelet Collocation Method (PAWCM) Reckinger2010 ; Reckinger2016 , which is the method that is applied here.

III.1 Wavelet-Based Grid Adaptation

The PAWCM numerical approach has been applied previously to the simulation of compressible RTI by Reckinger et al. Reckinger2016 , where validation and details of the numerical method are exhaustively outlined. These details are repeated only briefly here, and the reader is referred to Reckinger2016 for additional information.

Fundamentally, PAWCM uses the natural properties of the wavelet transform to locate areas of steep gradients and to provide direct control over the grid cell size used to resolve the gradients. Essentially, through PAWCM, a flow field variable is transformed into wavelet space, resulting in wavelet basis functions and coefficients that are localized in both wave and physical spaces. From there, the coefficients are passed through a thresholding filter where all of the coefficients with magnitudes above the parameter $\varepsilon$ are kept, and any of those below $\varepsilon$ are set to zero. The resulting thresholded decomposition can thus be written for a generic variable $f$ as

[TABLE]

where $\phi_{k}$ are scaling functions on the coarsest level, $c_{k}$ are the corresponding coarse-level wavelet coefficients, $\psi_{l}$ are the scaling interpolating functions on any arbitrary level, $d_{l}$ are the coefficients to which the thresholding is applied, $l$ and $k$ represent physical grid points, and $\alpha$ and $j$ represent the wavelet family and level of resolution, respectively Vasilyev2000a ; Nejadmalayeri2015 . The effect of setting any one of the coefficients $d_{l}$ to zero is the removal of a grid point at that level of resolution. These coefficients take on large values for large gradients, and small values in relatively uniform regions. The effective resolution is set by a base grid size and the limit put on $j$ (referred to as $j_{max}$ herein). This results in the error being $\mathcal{O}(\varepsilon)$ and the resolution in a single direction being $p\cdot 2^{(j_{max}-1)}$, where $p$ is the base resolution Vasilyev2000a ; Schneider2010a ; Nejadmalayeri2015 .

As outlined in Reckinger et al. Reckinger2016 , PAWCM has been implemented in a way that enables it to work with finite difference approaches to solving governing equations such as those in Eqs. (1)-(4). In solving these equations, fourth-order central differences have been applied spatially, and a third-order total variation diminishing explicit Runge-Kutta scheme has been applied in time. The PAWCM algorithm is highly parallelized, having successfully run on up to 5,000 cores, and is able to perform arbitrary domain decompositions using the Zoltan library. It has a tree-like data structure for easy MPI communications, as well as direct error control. As a result, the additional computational overhead introduced by the wavelet methodology is offset by the capability to use many processors and to achieve grid compression ratios greater than 90% Vasilyev2000a ; Schneider2010a ; Nejadmalayeri2015 .

Substantial discussion was provided in Reckinger et al. Reckinger2016 regarding the flow variables on which to adapt the grid in the DNS. Since the wavelet method is so flexible, it is possible to adapt the grid on any flow field variable that is calculable and of interest. In the present study, adaptation for the initial time steps was performed using the vorticity, the norm of the strain rate tensor, and the gradient of the species mass fraction $Y_{2}$, in addition to the velocity and mass fraction fields. This approach allowed the RTI to develop with sufficient accuracy prior to further refining the grid on more complex flow variables at later times to reflect the increasing complexity of the flow. In particular, at late times in the present study, adaptation was performed using the baroclinic torque to ensure that this dynamically important term was fully resolved for the analysis of the vorticity dynamics. Additional details on grid convergence and resolution can be found in Reckinger2016 .

III.2 Simulation Setup

In the present study, PAWCM is used to solve the governing equations outlined in Section II.1 for the background and initial conditions described in Section II.2. The simulations have been carried out in 2D and the total domain size was $16\lambda$ in the $x_{1}$ direction and $\lambda$ in the $x_{2}$ direction, where $\lambda$ corresponds to the wavelength of the applied perturbation. The maximum effective grid resolution resulting from the adaptive wavelet approach was $\Delta x^{*}=2.4\times 10^{-4}$, where $\Delta x^{*}=\Delta x/\lambda$ and $\Delta x$ is the grid cell size.

The Atwood number studied was 0.04, and the Mach numbers used were 0.3 (nearly incompressible), 0.6, 0.9, and 1.2. The perturbation Reynolds numbers, $Re_{p}$, investigated were 5,000, 10,000, and 20,000, giving a total of twelve simulations performed in the present study (i.e., four different values of $M$, and three values of $Re_{p}$ for each $M$). The highest Reynolds number is of particular interest because it has been shown to be the minimum perturbation Reynolds number necessary to reach the chaotic growth regime for the incompressible limit of this particular case Wei2012 . Each of the simulations were performed up to a non-dimensional time of $t^{*}=t/\sqrt{\lambda/g}=20$, corresponding to the time at which the bubble and spike had reached heights of roughly $\lambda$ (or $x_{1}^{*}=\pm 1$) for the $M=0.3$ case.

Boundaries in the $x_{2}$ direction are taken to be periodic. In the $x_{1}$ direction, at the top and bottom of the domain, shear-free slip boundary conditions were implemented with numerical diffusion buffer zones immediately before each boundary interior to the domain. The purpose of these “open” boundary conditions is to essentially mimic an infinite domain and to ensure that both the background stratification is preserved and that none of the shocks introduced by the RTI initialization are reflected back into the domain. In particular, the buffer zones ensure that any shockwaves are dissipated prior to reaching the boundaries Reckinger2016 .

As discussed in Reckinger2016 , some artificial thickening of the interface at $x^{*}_{1}=0$ and $t^{*}=0$ can be beneficial since the thicker interface can act as a buffer layer to absorb other numerical errors. In general, however, thicker interfaces have the potential to introduce asymmetries in the initial conditions which propagate as undesirable longer-time asymmetries during RTI growth. Based on these two competing considerations, the number of points across the interface was chosen to be 16, to both minimize the asymmetry and to gain some measure of beneficial buffering effects. Finally, it was found that higher resolutions led to better initial conditions. At a level of $j_{max}=7$, it was found that the asymmetry drops below machine precision, and thus this level of resolution was deemed sufficient for the present simulations.

IV Rayleigh Taylor Instability Growth and Characteristics

The PAWCM-enabled simulations performed here are designed to allow examination of stratification strength (as parameterized by $M$) and Reynolds number (as parameterized by $Re_{p}$) effects on RTI growth and characteristics. In the following, these two effects are investigated with a primary focus on the heights and velocities of bubbles and spikes formed during the RTI development. Here the “height” is denoted $h$ and refers to the absolute value of the respective distances from $x_{1}=0$ of the bubble and spike “tips” in the $x_{1}$ direction. The bubble and spike tips correspond to the $99\%$ and $1\%$ mass fraction values, respectively. Bubble and spike velocities, denoted $u_{h}$, are computed from the time derivatives of the bubble and spike heights. An analysis of the dynamics underlying the observed bubble and spike behaviors is outlined in Section V.

IV.1 Effects of Stratification Strength

Figure 2 shows RTI growth as a function of time for each of the four stratification strengths, where $Re_{p}=20,000$ in all cases. For each case, bubbles and spikes form soon after initializing the simulation and the RTI grows as $t$ increases. Small-scale features in each case become increasingly pronounced as the RTI evolves, and secondary vortices are most prominent for the weakest stratification (i.e., $M=0.3$). The corresponding bubble and spike growths decrease as the stratification strength increases; for the strongest stratification (i.e., $M=1.2$), the RTI growth is halted relatively early in its evolution.

Consistent with the fields in Figure 2, Figure 3(a) shows that the suppression of RTI growth compared to the incompressible case from Wei & Livescu Wei2012 occurs for all stratifications considered. For the two strongest stratifications (i.e., $M=0.9$ and 1.2), the bubble and spike each reach maximum heights before $t^{*}=20$ and stop growing.

The dependence of RTI growth on stratification strength can be investigated further by considering time series of the bubble and spike tip velocities, as shown in Figure 3(b). This figure indicates that bubble and spike velocities for the larger Mach numbers all trend towards zero, indicative of the complete suppression of RTI for strong stratifications. For $M=0.3$, however, there is a re-acceleration of the spike tip shortly after $t^{*}=15$.

In addition to these changes in the bubble and spike heights with varying stratification strength, Figure 3(a) also shows that spikes reach consistently greater heights than bubbles for all $M$. This asymmetry, particularly for low $M$, is not present in the purely incompressible case of Wei & Livescu Wei2012 , where it was found that for the low Atwood number case of 0.04, bubble and spike heights were close until after the re-acceleration regime. As indicated by Figure 3(b), the velocities at the tips of the spikes are consistently larger than those at the tips of the bubbles, although the difference between these velocities becomes significant only for the $M=0.3$ case after $t^{*}>10$.

It should be noted that full suppression of RTI growth for all but the lowest value of $M$ cannot be predicted based solely on considerations of the potential energy of the system. This is shown in Figure 3(b), where only the lowest value of $M$ reaches a plateau near the velocity predicted from either drag Birkhoff1960 ; Oron2001 or potential flow Goncharov2002 models (namely, $u_{h}/\sqrt{g\lambda}\approx 0.063$).

Based on the DNS results for $M=0.3$ to 1.2, the primary observations are that larger stratifications are associated with decreasing bubble and spike growth rates, resulting in a suppression of the RTI for all but the smallest value of $M$ studied here, and that smaller stratifications are associated with more asymmetric bubble and spike growth rates. This amounts to an anomalous asymmetry at low stratifications (i.e. $M=0.3$), since both zero and large $M$ limits are more symmetrical. The results concerning the suppression of the instability are in general agreement with those from prior studies Lafay2007 ; Gauthier2010 ; Wei2012 ; Reckinger2016 ; Gauthier2017 . To better understand the dynamics leading to RTI suppression and the development of bubble and spike asymmetries, an analysis of the underlying vorticity dynamics is performed in Section V.

IV.2 Effects of Reynolds Number

Figure 4 shows RTI growth for the weakest (i.e., $M=0.3$) and strongest (i.e., $M=1.2$) stratifications for perturbation Reynolds numbers $Re_{p}=5,000$, $10,000$, and $20,000$. There is little qualitative dependence of the bubble and spike heights on $Re_{p}$, indicating that these large-scale characteristics of RTI growth are already in an asymptotic limit for $Re_{p}=5,000$. This is consistent with the results from Wei & Livescu Wei2012 , where it was found that there is little difference in the RTI growth rates before the onset of the very late chaotic development for values of $Re_{p}$ above roughly $1,500$.

Despite the relative similarity of the large-scale structure for the three values of $Re_{p}$ examined here, however, there is substantial dependence of small-scale structure on $Re_{p}$. In particular, Figure 4 shows that an increasing amount of small-scale detail emerges as $Re_{p}$ increases, corresponding to the occurrence of viscous dissipation at increasingly smaller scales. This increase in scale range with increasing $Re_{p}$ results in the formation of secondary vortices for $M=0.3$. Even though there is also increasing small scale structure for $M=1.2$ with increasing $Re_{p}$, the formation of secondary vortices is less pronounced for this higher stratification due to the overall suppression in the RTI growth.

From a quantitative perspective, Figure 5 shows bubble and spike heights and velocities for each of the four values of $M$ examined in the present study. For the bubble and spike heights shown in Figure 5(a), there is little or no dependence on $Re_{p}$ for any stratification strength. However, for the velocities in Figure 5(b), there is a clear trend towards faster initial accelerations as $Re_{p}$ increases. The bubbles and spikes also reach larger maximum velocities as $Re_{p}$ increases. However, at very early times in the evolution for each $M$, during diffusive growth, bubble and spike velocities are largest for small $Re_{p}$, eventually crossing over in each case at $t^{*}\approx 5$ such that the higher $Re_{p}$ cases have greater velocities at later times. This result is consistent with the crossover in speeds observed by Wei & Livescu Wei2012 and, to a somewhat lesser extent, by Gauthier Gauthier2017 .

These trends are consistent for all stratification strengths, although the differences with $Re_{p}$ become more pronounced as $M$ increases. For example, the peak bubble and spike velocities for $M=1.2$ are reached at roughly $t^{*}=12$ when $Re_{p}=20,000$ and at roughly $t^{*}=14$ when $Re_{p}=5,000$. After reaching the peak values, however, the bubble and spike velocities become substantially less dependent on $Re_{p}$. For the case with smallest $M$, the results approach the nearly incompressible limit where, as shown by Wei & Livescu Wei2012 , no dependence on $Re_{p}$ is observed above $Re_{p}\approx 1,500$ during the times examined here (before the onset of late time chaotic development regime).

Taken together, these results indicate that, for the values of $Re_{p}$ examined here, there is little dependence of the global RTI growth on $Re_{p}$ during the later stages of the instability at higher stratifications and through the early re-acceleration stage for $M=0.3$. However, the early time evolution, small scale structure, and the appearance of secondary vortices are all substantially affected by $Re_{p}$. Given the increasing effect of $Re_{p}$ with increasing $M$, it may be the case that $Re_{p}$ effects become increasingly pronounced for even stronger stratifications than the $M=1.2$ case examined here; exploring such more strongly stratified scenarios is left as a direction for future research.

V Vorticity Dynamics for Compressible Rayleigh Taylor Instability

Properties and dynamics of the vorticity vector, $\omega_{i}=\epsilon_{ijk}\partial u_{k}/\partial x_{j}$, where $\epsilon_{ijk}$ is the alternating tensor, have been widely studied to understand flow behavior in a variety of contexts. For compressible flows more specifically, vorticity has been studied in shock-driven Larsson2009 ; Larsson2013 ; Livescu2016 , reacting Steinberg2009 ; Hamlington2011 , and various types of buoyant Gauthier2008 flows, revealing the dynamical importance of variable density effects such as dilatation and baroclinic torque. In the case of RTI, however, only Gauthier Gauthier2017 has examined vorticity dynamics in the fully compressible regime, and for only one value of the initial stratification strength.

In the following sections, properties of the vorticity during RTI growth are outlined as a function of stratification strength, and terms in the non-dimensional compressible vorticity transport equation are subsequently examined to understand the underlying dynamics. The role of baroclinic torque, in particular, in the suppression of RTI growth for strong stratifications and in the formation of bubble and spike asymmetries for weak stratifications is outlined. It should be noted that the importance of baroclinic torque in RTI growth is not new or surprising and has been highlighted in several previous studies Lafay2007 ; Schneider2016a ; Gauthier2013 ; Gauthier2017 ; Wei2012 . The primary contribution of the current work is in explaining how the baroclinic torque varies with initial stratification strength, as well as how RTI suppression and asymmetry arise from a dynamical perspective.

V.1 Vorticity Evolution for Compressible RTI

In the 2D simulations, $\omega_{3}^{*}$ is the only nonzero component of the vorticity, and Figure 6 shows the temporal evolution of this component for each of the stratification strengths. For each stratification strength, the vorticity field initially develops as a vortex pair with generally positive vorticity for $x_{2}<0.5$ and negative vorticity for $x_{2}>0.5$. These initial vortex pairs evolve by moving downwards slowly in the domain, while the Kelvin-Helmholtz instability on the sides of the bubbles and spikes sheds further vortex pairs. The overall spatial extent of vorticity production is greatest for weak stratification (i.e., $M=0.3$), with “fronts” of non-zero vorticity magnitude that propagate upwards and downwards in an analogous way to the propagation of bubbles and spikes, respectively, as shown in Figure 2. The vorticity evolution at $M=0.3$ is reminiscent of the overall picture in the incompressible case, with induced vortical velocity helping the instability growth and leading to re-acceleration and late time chaotic development. However, at higher $M$ values, no additional vortex pairs are generated. The overall magnitude of the vorticity is also shown in Figure 6 to decrease with increasing $M$.

The overall $M$-dependence of the vorticity magnitude is also explored in Figure 7 using the vorticity averaged over the half domain, denoted $\overline{\overline{\omega}}_{3}$, where the half-domain averaging operator is defined for an arbitrary quantity $f$ as

[TABLE]

Figure 7 shows that $\overline{\overline{\omega}}_{3}$ generally increases at early times at a rate that is larger with decreasing stratification. After the initial growth of $\overline{\overline{\omega}}_{3}$ shown in Figure 7, the average vorticity decreases with time for all but the weakest stratification (i.e., $M=0.3$). This result mirrors the suppression of RTI growth for all but the weakest stratification, seen in Figure 3.

V.2 Non-Dimensional Compressible Vorticity Transport Equation

The dynamics governing vorticity evolution in compressible RTI can be understood from the non-dimensional vorticity transport equation, which reveals the explicit dependence of the dynamics on $A$, $M$, and $Re_{p}$. A similar equation was derived using the Boussinesq approximation by Schneider & Gauthier Schneider2016a , although any explicit dependence on the initial stratification strength was omitted in the derivation. Here, the non-dimensional transport equation is derived for the fully compressible case, permitting the explicit identification of dependencies on stratification strength $M$.

By taking the curl of the momentum equation in Eq. (2), the transport equation for the 3D vorticity vector is obtained for a variable density, variable viscosity compressible flow as

[TABLE]

where $D/Dt\equiv\partial/\partial t+u_{i}\partial/\partial x_{i}$ is the Lagrangian derivative and $v\equiv 1/\rho$ is the specific volume, which is used here instead of $\rho$ to simplify the derivation. The first term on the right-hand side of Eq. (27) represents vortex stretching, the second term represents dilatation, the third term is the baroclinic torque, and the last term is viscous diffusion, where the viscous stress tensor $\tau_{kl}$ has been expressed in terms of the deviatoric strain rate tensor $S^{\prime}_{kl}$ [see Eq. (9)].

The last term in Eq. (27), representing viscous diffusion, can be separated into an essentially incompressible term that is present for all $M$, regardless of whether viscosity $\mu$ is spatially and temporally varying, and into a term that is only present when $\mu$ is non-constant. In the present simulations, $\mu=\nu\rho$, where $\nu$ is a constant given in terms of problem parameters as in Eq. (17), and $\rho$ is the spatially and temporally varying density. Expansion of the diffusion term in Eq. (27) then gives the vorticity transport equation for a variable density, variable viscosity flow as

[TABLE]

The first two terms on the right-hand side of this equation are present even in constant density, constant viscosity flows, while the last three terms are only nonzero when $v$ (and, by extension, the density) is non-constant.

Using the characteristic time scale $\sqrt{\lambda/g}$ to define the non-dimensional vorticity $\omega_{i}^{*}\equiv\omega_{i}\sqrt{\lambda/g}$, and using $\rho_{I}g\lambda$ as the characteristic pressure, Eq. (28) can be written in non-dimensional form as

[TABLE]

Based on the above equation, both diffusive terms thus scale in an identical fashion with $A$ and $Re_{p}$. It should be noted, however, that the stratification strength $M$ does not appear explicitly in Eq. (29), although it is present implicitly in the baroclinic torque term [i.e., the third term on the right-hand side of Eq. (29)]. To reveal this dependence, this term can be rewritten by defining new perturbation variables $v^{\prime*}$ and $p^{\prime*}$ that express $v^{*}$ and $p^{*}$ relative to their respective $A=0$ background stratifications as

[TABLE]

where $\textbf{x}^{*}=[x^{*}_{1},x^{*}_{2},x^{*}_{3}]$, $v_{M}^{*}\equiv 1/\rho_{M}^{*}$, and $\rho_{M}^{*}$ and $p_{M}^{*}$ correspond to the $A=0$ profiles of $\rho^{*}_{0}$ and $p_{0}^{*}$ from Eqs. (23) and (24), respectively. The $A=0$ profiles are used for normalization purposes to avoid discontinuities in the derivatives of the background profiles that arise when $A$ is nonzero (particularly for the first derivative of $p_{0}^{*}$). The resulting $A=0$ profiles are, nevertheless, not substantially different than the $A=0.04$ profiles (see Figure 1) and serve the purpose of explicitly revealing the dependence of the baroclinic torque on $M$.

The decompositions in Eq. (30) are designed to yield $v^{\prime*}=v^{*}$ and $p^{\prime*}=p^{*}$ in the limit as $M\rightarrow 0$, as well as $\partial v^{\prime*}/\partial x^{*}_{i}=\partial v^{*}/\partial x^{*}_{i}$ and $\partial p^{\prime*}/\partial x^{*}_{i}=\partial p^{*}/\partial x^{*}_{i}$ in the same limit. The baroclinic torque term in Eq. (29) depends only on gradients of $v^{*}$ and $p^{*}$, and the equivalency of the perturbation and total gradients can be shown for $M\rightarrow 0$ as

[TABLE]

where $v_{M}^{*}\rightarrow 1$ and $M^{2}p_{M}^{*}\rightarrow 1$ as $M\rightarrow 0$. As $M$ becomes large and background stratification becomes increasingly strong, the magnitude of $\partial v^{\prime*}/\partial x^{*}_{i}$ becomes increasingly small and $\partial p^{\prime*}/\partial x^{*}_{i}$ approaches the background stratification everywhere, as shown in Figure 8.

Using Eqs. (31) and (32), it can be shown that the baroclinic torque on the right-hand side of Eq. (29) can be written as

[TABLE]

where it is assumed that $v_{M}^{*}$ and $p^{*}_{M}$ depend only on $x_{1}^{*}$. The first term on the right in Eq. (33) represents the baroclinic torque that is independent of the initial background stratification, and this is the only remaining term in the limit as $M\rightarrow 0$. The second and third terms represent the baroclinic torques associated with the initial stratified background pressure and specific volume fields, respectively. It should be noted that the present analysis is specific to the isothermal forms for $v_{M}^{*}$ and $p^{*}_{M}$ obtained from Eqs. (31) and (32), and that the scaling may differ for different initial background conditions (e.g., isentropic or isobaric conditions).

After substituting Eq. (33) into Eq. (29), the non-dimensional 3D vorticity transport equation is obtained for a compressible flow with initial background stratification as

[TABLE]

where, once more, the first four terms are present even in the limit as $M\rightarrow 0$ and the fifth and sixth terms are only significant for nonzero $M$. The corresponding transport equation for the vorticity magnitude $\omega^{*}\equiv(\omega_{i}^{*}\omega_{i}^{*})^{1/2}$ is given by

[TABLE]

where $\widehat{\omega}^{*}_{i}\equiv\omega^{*}_{i}/\omega^{*}$ is the vorticity unit vector (where the magnitude of $\widehat{\omega}^{*}_{i}$ is unity by definition). This expression is valid for any $M$, $A$, and $Re_{p}$ provided that $v_{M}^{*}$ and $p_{M}^{*}$ are given by Eqs. (31) and (32) and that $\nu$ is constant. In the above expression, $\mathcal{T}^{*}_{1}$ represents production and destruction of $\omega^{*}$ due to vortex stretching, $\mathcal{T}^{*}_{2}$ represents diffusion of vorticity by constant viscosity, $\mathcal{T}^{*}_{3}$ represents dilatational effects, $\mathcal{T}^{*}_{4}$ represents stratification-independent baroclinic torque, $\mathcal{T}^{*}_{5}$ represents baroclinic torque production associated with the background pressure field, $\mathcal{T}^{*}_{6}$ represents baroclinic torque production associated with the background specific volume field, and $\mathcal{T}^{*}_{7}$ represents diffusion associated with variable viscosity. In the limit as $M\rightarrow 0$, both $\mathcal{T}^{*}_{6}$ and $\mathcal{T}^{*}_{7}$ terms go to zero.

V.3 Effects of Stratification Strength on the Dynamics of the Vorticity

Figure 9 shows fields of $\mathcal{T}^{*}_{2}$-$\mathcal{T}^{*}_{7}$ from Eq. (35) for the four different stratification strengths (with $Re_{p}=20,000$ in all cases) at a late stage ($t^{*}=20$) in the 2D simulations. The vortex stretching term $\mathcal{T}^{*}_{1}$ is identically zero in 2D and is thus not shown here.

For small $M$, Figure 9 shows that $\mathcal{T}^{*}_{4}$, representing the perturbation baroclinic torque, is primarily positive (indicating vorticity production) and roughly an order of magnitude larger than $\mathcal{T}^{*}_{5}$ (baroclinic torque due to the background pressure) and $\mathcal{T}^{*}_{6}$ (baroclinic torque due to the background specific volume). The dilatation term $\mathcal{T}^{*}_{3}$ has a similar magnitude to $\mathcal{T}^{*}_{5}$ and $\mathcal{T}^{*}_{6}$, but is smaller than $\mathcal{T}^{*}_{4}$. The constant viscosity diffusion term, $\mathcal{T}^{*}_{2}$, is much larger than the variable viscosity contribution, $\mathcal{T}^{*}_{7}$, and reaches peak magnitudes similar to, but still smaller than, $\mathcal{T}^{*}_{4}$.

Figure 9 thus indicates that the primary dynamical effects for low $M$ are the perturbation baroclinic torque (i.e., $\mathcal{T}^{*}_{4}$) and constant viscosity diffusion (i.e., $\mathcal{T}^{*}_{2}$), although the former dominates the latter, resulting in the growth of the instability for low $M$. The relative magnitudes of these terms are shown in Figure 10, where the terms $\mathcal{T}^{*}_{i}$ from Eq. (35) are averaged over half of the domain along the $x_{2}$ direction to give $\overline{\mathcal{T}}^{*}_{i}$ as a function of $x_{1}$, with the average defined as

[TABLE]

Figure 10 also shows results for the averages of $\overline{\mathcal{T}}^{*}_{i}$ over $x_{1}$ for $x_{1}<0$ and $x_{1}>0$. For the weakest stratification examined here, Figure 10(a) shows that the enstrophy is created on average due almost entirely to the perturbation baroclinic torque. There is only a relatively small enstrophy destruction contribution due to the constant viscosity diffusion.

Although the perturbation baroclinic torque $\mathcal{T}^{*}_{4}$ can become locally negative due to density inversions (i.e. negative density gradients) created by vortical motions, Figures 9 and 10 show that this term remains mostly positive for all but the strongest stratification, due to the presence of the instability. Nevertheless, $\mathcal{T}^{*}_{4}$ does decrease in magnitude as $M$ increases and, in particular, Figure 10(d) shows that this term can contribute to the destruction of vorticity magnitude for sufficiently large stratification. This is consistent with Figure 8, which shows that $\partial v^{\prime*}/\partial x^{*}_{i}$ approaches zero, while $\partial p^{\prime*}/\partial x^{*}_{i}$ becomes close to 1, as $M$ increases. The reduced vorticity production at larger stratifications corresponds to the suppression of the instability growth.

For all values of $M$ considered here, Figures 9 and 10 show that $\mathcal{T}^{*}_{4}$ has the largest magnitude in general, but, since it decreases with $M$, becomes more similar in magnitude to the other terms at the largest stratification considered. At $M=1.2$, the vorticity production is much smaller, consistent with the overall suppression of the instability. For large stratifications, the reduced vorticity magnitude also translates into lower self-propagating velocity for the vortex pairs generated at the bubble/spike interface. In turn, this results in part of the fresh fluid brought towards the bubble/spike peaks by the induced vortical velocity returning back to the mixing layer. Thus, at $M=0.9$, Figure 10(c) shows density inversions (i.e., negative $\overline{\mathcal{T}}^{*}_{4}$ or stabilizing regions) near the edges of the layer, while at $M=1.2$ in Figure 10(d) these regions can occur throughout the layer.

Variations in the magnitudes of $\mathcal{T}^{*}_{5}$ and $\mathcal{T}^{*}_{6}$ with $M$ shown in Figures 9 and 10 are somewhat more complicated. In particular, the peak magnitudes of both terms increase from $M=0.3$, but start decreasing again at larger $M$ and become smaller for $M=1.2$. Term $\mathcal{T}^{*}_{5}$ reaches its peak magnitude at slightly smaller Mach number than $\mathcal{T}^{*}_{6}$ ($M\sim 0.6$ versus $M\sim 0.9$). At large $M$, the gradient contributions to both terms become uniform, with values of 0 and 1, respectively (see Figure 8). The prefactors $(M^{2}p_{M}^{*}-1)$ and $M^{2}v_{M}^{*}$ depend on the initial background stratification and become large in the far-field, but are small near the centerline, even for large $M$. Therefore, as the instability growth is suppressed at large stratifications, $\mathcal{T}^{*}_{5}$ and $\mathcal{T}^{*}_{6}$ are more confined to the region close to the centerline and never reach regions with large prefactor values.

Perhaps most significantly, both $\mathcal{T}^{*}_{5}$ and $\mathcal{T}^{*}_{6}$ exhibit asymmetries that affect the overall growth of the instability. Aside from local inversions, Figures 9 and 10 show that $\mathcal{T}^{*}_{5}$ presents a top-bottom asymmetry with respect to the $x_{1}=0$ initial location of the instability (i.e., it is positive on the spike side and negative on the bubble side). On the other hand, Figure 9 shows that $\mathcal{T}^{*}_{6}$ presents a left-right asymmetry with respect to the interface between the heavy and light fluid, with negative values inside the spike and positive values inside the bubble regions. Conversely, the dilatation term $\mathcal{T}^{*}_{3}$ shows the opposite left-right asymmetry, with positive values inside the spike and negative values inside the bubble regions. It should be noted, however, that $\mathcal{T}^{*}_{5}$ becomes larger with respect to $\mathcal{T}^{*}_{4}$ as $M$ increases. The term $\mathcal{T}^{*}_{4}$ is itself also asymmetric, as shown in Figure 10, and it is likely that this asymmetry is the underlying cause of the differences in bubble and spike growth rates, particularly for small $M$. Moreover, weak asymmetry in bubble and spike growth rates is observed in the incompressible limit, indicating that $\mathcal{T}^{*}_{5}$ may be a contributor to, but not the sole cause of, the asymmetry, since this term approaches zero as $M\rightarrow 0$.

As explained above, the bubble-spike asymmetry is small in the incompressible case before the chaotic stage, it becomes noticeable at $M=0.3$, and then decreases again at large stratifications. The history of the top-down asymmetry in the vorticity generation can also be seen from the time evolutions of $\overline{\mathcal{T}}_{i}$ in Figure 11. This figure does not identify the left-right asymmetry, which has a more dynamical effect, as it influences the vortical motions separately within the bubble and spike regions. However, it does show that $\overline{\mathcal{T}}^{*}_{4}$ begins symmetrical and develops the top-down asymmetry at some later time. On the other hand, $\overline{\mathcal{T}}^{*}_{5}$ is asymmetric from the beginning, such that it represents the source of this asymmetry. This is consistent with the incompressible flow results, where terms $\overline{\mathcal{T}}^{*}_{3}$, $\overline{\mathcal{T}}^{*}_{5}$, and $\overline{\mathcal{T}}^{*}_{6}$ are zero, and $\overline{\mathcal{T}}^{*}_{4}$ remains relatively symmetrical until later times. Again, at large stratifications, the overall reduction in vorticity production and suppression of the instability prevents the bubble/spike asymmetry from becoming more pronounced.

Finally, Figures 9–11 show that the variable viscosity diffusion term, $\mathcal{T}^{*}_{7}$, is negligible for all $M$ considered here and all times. By contrast, the magnitude of the constant viscosity diffusion term, $\mathcal{T}^{*}_{2}$, remains more uniform with increasing $M$, although it does become more consistently negative as $M$ increases, as shown most clearly in Figure 11. This indicates that constant viscosity diffusion is the primary term leading to destruction of vorticity magnitude, and this term begins to rival the magnitude of the perturbation baroclinic torque term (i.e., $\mathcal{T}^{*}_{4}$) for large $M$.

Taken together, these results are indicative of larger vorticity production within the spike region, as compared to the bubble region, due to compressibility and stratification effects. Because the bubble and spike vertical axes are maintained throughout the flow evolution for the single mode case, the vorticity field itself retains a similar symmetry. This results in an induced vortical velocity along the bubble/spike axes, which helps the instability grow, similar to the incompressible case Wei2012 . However, for the compressible case, the dilatation term and baroclinic contributions sum up to a bubble/spike asymmetry even at low Mach numbers. At higher Mach numbers, due to the overall suppression of the instability, these contributions also decrease and the asymmetry becomes small again. Overall, the primary dynamical balance is between constant viscosity diffusion, which leads to the destruction of vorticity magnitude, and perturbation baroclinic torque, which leads to vorticity magnitude production.

VI Conclusions

In the present study, wavelet-based adaptive mesh refinement has been used to perform DNS of 2D single-mode compressible low Atwood number RTI for four different isothermal stratification strengths, corresponding to Mach numbers from 0.3 to 1.2, and for three different perturbation Reynolds numbers from $5,000$ to $20,000$. The simulation results have been examined to understand the effects of stratification strength and Reynolds number on the characteristics, dynamics, and rate of RTI growth. In the present context, compressibility is controlled through the values of the background pressure at the interface between the heavier and lighter fluids, which also affects the background stratification strength, and would be considered flow, as opposed to fluid, compressibility. In this context, the incompressible limit is reached as the speed of sounds goes to infinity by increasing the interface pressure and temperature, such that the interface density remains constant. The practical setup corresponds to an enclosed fluid system that is uniformly heated (i.e., heating at constant volume).

For weak stratifications, RTI growth was found to undergo a re-acceleration after reaching a plateau in the growth rate that approximately matched predictions from potential flow theory. As the stratification strength increased, however, this re-acceleration was found to no longer occur, and the RTI growth was suppressed; this suppression occurred in the present study for all Mach numbers greater than 0.3. For weak stratifications, the bubble was found to grow at a slower rate than the spike, but this asymmetry progressively weakened as the stratification strength increased. The Reynolds number was found to have little impact on RTI growth for the range of Mach numbers and for the simulation length examined here. However, small-scale structure was found to become more pronounced as the Reynolds number increased. At very early times, during the diffusive stage, the growth rates were larger at smaller Reynolds numbers, but the instability became faster during the linear and weakly nonlinear stages at higher Reynolds numbers, consistent with prior studies of Reynolds number effects Wei2012 ; Gauthier2017 .

To determine the origins of the observed results, the dynamics of the vorticity magnitude were examined in detail. A non-dimensional compressible vorticity transport equation was derived to explicitly show dependencies on the Mach, Atwood, and Reynolds numbers, and the effects of stratification strength were studied for each of the terms in the transport equation. This analysis showed that incompressible baroclinic torque was the dominant driver of RTI growth for the range of stratifications considered, and its decrease at higher stratifications corresponded to the overall instability suppression. Asymmetries in the RTI growth were found to be the result of compressibility effects, as a consequence of the dilatation term and background stratification contributions to the baroclinic torque. However, for strong stratifications, since the instability did not evolve far from the centerline, the latter contributions remain small and the bubble/spike asymmetry does not become pronounced.

In total, the simulations and analysis performed in this study have enabled the three questions posed in Section I to be fully addressed. However, much work remains to be done. In particular, the present analysis of vorticity dynamics should be extended to multi-mode initial perturbations, to different stratification types (e.g., isopycnic and isentropic stratifications), and to 3D cases where vortex stretching effects in the vorticity dynamics are nonzero. It would also be of interest to explore longer simulation times for the weakly stratified cases to determine whether the chaotic development regime noted by Wei & Livescu Wei2012 is recovered in the context of fully compressible simulations.

VII Acknowledgements

This work was made possible in part by funding from the LDRD program at Los Alamos National Laboratory through project number 20150568ER. SAW was supported by Los Alamos National Laboratory, under Grant No. 316898. Computational resources were provided by the LANL Institutional Computing (IC) Program. Dr. Oleg Vasilyev provided inspiration and early direction for this work.