## Abstract

We present a tentative review of compressibility effects in Rayleigh–Taylor instability-induced flows. The linear, nonlinear and turbulent regimes are considered. We first make the classical distinction between the static compressibility or stratification, and the dynamic compressibility owing to the finite speed of sound. We then discuss the quasi-incompressible limits of the Navier–Stokes equations (i.e. the low-Mach number, anelastic and Boussinesq approximations). We also review some results about stratified compressible flows for which instability criteria have been derived rigorously. Two types of modes, convective and acoustic, are possible in these flows. Linear stability results for perfect fluids obtained from an analytical approach, as well as viscous fluid results obtained from numerical approaches, are also reviewed. In the turbulent regime, we introduce Chandrasekhar’s observation that the largest structures in the density fluctuations are determined by the initial conditions. The effects of compressibility obtained by numerical simulations in both the nonlinear and turbulent regimes are discussed. The modifications made to statistical models of fully developed turbulence in order to account for compressibility effects are also treated briefly. We also point out the analogy with turbulent compressible Kelvin–Helmholtz mixing layers and we suggest some lines for further investigations.

## 1. Introduction

The superposition of a heavy fluid above a lighter one in a constant acceleration field is a potentially unstable hydrodynamic configuration called the Rayleigh–Taylor instability (RTI; Rayleigh 1883; Taylor 1950; Chandrasekhar 1961). Such instabilities are present in various physical situations and play a prominent role, for instance, in inertial-confinement fusion (Lindl 1998; Atzeni & Meyer-Ter-Vhen 2004) and in supernova explosions in astrophysics (e.g. Zingale *et al.* 2005). In addition to these applications, the RTI remains a basic problem in fluid dynamics, despite the fact that it has been studied for a number of decades. This instability belongs to the large class of buoyant instabilities which, in its simplest form, occurs between two inviscid incompressible fluids. However, the presence of transport coefficients (viscosity, thermal conduction and species diffusion), surface tension, compressibility, etc. modifies the phenomenology of the instability. In situations where both the interface between the two fluids is very corrugated and irregular such that a large number of Fourier modes are present, and the viscosity coefficient time scale is much larger than the instability time scale, the RTI may generate a turbulent mixing layer. Beyond the result of an experiment or a numerical simulation, it is worth analysing the mixing formation owing to nonlinear interactions, the structure of the mixing and the influence of various parameters such as viscosity, thermal conduction, species diffusion and compressibility effects. In some cases, such as in inertial-confinement fusion and in astrophysics, variable fluid density is expected to play a particularly significant role. It is well known that compressibility effects have two different origins. The first one, called ‘static compressibility’, is due to the variable density of the fluid. In an acceleration field, it results in stratification and leads to density gradient length scales, *L*_{ρH,L}, for the heavy and light fluids, respectively. The second one, called ‘dynamic compressibility’, is essentially an effect of the finite speed of sound. For a perfect gas equation of state (EOS), this effect is governed by the adiabatic indices *γ*_{H,L} of the two fluids. In RTI studies, the word ‘compressible’ is often used with both meanings, leading to some misunderstandings. The aim of this paper is to present a tentative review of compressibility effects in RTI-induced flows, according to the turbulent mixing and beyond (TMB) conference objectives of a better understanding of turbulence mixing and turbulence in unsteady flows. We note that compressibility effects in some buoyant flows have already been discussed in the proceedings of the TMB conference (Gauthier *et al.* 2008).

This paper is organized as follows. In §2, we review the governing equations and Kovásznay’s approach where a small amplitude motion in a compressible flow may be decomposed into three ‘modes’, i.e. the vorticity mode, the acoustic (pressure) and the entropy mode. From this discussion, the question of the limits of the Navier–Stokes (NS) equations then naturally arises. For weakly compressible flows, simplified sets of equations are available, where the acoustic contribution or the stratification has been removed, and are substantially simpler than the full original set. Before looking at the linear regime of RTI in §3, we investigate the local instability criterion for stratified flows of perfect fluids. We identify convective- and acoustic-type modes and present two examples of the latter, one in isentropic stratified atmospheres and another in imploding spherical shells. The linear regime of inviscid fluid flow is investigated and it is shown that there exist several ‘compressible RTI’, depending on the thermodynamic state, i.e. whether it is isentropic, isothermal or general, of the hydrostatic equilibrium state and the nature of the perturbations. Results for the linear regime of the RTI for compressible perfect fluids, for which analytical results exist, are described in greater detail. Regarding the behaviour of viscous fluids, the few results obtained by numerical means are reviewed. The turbulent regime of compressible RTI is investigated in §5 and we introduce Chandrasekhar’s observation on isotropic compressible turbulence, which states that the largest structures in the density fluctuations are determined by the initial conditions and represent permanent features of the flow. The nonlinear regime is investigated through numerical simulations by solving the full Euler or NS equations but most simulations are carried out in the weakly compressible limit, i.e. for weakly stratified equilibrium states and for the same value of the adiabatic indices for the two fluids. From a numerical point of view, it turns out that, despite the present computational resources available, a fully resolved simulation, i.e. down to the Kolmogorov scale, of the three-dimensional NS equations for compressible fluids, starting from rest through to fully developed turbulence, remains a challenging configuration. We conclude by presenting some features of the numerical simulations and by showing that the study of turbulent compressible Kelvin–Helmholtz mixing layers may be beneficial to the investigation of RT mixing layers.

## 2. The NS equations

### (a) The complete NS equations

The RTI concerns the acceleration of the thin layer separating two fluids of different densities. In immiscible fluids, this layer reduces to an interface, whereas for miscible fluids their densities vary abruptly over a finite volume. For such miscible fluids, one uses the complete NS equations within the single fluid approximation. These equations are written as follows:
2.1
2.2
2.3
2.4and
2.5
where *u*_{i} (*i*,*j*=1,2,3) are the three components of the velocity and *p*, *ρ* and *T* are the pressure, density and temperature of the single fluid. The specific heats at constant volume and pressure are denoted *C*_{v,p;H,L} and their ratios *γ*_{H,L}, where the indices H, L stand for the heavy and light fluids. The expression of the mixing specific heats reads *C*_{v,p;m}(*c*)=*cC*_{v,p;H}+(1−*c*)*C*_{v,p;L}, where *c* is the concentration defined below. The stress tensor *σ*_{ij}=*μ*(∂_{j}*u*_{i}+∂_{i}*u*_{j}−2/3*δ*_{ij}∂_{ℓ}*u*_{ℓ}) is defined within the Stokes approximation, where *μ* is the dynamic viscosity coefficient, and is the rate-of-deformation tensor. The reference value of the adiabatic index is *γ*_{r} and the derivative with respect to the concentration is denoted *d*_{c}. The motion takes place in a three-dimensional closed rectangular domain *Ω*=*L*_{x}×*L*_{y}×*L*_{z} (*L*_{z}≡*h*_{H}−*h*_{L}, where |*h*_{H}| and |*h*_{L}| are the thicknesses of the heavy and light fluid layers, respectively). The *z*-axis is directed upward so that the gravity, represented by the vector * g*=(00−

*g*), is negative along this direction. Previous results obtained by different authors will be given in this coordinate system. The heavy and light fluids are initially located in the upper and lower side of the domain such that 0≤

*z*≤

*h*

_{H}and

*h*

_{L}≤

*z*≤0. It is convenient to use the following units:

*L*

_{y}(the horizontal width of the domain) for length, (

*L*

_{y}/

*g*)

^{1/2}for time, the reference density is

*ρ*

_{r}=(

*ρ*

_{H}(0

^{+})+

*ρ*

_{L}(0

^{−})) and the reference temperature is a uniform temperature

*T*. The expressions of the Reynolds, Schmidt and Prandtl numbers are , and

*Pr*=

*μγC*

_{v,r}/

*κ*, where

*C*

_{v,r}is some reference value of the specific heat at constant volume. The symbols and

*κ*stand for the coefficient of diffusion of species and the thermal conductivity coefficient, respectively.

For miscible fluids, one usually uses the classical ‘partial pressures–partial densities’ mix model (Le Creurer 2005), which reads *p*=*p*_{H}+*p*_{L}=*ρT*(1+*A*−2*Ac*), *ρ*=*ρ*_{H}+*ρ*_{L} and *T*=*T*_{H}=*T*_{L}, where the Atwood number is *A*=(*ρ*_{H}(0^{+})−*ρ*_{L}(0^{−})/(*ρ*_{H}(0^{+})+*ρ*_{L}(0^{−}). The expression of the partial pressures reads , where *ρ*_{H,L} and *T*_{H,L} are the partial densities and the temperatures of the fluids. We also introduced the fluid concentration *c*=*ρ*_{H}/*ρ* and the reference concentration is equal to *c*_{r}=(1−*A*)/2.

For immiscible fluids, the set of governing equations is obtained by letting the Schmidt number, *Sc*, tend to infinity, and each fluid obey system (2.1)–(2.3) with EOS (2.5). Equation (2.4) is no longer useful. Moreover, for inviscid flows, the Reynolds number, *Re*, is set to infinity.

We also define a parameter *Sr*, called stratification parameter, by *Sr*=*L*_{y}/*L*_{sh}, where is the ‘scale height’ of the hydrostatic equilibrium state (with , where the molar weights are and is the perfect gas constant). The initial one-dimensional unstable equilibrium state (, , , and ) is found by assuming the piecewise hydrostatic equilibrium in both the heavy and the light fluids with a constant temperature. The heavy and light fluid density profiles then read , where are the density gradient length scales, which measure the stratification of the two equilibrium states. Let us note that the stratification of the two fluids is not uniquely determined by the stratification parameter *Sr*, as it also depends on the Atwood number *A*. The speed of sound of the two fluids, , also depends on the stratification parameter *Sr* and the Atwood number *A*. For unstratified configurations, the parameter *Sr* tends to zero and density gradient length scales and speeds of sound tend to infinity. A ratio of the form *L*/*L*_{sh}, where *L* is a characteristic length scale, appears in most compressible RT studies. It sometimes appears only as a condition for a limiting case but it often appears as a Mach number based on the free fall velocity. It seems more appropriate to us to interpret this ratio as a stratification parameter, as no fluid motion is involved in the equilibrium state. This characteristic length scale *L* is usually a transverse length scale, based on the average wavelength of the perturbation, or the horizontal size of the domain. This stratification parameter will appear in various forms in the present review paper.

### (b) Quasi-incompressible limits of the NS equations

As already stated in the introduction, compressibility effects have mainly two different origins, the ‘static compressibility’ or the stratification, and the ‘dynamic compressibility’, which is an effect of the finite value of the speed of sound, calculated from the EOS. The unstratified limit, *Sr*≪1, is recovered for small thicknesses *L*≪*L*_{sh} or weak acceleration *g*≪1, and the incompressible limit is recovered by letting the adiabatic index *γ* tend to infinity. On the other hand, compressibility effects can arise when the system departs from incompressible behaviour or from qualitatively different solutions. The question is how to define or characterize these effects. In previous studies, Chu & Kovásznay (1957) and Monin & Yaglom (1962) showed that small amplitude motion in a compressible flow may be decomposed into three basic modes, i.e. vorticity, acoustic and entropy motions. In the nonlinear regime, there are interactions between these three basic modes, some interactions being more efficient than other ones. For example, the acoustic–acoustic interaction may only generate the acoustic mode, whereas the interaction between acoustic and entropy modes generates the three modes. In some situations, the acoustic mode carries a small amount of energy and may be neglected, but the stratification of the initial equilibrium state may not be negligible. Moreover, the presence of acoustic waves is very restrictive for the numerical simulations. Indeed, the time step is subject to the Courant–Friedrichs–Lewy (CFL) condition for compressible flows (Courant *et al.* 1967). It is then tempting to get rid of the acoustic waves while keeping the stratification. This is the domain of the quasi-incompressible limits of the NS equations and essentially two approximations have been developed—lying between the full compressible equations and the Boussinesq approximation, the low Mach number (Rehm & Baum 1978; Paolucci 1982) and the anelastic approximations (Gough 1969). These approximations are obtained from an asymptotic analysis of the small parameter which is a function of the Mach number. This analysis provides three criteria to determine whether or not a flow has to be calculated with the full NS equations. Let us also recall that flows with variable density are not necessarily ‘compressible’ (e.g. the mixing of two incompressible flows). On the other hand, flows are characterized by dimensionless numbers reflecting the order of magnitude of the characteristic time and length scales. The Mach number *M*=*u*/*c*_{s}, where *u* is a typical velocity and *c*_{s} is the isentropic speed of sound, gives the magnitude of compressibility effects. The Froude number , where *L* is a typical length scale and *g* the external uniform acceleration, gives the ratio of inertial forces with respect to gravity or buoyancy forces. The first criterion of the quasi-incompressible limits, *γM*^{2}≪1, gives the order of magnitude of the acoustic mode while the second criterion focuses on the quantity *γM*^{2}/*Fr*^{2}=*L*/*L*_{sh}, which measures the stratification of the flow. For weakly stratified flows, typical length scales are much smaller than the scale height (*L*/*L*_{sh}≪1), and one then obtains the low Mach number approximation. In contrast, for stratified flows, the relevant limit of the NS equations is the anelastic approximation. The third criterion ensures that the hydrodynamic approximation is valid; this one may be written as *γM*^{2}/(*RePr*)=*γM*^{2}/*Pe*≪1, where *Re*, *Pr* and *Pe* are the Reynolds, Prandtl and Péclet numbers, respectively. The sets of partial differential equations obtained from such asymptotic expansions have essentially the same mathematical properties as the NS equations for incompressible flows, i.e. the hyperbolic part of the full NS has been removed. Propagation of acoustic waves is no longer possible, however, and as a consequence the time step in numerical computations is no longer subject to the CFL condition. Instead, the constraint for the time step is only based on the fluid velocity. For vanishing stratifications, the anelastic approximation leads to the Boussinesq approximation in which the velocity field is solenoidal. More rigorous derivations, using two-parameter asymptotic expansions, have been devised (e.g. Bois 2006). The low Mach number and the anelastic approximations are widely used in astrophysics and in meteorology, for example, but to the authors’ knowledge they have not been used in investigations of the RTI.

## 3. Linear stability of stratified flows

### (a) Stratified compressible fluid and local instability criteria

We start the stability analysis of compressible flows with the study of stratified equilibrium states although, strictly speaking, the RTI between two perfect fluids deals with a density jump at the interface. However, when the density gradient is finite and the profiles are continuous, stability results on stratified flows also apply. With such simpler configurations, it is possible to obtain several general results, and, in particular, local instability criteria, which are useful estimates of unstable flow regions. These criteria also provide useful approximate guides even for the stability of time-dependent viscous flows (Dufour *et al.* 1984; Gamaly 1993). Chandrasekhar first showed from the properties of the Sturm–Liouville equation, which governs the pressure perturbation in an incompressible flow, that the necessary and sufficient condition that a stratified heterogeneous fluid be stable is that *∇ρ* should be negative everywhere; moreover, if *∇ρ* should be positive anywhere inside the fluid, the stratification is unstable (Chandrasekhar 1961, ch. X, §92; for arbitrary acceleration this instability condition reads *∇p**∇ρ*<0). For compressible fluids, Landau & Lifshitz (1959, ch. 1, §4) have proved on empirical grounds that the condition for stability is *∇S*>0. Later, Gamaly et al. (1976, 1980) generalized this condition proving that, within the linear approximation of the Euler equations, the motion is stable in the region where the following condition is satisfied:3.1In a multi-dimensional configuration, the instability condition is written as3.2where *p*, *ρ* and *S* are the pressure, density and specific entropy of the fluid, and *c*_{s} is the isentropic speed of sound. The *i*th component of the gradient is ∇_{i}. Criterion (3.2) reduces to Chandrasekhar’s result for incompressible flows, where the speed of sound tends to infinity.

### (b) General local instability criteria in a stratified perfect compressible fluid

Sitt (1980) generalized this result to an arbitrary type of thermodynamic perturbation for a stratified perfect compressible fluid. He considers a Lagrangian description valid for the time-dependent flows found in inertial-confinement fusion (Dufour *et al.* 1984); however, for the purpose considered here, this point is not relevant. One simply considers a uniformly accelerated one-dimensional flow in plane geometry. The mean flow is denoted with the superscript (*o*). The pressure and density of the equilibrium state are *p*^{(o)}, *ρ*^{(o)} and the EOS is *p*^{(o)}=*p*(*ρ*^{(o)},*F*^{(o)}), where *F* stands for any thermodynamic function such as temperature and entropy. The acceleration is written as3.3where is the density gradient length scale and the iso-*F* speed of sound. The Euler equations are linearized and we seek solutions for a variable *Ψ* of the form , where *σ* is a complex number and . As a result, the pressure perturbation , for the mode *k*, satisfies the following second-order differential equation:3.4The EOS iso-*F* for the perturbations of density and pressure is written . We also suppose that *Σ*^{(o)} is constant and we define *ζ*≡−(*σ*^{2}+*gk*^{2}*Σ*^{(o)}/*σ*^{2}). It appears that the eigenvalue problem (3.4) is a self-adjoint Sturm–Liouville problem. Consequently, it admits two infinite, discrete sets of perturbation modes, the growth rate of which are given by the following expression:3.5The modes are acoustic-type modes. Indeed, in the particular case where *Σ*^{(o)} and are constants, one can show that , where is some constant. Generally, the growth rate has a real and an imaginary part leading to the possibility of overstability (Chandrasekhar 1961, ch. I, §2, p. 3). The modes are convective-type modes. If *g* *Σ*^{(o)}>0, there is a marginal stability with an oscillatory behaviour and exponentially growing if *g* *Σ*^{(o)}<0. This latter condition provides a necessary and sufficient condition for convective instability for a stratified perfect compressible fluid against iso-*F* perturbations3.6As particular cases, we have the following criteria for instability, written in a three-dimensional geometry:

— the isothermal criterion

*∇*_{i}*p*^{(o)}*∇*_{i}*T*^{(o)}>0;— the isentropic one

*∇*_{i}*p*^{(o)}*∇*_{i}*S*^{(o)}>0; and— the incompressible one, which may be obtained from the previous criterion by taking the limit of infinite adiabatic index

*γ*. One then obtains*∇*_{i}*p*^{(o)}*∇*_{i}*ρ*^{(o)}<0.

When the speed of sound is uniform, one obtains a biquadratic equation for the growth rate *σ* (Gamaly *et al.* 1976, 1980; Sitt 1980). For small wave numbers, *k*≪*χ*_{ρ}, the expression of this growth rate is3.7This is the Brunt–Väisälä frequency, which is the natural frequency of an equilibrium state defined by a density profile, *ρ*^{(o)}, and a pressure profile, *p*^{(o)}.

#### (i) Adiabatic stratified atmospheres

Scannapieco (1981) studied a particular configuration, namely the planar exponentially stratified adiabatic atmosphere3.8where *ρ*_{o}=*ρ*^{(o)}(*z*=0) and *H* denotes here the density gradient length scale (equation (3.3)). Stability analysis is performed through a normal mode analysis leading to the second-order differential problem (3.4). He then makes the approximation that over an interval *d* such that |*H*|/*d*≪1, the isentropic speed of sound is uniform. With this restrictive hypothesis, problem (3.3) becomes a constant coefficient second-order differential equation, which can be solved explicitly. This equation admits three physically distinct types of modes that are identified as acoustic, gravity and Lamb modes (Cox 1980; Scannapieco 1981). The behaviour of these three modes depends on the value of a stratification parameter, denoted here by , and defined as , where *T* is a reference value defined by the uniform speed of sound *c*_{s}. It turns out that the acoustic mode is stable for all values of , whereas the gravity mode is growing for and . The Lamb mode, which is an acoustic-type mode, is growing-oscillating, i.e. overstable, except for the particular value . The Lamb modes are an example of the previous acoustic-type modes given by equation (3.5).

#### (ii) Unstable acoustic-type waves in imploding spherical shells

Prior to the works just quoted, Book (1978) found unstable acoustic-type waves in imploding spherical shells. He solved the Euler equations for a perfect gas in spherical geometry, with an initial density *ρ*(*r*) in the region defined by *r*_{−}≤*r*≤*r*_{+}. The class of solutions of this system is known as homogeneous self-similar motion. Both the homentropic motion and the uniform density cases have been solved. In such cases, the eigenvalue problem takes the form of a Sturm–Liouville problem with non-constant coefficients whose solutions are given in terms of confluent geometric functions. From the properties of the spectrum, he found that these acoustic waves are unstable, and moreover long-wavelength modes have the fastest relative amplification. These modes are obviously related to the unstable acoustic-type waves formally described by Sitt (1980).

## 4. Linear regime of the RTI

### (a) Linear regime of the RTI for perfect compressible fluids

A wide range of hydrostatic equilibrium states are solutions of the NS/Euler equations, depending on the type of EOS, the functional form of the thermal conductivity and the boundary conditions. Different equilibrium states may lead to different stability properties. So far, only the isothermal and isentropic equilibrium states have been investigated. Let us review these contributions. The first attempt of Vandervoort (1961) was not validated (Plesset & Prosperetti 1982; Shivamoggi 1982). Mitchner & Landshoff (1964) provided an expression for the growth rate at the limit of small wavelength disturbances. Plesset & Hsieh (1964) obtained a dispersion relation, under the condition , where *λ* is the wavelength of the perturbation. As the dimensionless number, *g* *λ*/*c*^{2}_{H,L}, is of the form of a stratification number, their results hold only for weakly stratified equilibrium states. The dispersion relation is written under the form of a correction of the classical RT growth rate , depending on the difference in the speeds of sound, *c*^{2}_{H,L}, in the two fluids. Blake (1972) first gave, though with no details, the dispersion relation for the isothermal case and found that compressibility stabilizes the instability. Mathews & Blumenthal (1977) also established the dispersion relation for isothermal equilibrium state and isothermal perturbations. Baker (1983) also gave the dispersion relation and found that ‘compressibility can either enhance or decrease the growth rate depending principally upon whether the lighter fluid is more or less compressible’. Bernstein & Book (1983) also investigated the isothermal case and found that the growth rate decreases with increasing adiabatic index *γ*. The variation is expected to be ‘a few tens of percent’ for materials of interest in inertial confinement fusion. Lezzi & Prosperetti (1989) investigated the behaviour of isentropic perturbations from an isentropic hydrostatic equilibrium state. They provided a careful parameter study the main result of which is that increasing the compressibility of the heavy fluid (decreasing *γ*) stabilizes the instability whereas increasing compressibility of the light fluid destabilizes the flow. They also found that ‘compressibility has a stabilizing effect at small wavelengths and a destabilizing effect at long wavelengths. The magnitude of these effects is, however, small in most circumstances.’ They also proved that the growth rate is real. More recently, Livescu (2004) analysed the system of two immiscible compressible ideal fluids in the presence of surface tension. He recognized three parameters, the equilibrium pressure at the interface and the adiabatic indices *γ*_{H,L}. As this equilibrium pressure increases, the flow is less compressible and the growth rate increases, while as *γ* increases, the fluid is less compressible and the growth rate decreases. He also showed an interesting finite size domain effect: for very small sizes of the heavy fluid with respect to the wavelength of the initial perturbation (*L*_{H}≪1), the growth rate may become larger than the constant density incompressible growth rate, provided that the adiabatic index *γ*_{L} is close to 1, and the stratification smaller than a critical value.

These results are summarized in table 1, where LL denotes a case where the first symbol stands for the thermodynamic state of the hydrostatic equilibrium and the second one for the perturbation type, with L=S, T or G for the isentropic, isothermal and general cases, respectively. Particular perfect fluid RTI configurations are obtained from equations (2.1)–(2.5), by taking limiting values of the transport coefficients, whereas G cases may be obtained, for example, with variable transport coefficients. Three perfect fluid cases (SS, TT and TS) have been solved analytically. Written in a dimensionless form, the dispersion relations depend on the Atwood number *A*, a stratification parameter *Sr*—attached to the hydrostatic density profile—the compressibility parameters *γ*_{H,L} (except for the TT case)—attached to the fluid—and the thicknesses of the fluid layers, *h*_{H,L}. The SS case is obtained in the limit of infinite Reynolds numbers. Compressibility effects are studied only through the EOS’s parameters as the stratification is not recognized as a parameter (Lezzi & Prosperetti 1989). The TT case (Blake 1972; Mathews & Blumenthal 1977; Baker 1983; Bernstein & Book 1983; Newcomb 1983; Book 1986; Ribeyre *et al.* 2004) is obtained for *Re* *Pr*≠0 and we have the following inequality for the linear growth rate: . The TS case (Bernstein & Book 1983; Livescu 2004) is obtained for a finite value of the thermal conductivity, and vanishing transport coefficient values for the perturbations. We have, for infinite domains, the following order relationship:4.1

From these previous studies, it appears that (i) the growth rate strongly decreases as the stratification increases; (ii) the growth rate increases as the adiabatic indices *γ*_{H,L} increase; (iii) stratification and compressibility effects are more important at small wavenumbers; (iv) growth rates are larger when the light fluid is more compressible than the heavy one (*γ*_{L}<*γ*_{H}); and (v) compressibility effects are larger at small Atwood numbers.

#### (i) Method of solution

The normal mode analysis is applied in order to obtain the dispersion relation, i.e. a relationship between the growth rate *σ* and the wavenumber *k*. A two-dimensional solution may be sought, as the linear stability problem depends only on the norm of the wavenumber of the perturbation, . One may also use the hypothesis of irrotational velocity vector field (Landau & Lifshitz 1959; Plesset & Hsieh 1964; Lezzi & Prosperetti 1989). The first step is to solve the linear equation for the perturbation in each fluid before imposing a matching condition at the interface, which gives the dispersion relation. The pressure perturbations are solution of a second-order ordinary differential equation of the form4.2where the function *C*_{i}(*z*) depends on the particular hydrostatic equilibrium state T, S or G. For a non-uniform speed of sound *c*_{i}, the solution is a linear combination of confluent hypergeometric functions. For particular cases where the speed of sound is uniform, the solution is a combination of exponential functions. The matching condition at the interface may be imposed either by assuming the continuity of the vertical velocity and the pressure at the interface (Mathews & Blumenthal 1977), or equivalently by integrating the equation over a small domain around the interface (Chandrasekhar 1961; Livescu 2004). The dispersion relation is usually written in the form below. For example, for the TT case in an unbounded domain (Mathews & Blumenthal 1977), one has4.3where the isothermal speeds of sound are . The expression of a physical variable is , where *A*_{H/L} and *B*_{H/L} are four constants and the four generalized wavenumbers read4.4The SS and TS cases (Lezzi & Prosperetti 1989; Livescu 2004) may be recast in a similar form. An attempt at solving the three-layer configuration (Hoshoudy 2007) has to be confirmed (Livescu 2008). In that case, one has two interfaces where matching conditions have to be imposed.

### (b) Linear regime of the RTI for viscous compressible fluids

There are actually very few results available regarding the linear regime of the RTI for viscous compressible fluids. One of these deals with immiscible fluids (Livescu 2004) and one with miscible fluids (Lafay *et al.* 2007). Livescu (2004) numerically solves the linear system obtained from equations (2.1)–(2.3) and (2.5), for finite values of the Reynolds number, *Re*, but infinite values of the Schmidt number, *Sc*, corresponding to immiscible fluids. It results in a single fourth-order ordinary differential equation for the vertical velocity, which reads4.5where *D* is the first derivative operator and the coefficients *A*_{i} depend on the hydrostatic equilibrium state. The previous equation is solved with a Runge–Kutta method and the matching condition at the interface is solved with Broyden’s method. As expected, finite Reynolds number introduces a cut-off wavenumber in the dispersion curve and, consequently, a wavenumber at which the growth rate is maximum. The author notes a qualitative difference between the compressible and incompressible cases at small Atwood numbers. For compressible configurations, the critical wavenumber increases with the Atwood number, as opposed to the incompressible case. Lafay *et al.* (2007) numerically solve the linear system derived from the nonlinear system (2.1)–(2.5) for finite values of the Reynolds number, *Re*, the Schmidt number, *Sc*, and Prandtl number, *Pr*, corresponding to miscible fluids. It results in a ninth-order differential equation written under the form of a generalized linear eigenvalue problem4.6where and the matrices *A* and *B* are 5×5 block matrices for two-dimensional geometry. These matrices depend on the equilibrium state and on the first- and second-derivative operators. This generalized linear eigenvalue problem is solved with an auto-adaptive domain decomposition Chebyshev method (Gauthier *et al.* 2005). Some results have been obtained from this approach. In particular, they show that stratification effects are significant at all unstable wavenumbers and compressibility effects are more important at large wavenumbers. Examples of dispersion curves for various values of the stratification and the adiabatic indices are given in figure 1. The influence of the transport coefficients has also been investigated. Let us recall that a cut-off wavenumber requires the presence of a surface tension (Chandrasekhar 1961) or a diffusion of species (Duff *et al.* 1962). Indeed, viscosity without diffusion of species simply dampens the high wavenumbers. The influence of thermal conduction is roughly the same. It is also worth noting that the cut-off wavenumber may be increased by increasing either the Reynolds, Schmidt or Prandtl number.

## 5. Turbulent regime

### (a) Influence of initial conditions in isotropic compressible turbulence

The influence of initial conditions on turbulent RT mixing layers is a major concern and has been investigated by several authors, at least for weakly compressible flows (e.g. Ramaprabhu *et al.* 2005). On the other hand, the homogeneous and isotropic compressible turbulence has been studied within the framework of a statistical approach (Monin & Yaglom 1962, §20, p. 317). The *k*^{−3/2}-spectrum of acoustic turbulence has been given by Zakharov & Sagdeev (1970) and the acoustic-vortical turbulence by Moiseev *et al.* (1977). For the extension to anisotropic turbulence and a recent review, see Kuznetsov & Krasnoselskikh (2008). The numerical simulation approach has been followed by several authors (Passot & Pouquet 1987; Kida & Orszag 1992; Pais 1997; Porter *et al.* 1999; Samtaney *et al.* 2001; Xinliang *et al.* 2002; Pirozzoli & Grasso 2004) and investigations with the non-extensive statistical mechanics approach have been carried out by Shivamoggi (2003). Here, we would like to point out Chandrasekhar’s (1951) observation on the fluctuations of density in isotropic viscous turbulence in the fully compressible regime. The starting point is to define the correlation of density fluctuations *δρ* as5.1where the mean density is constant and uniform and the prime stands for the quantity taken at the point (*i*=1,2,3). From the continuity equation (2.1), one derives the evolution equation for the correlation of density fluctuations5.2In homogeneous and isotropic turbulence, it may be proved that the two correlations introduced in equation (5.2) are written as5.3

If the defining scalar *L*(*r*,*t*) of the previous correlation goes to zero faster than *r*^{−3}, it follows from equation (5.2) that the integral5.4is a constant. The quantity *I* is thus an invariant of the flow. Let *Π*(**k**) and *Λ*_{j}(**k**) be the Fourier transform of the two correlations *ϖ*(**r**) and . If the expansion of *Π*(**k**) is of the form5.5it follows from equation (5.2), written in the Fourier space, that *Π*_{o}=*c*^{st}. This means that the behaviour of the spectrum of for is given by 4*πΠ*_{o}*k*^{2}, where *Π*_{o} is a constant independent of time (Chandrasekhar 1951). In other words, the largest structures in the density fluctuations are determined by the initial conditions and represent permanent features of the flow. Chandrashekar makes an analogy with the Loitsiansky invariant. Indeed, Loitsiansky has shown, within the framework of incompressible isotropic turbulence, that the quantity , where *B*_{LL} is the double longitudinal velocity correlation, is constant in time, provided that the triple velocity correlation decreases faster than *r*^{−4} as *r* tends to the infinity. This point has been the object of controversy over the last decades. There is recent numerical evidence (Ishida *et al.* 2006) for the suppression of long-range velocity correlations in isotropic turbulence, by a screening effect owing to vorticity structure; so, the Loitsiansky integral is indeed an invariant of the turbulent flow. On the other hand, in a compressible flow, the pressure is a local quantity as opposed to incompressible flows; so, the behaviour of the correlation (5.3) deserves a careful analysis, which, to the authors’ knowledge, has not been carried out. Therefore, the conclusion regarding the behaviour of the largest scales of the fluctuations in density in isotropic compressible turbulence has to be taken with some caution. However, would Chandrashekar’s observation turn out to be true, it would reveal a specific behaviour of RT turbulent mixing layers in the fully compressible regime.

### (b) Numerical simulations

#### (i) The direct numerical simulation approach

Most of the available RT mixing layer numerical simulations have been carried out with weak stratification parameter values, and either the Euler or the NS equations are solved. Moreover, the adiabatic indices of the two fluids are usually the same and equal to *γ*=5/3 or 7/5. In these cases, both static and dynamic compressibility effects are expected to be small (Livescu 2004; Jin *et al.* 2005; Lafay *et al.* 2007). Table 2 gives some references on these numerical simulations. A more complete review of RT numerical simulations for both incompressible and compressible flows may be found in Le Creurer & Gauthier (2008). Sin’kova *et al.* (1999) performed simulations using the numerical code TREK in the limit where the stratification parameter is approximately equal to 0.2×10^{−2} with adiabatic indices equal to 7/5. Belotserkovskii & Oparin (1999) also carried out RT simulations using the equations for compressible fluids in the incompressible limit (Inogamov & Oparin 1999). In particular, they showed that the two-dimensional perturbation is unstable, i.e. the flow acquires a three-dimensional character.

Jin *et al.* (2005) provided a study of compressibility effects through single mode two-dimensional numerical simulations. The stratification parameter, , where *λ* is the perturbation wavelength, close to the parameter *Sr* used in equations (2.1)–(2.5), ranges from some 10^{−2} to 0.50 and the adiabatic indices of the two fluids differ significantly (1.1≤*γ*_{H,L}≤4). They found that pressure differences and interface shape, however, display significant dependence on the EOS parameters even for the weakly compressible flows. For three-dimensional multimode mixing, we expect accordingly that density stratification rather than drag will provide the leading compressibility effect. Such results are in agreement with linear theory results (Livescu 2004; Lafay *et al.* 2007). However, regarding three-dimensional mixing layers, large differences between the adiabatic indices of the heavy and light fluids lead, within the linear theory, to large differences in growth rates. Lafay *et al.* (2007) also study compressibility effects through single mode two-dimensional numerical simulations. They show that the conclusions drawn in the linear regime usually hold in the nonlinear single mode regime. George & Glimm (2005) noticed that, in mixing zones arising from stratified equilibrium states, thickening self-similarity is lost as two length scales, *L*_{H,L}, appear. This loss of self-similarity is associated with negative time-dependent Atwood numbers, *A*(*t*), i.e. with density reversal. The turbulent mixing thickness is then calculated by integration: *ds* *ds*_{1}. As a result, self-similarity is recovered. However, they find a strong increase in mixing rates *α*_{eff} when the stratification increases. This effect is measured through the dimensionless parameter , where *λ* is an average initial wavelength and simulations are carried out with stratification parameter ranging from *M*^{2}=0.008 to *M*^{2}=1.

#### (ii) The large-eddy simulation approach

Direct numerical simulations (DNSs) of turbulent flows are severely limited in Reynolds numbers, even with the availability of large parallel computers and a fortiori for compressible flows. In RT flows, this means that the self-similar regime may be difficult to attain (Mattner *et al.* 2004). This is the motivation for investigators to use the large-eddy simulation (LES) approach, in which large scales are simulated and small scales are modelled. There are numerous LES-type methods, from the simplest one, the MILES (monotone-integrated large-eddy simulation) method where no subgrid-scale model is added to the equations, to more sophisticated methods (Meneveau & Katz 2000). The idea behind the MILES approach is that the numerical dissipation resulting from monotonic numerical schemes is close to the dissipation of the subgrid-scale model. A quantitative comparison between MILES simulations based on various numerical schemes and coming from different institutions has been reported in Dimonte *et al.* (2004). These simulations have been carried out in a weakly compressible limit, i.e. the stratification parameter is small and the adiabatic indices are equal to 5/3. The numerical results are in reasonable agreement with previous DNSs, experiments and theory, except for the bubble amplitude, which is smaller by a factor of 2 than the value found in the experiments.

Mattner *et al.* (2004) report on LES RT simulations by using the stretched-vortex Lundgren subgrid-scale model and mixing models (Pullin 2001). Moreover a LES–DNS comparison is carried out, where the spatial resolution is divided by four in each direction for the LES calculation, and gives similar averaged results. Mellado *et al.* (2005) also carried out LES of RT-mixing layers by using a dynamic mixed model. They argue that the turbulent Mach number (defined as , where *K* is the turbulent kinetic energy and 〈*c*〉 the average speed of sound) is limited in such mixing layers, and for the particular case of perfect gases, the bound is between the subsonic values 0.25 and 0.60. Consequently, the authors conclude that compressibility effects in RTI-induced flows are small. This conclusion is partly supported by the weakness of the acoustic part of the density fluctuations. Indeed, decomposition of the density fluctuations into the entropy part and the acoustic part shows that the latter is less than 10 per cent. However, this conclusion is in contradiction with previous findings in the linear (Livescu 2004; Lafay *et al.* 2007) and nonlinear (Jin *et al.* 2005; Lafay *et al.* 2007) regimes, where both effects of stratification and EOS effects (*γ*’s variation for perfect gases) may be large. Recently, Olson & Cook (2007) gave an example of a strong compressible behaviour in the turbulent regime. Indeed, they show that an RTI may launch a shock wave into the upper fluid provided that its height is large enough. In this numerical experiment, rising bubbles act as pistons, which generate shocklets. These shocklets coalesce into a strong shock wave, which increases in strength as it propagates. They also show that the turbulent Mach number, *M*_{t}(*z*,*t*), reaches maximum supersonic values just behind the shock wave at large times, although *M*_{t}(0,*t*) at the initial location of the mixing layer remains small. These simulations have been carried out with a specific LES strategy. Indeed, Cook (2007) suggests to modify transport coefficients (shear and bulk viscosity, thermal conductivity and species diffusivity) to damp out high wavenumber modes instead of filtering the governing equations. On the other hand, let us note that a two-point closure model for isotropic compressible turbulence has been devised (Fauchet & Bertoglio 1999), which may be used as a subgrid-turbulence model.

### (c) Compressible turbulence modelling

Large-scale three-dimensional numerical simulations of turbulent flows produce a large amount of data that has to be analysed. One way to reduce these data is to average in the homogeneous transverse directions in order to obtain one- or two-dimensional profiles. For RT flows, both transverse directions are homogeneous so that the averaging process leads to one-dimensional profiles. By doing so, one gets close to the philosophy of turbulence modelling in which an ensemble-average process is first carried out. The second step, the closure, or modelling step, consists of calculating correlations between fluctuating variables of order (*n*+1) versus correlations of order *n* through phenomenological algebraic or differential equations balanced by numerical constants. This method produces engineering models used to quickly calculate turbulent flows. It also leads to a worthy viewpoint on the turbulent flow that may further our understanding. In compressible flows, one classically uses the mass weighted averaging, or Favre averaging, defined as . Such averaging provides equations with a form similar to the instantaneous equations, only with additional terms, and very similar to the ones obtained from a Reynolds-averaging process carried out on the equations for incompressible flows. Modelling is then usually performed with the same rules used for incompressible flows. This is a great advantage in terms of simplicity of modelling. However, only the variable density effects are taken into account, and these in a simple way. Note that the Reynolds-averaging is also justified for compressible flows but produces much more complicated average equations. A comparison between the two approaches has been carried out by several authors (see the example of the isothermal turbulent binary mixture flow, Ould Sid Ahmed & Lili 2002). Compressibility effects are usually much more complex than a simple ‘variable density’ as acoustics and thermodynamics are involved (see the Kovásznay modes in §2*b*). The objective of this section is not to review the various turbulence models developed for compressible flows. The reader is referred, for example, to the book of Chassaing *et al.* (2002, ch. 10). Instead, we would like to show first that the knowledge accumulated over the last two decades on turbulent compressible Kelvin–Helmholtz mixing layers may be beneficial to the investigation of RT mixing layers (Chassaing *et al.* 2002). Second, we detail some evolution equations for typical compressible correlation, without referring to closure assumptions. Turbulent shear mixing layer experiments have shown a strong growth rate reduction for increasing convective Mach numbers. For a shear mixing layer, the freestreams of which have velocities *U*_{1} and *U*_{2}, the convective Mach number is defined by5.6

This reduction occurs at subsonic convective Mach numbers, *M*_{c}≈0.5, and increases for larger *M*_{c}. At supersonic Mach numbers, the growth rate of the thickening is only 40 per cent, the value of the incompressible growth rate. First, the linear stability theory (Sandham & Reynolds 1989) exhibits this strong reduction of growth rate. The physical processes described by the linear theory are also observed in the later large-scale nonlinear development. As a result, linear growth rate correctly predicts the experimentally observed trends in thickening rate that are due to velocity ratio, density ratio and Mach numbers. Second, the turbulent stage has been investigated through experiments, numerical simulations and modelling. In particular, the modelling of these mixing layers has revealed that both the stratification and finite speed of sound effects had to be taken into account to correctly represent the growth rate reduction for increasing convective Mach numbers. Regarding modelling, the first step is to recognize the specific compressible contributions in the average equations. There are essentially three terms in the Favre-average equations, which may contain some compressibility effects. The first compressible contribution is the pressure-dilatation correlation , which arises in the exact average equation for the internal energy and the turbulent kinetic energy. This term is interpreted both as a source term and as a contribution to the dissipation. The underlying physical phenomena could be the presence of shocklets in the mixing layer such as those identified in isotropic compressible turbulence by Passot & Pouquet (1987). The second contribution is the dissipation term , where , which also arises in the same equations. This term contains a hidden compressible part, which has been exhibited (Zeman 1990; Chassaing *et al.* 2002, ch. 10). Indeed, the dissipation may be expressed as5.7where , and is the non-homogeneous contribution. The compressible part of the dissipation has been isolated in the fluctuations of the divergence of the velocity. These contributions have received specific modelling in shear flow studies. The important point here is that these modellings are not sufficient to explain the growth rate reduction of the turbulent Kelvin–Helmholtz mixing layers for increasing convective Mach numbers, in configurations with high density ratios. Indeed, it has been argued that the stratification has to be taken into account. This may be achieved by considering variable phenomenological constants in the closure relations (Huang *et al.* 1994; Guézengar *et al.* 2000).

In some situations, for example, in RT and Richtmyer–Meshkov mixing layers and also in shear flows, the correlation of the mass flux times the mean pressure gradient, i.e. , is the dominant source term. It is then preferable to determine the turbulent mass flux evolution through a differential equation (Grégoire *et al.* 2005) instead of an algebraic relation. The exact evolution equation is obtained from the NS equations and the Favre decomposition. This equation is somewhat difficult to interpret and one may use the expansion of the density5.8to neglect triple and higher-order correlations. The resulting approximate evolution equation for the turbulent mass flux then reads5.9The three first terms of this equation are source terms owing to gradients of mean density, velocity and pressure, respectively. The following three terms are of diffusion type, whereas the three last ones represent the dissipation. The exact equation for the root mean square (r.m.s.) value of the density fluctuations reads5.10In this equation, the first term of the right-hand side (r.h.s.) represents the diffusion of the r.m.s. value of the density fluctuations, the second and third terms are mean-flow production terms, whereas the two last ones are the dissipation.

We note that Ahmadi and co-workers (see Ahmadi 2004) have been developing turbulence models for compressible flows based on the second law of thermodynamics. This method provides models that resemble classical models. Some source terms of the dissipation rate equation are, however, different and there are, some, constraints on the phenomenological coefficient values.

Over the last decades, the study of RT-compressible turbulent mixing layers has lagged behind the general evolution of compressible turbulent flow investigations. These turbulent flows were studied and calculated with closure turbulence models, in the first instance with first- and then with second-order type. Their complexity has increased together with the number of phenomenological constants. Modelling requires experimental or numerical data to suggest and adjust closure relations and to calibrate numerical constants. The second development in turbulent mixing layer investigations was naturally numerical simulations, but true DNSs require highly accurate schemes (Jameson 2000), as numerical dissipative and dispersive errors act as a crude subgrid-scale model. This difficulty arises in incompressible flow simulations as well, but numerical methods for incompressible NS equations are more mature than the ones for compressible flows. An effort to design high-order numerical schemes for subsonic and supersonic flows is currently made by several investigators and should become more and more popular as we are entering the petascale era. This approach provides reference data from which turbulence structure and mixing structure may be investigated. Turbulence modelling may also benefit from these data. DNS also provides data that cannot be obtained from experiments although it is limited to weak and intermediate Reynolds number flows. By contrast, very large Reynolds numbers may be reached with the LES approach, as small scales are modelled. However, as the mixing process involves small scales, the picture of the mixing obtained from LES simulations partly reflects the subgrid-scale model properties. As a result, DNS and LES approaches only give a partial picture as they have their own limitations and they both contribute to the understanding of the basic mechanisms in turbulent mixing. Turbulence modelling has essentially become an engineering tool.

## 6. Summary and concluding remarks

We have reviewed the compressibility effects in the RTI by stressing the two different roots of these effects. The first one, called ‘static compressibility’, is the stratification of the initial hydrostatic equilibrium state. The second one, called ‘dynamic compressibility’, is essentially an effect of the finite speed of sound. On the other hand, small amplitude motions in a compressible flow may be decomposed into three modes, i.e. vorticity, acoustic and entropy (Chu & Kovásznay 1957). In some situations, the acoustic mode may be neglected, but the stratification may not be negligible. The quasi-incompressible limits of the NS equations, the low Mach number (for unstratified equilibrium states) and the anelastic (for stratified equilibrium states) approximations have been developed, thus removing the acoustic waves and the associated restrictive CFL condition in numerical simulations.

The known results regarding the stability of stratified flows have been recalled (Gamaly *et al.* 1976; Sitt 1980; Scannapieco 1981). Two types of modes may be identified: convective modes or gravity modes and acoustic-type modes. They have been exhibited in two examples one in isentropic stratified atmospheres (Scannapieco 1981) and the other one in imploding spherical shells (Book 1978).

The linear regime of the compressible RTI has also been reviewed. The various linear problems of RTI for perfect compressible fluids have been listed according to the thermodynamic state (isentropic, isothermal and general) of the equilibrium state and the perturbations, respectively. The configurations previously solved and the ones that have to be solved appear clearly. It thus opens the door for further contributions. Stability analysis of perfect fluids is usually handled with analytical methods (Livescu 2004) whereas, for viscous fluids, numerical methods are required (Lafay *et al.* 2007). Generally, stability analysis of RTI always led to a couple of real eigenvalues, which cross the imaginary axis. The associated behaviour is a controlled departure from incompressible type behaviour owing to stratification and dynamic compressibility. This departure is essentially localized at small wavenumbers and at small Atwood numbers. On the other hand, some unsuccessful attempts have been made in the past to find more complex behaviour, such as oscillatory motion. These attempts should be pursued by systematic investigations of the spectrum of the linear RT operator of the NS equations.

Regarding the numerical simulation of the turbulent mixing regime, Chandrasekhar (1951) pointed out that, in the case of isotropic compressible turbulence, the largest structures in the density fluctuations are determined by the initial conditions and represent permanent features of the flow. However, this remark still has to be checked through numerical simulations. If it were to turn out to be true, it would reveal a specific behaviour of RT mixing layers in the compressible regime. Some insights have been made in compressibility effects in three-dimensional mixing layers. It appears that stratification is the most important effect with respect to variations of the EOS parameters (Jin *et al.* 2005). In particular, self-similarity, versus the variable *A* *g* *t*^{2}, is lost in the compressible regime, but may be recovered if one introduces a local time-dependent Atwood number. A strong increase of mixing rates is also found for strongly stratified configurations (George & Glimm 2005). These findings have to be confirmed. An example of typical compressible behaviour has been exhibited where shocklets generated by RT bubbles coalesce into a shock wave (Olson & Cook 2007). More generally, the analysis in terms of average equations is useful for the analysis of the data obtained from large-scale numerical simulations. This approach may also lead to specific statistical models for the fully turbulent regime, which allows the rapid calculation of turbulent mixing layer evolution. In that respect, the knowledge of turbulent compressible Kelvin–Helmholtz mixing layers may be beneficial to the investigation of RT mixing layers.

## Acknowledgements

One of us (S.G.) thanks J.-M. Clarisse (CEA/Bruyères le Châtel, France) for helpful discussions over the last years. We also thank M. R. Connolly (Cavendish Laboratory, Cambridge, UK) for his help in the final step of this work.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond’.

- © 2010 The Royal Society