A tentative review is presented of various approaches for numerical simulations of two-fluid gaseous mixtures at high density ratios, as they have been applied to the Rayleigh–Taylor instability (RTI). Systems exhibiting such RTI behaviour extend from atomistic sizes to scales where the continuum approximation becomes valid. Each level of description can fit into a hierarchy of theoretical models and the governing equations appropriate for each model, with their assumptions, are presented. In particular, because the compressible to incompressible limit of the Navier–Stokes equations is not unique and understanding compressibility effects in the RTI critically depends on having the appropriate basis for comparison, two relevant incompressible limits are presented. One of these limits has not been considered before. Recent results from RTI simulations, spanning the levels of description presented, are reviewed in connection to the material mixing problem. Owing to the computational limitations, most in-depth RTI results have been obtained for the incompressible case. Two such results, concerning the asymmetry of the mixing and small-scale anisotropy anomaly, as well as the possibility of a mixing transition in the RTI, are surveyed. New lines for further investigation are suggested and it is hoped that bringing together such diverse levels of description may provide new ideas and increased motivation for studying such flows.
Molecular mixing in response to stirring by turbulence is an important process in many practical applications. When the microscopic densities of the fluids participating in the mixing are very different, these flows have been referred to as variable density (VD) flows in contrast to the Boussinesq approximation in which the densities are commensurate . In general, turbulent density fluctuations have been studied in conjunction with compressibility effects (e.g. in aeronautics) or as a result of temperature changes (e.g. in combustion). Such flows have been the subject of numerous fundamental turbulence studies, and many modelling strategies now exist . Even though all these effects can be generically described as ‘VD’, here, this term is reserved for density variations arising owing to mixing between fluids with vastly different molecular weights. In this case, fundamental turbulence studies as well as specific engineering models are scarce .
In VD flows, even if the fluids participating in the mixing are incompressible, then the density and specific volume change with the mixture composition, and the velocity field is not divergence free. VD mixing occurs in atmospheric and oceanic flows, astrophysical flows, combustion and many flows of chemical engineering interest. Many of these flows are driven by acceleration (e.g. gravity in geophysical and astrophysical flows) which, because the density is not uniform, leads to large differential fluid accelerations. If the acceleration is constant and the fluid configuration is unstable (i.e. the density gradient points opposite to the acceleration), then a fluid instability is generated in which small perturbations of the initial interface between the two fluids grow, interact nonlinearly and lead to turbulence. This instability is known as the Rayleigh–Taylor instability (RTI) and is of fundamental importance in a multitude of applications, from fluidized beds, oceans and atmosphere, to inertial, magnetic or gravitational confinement fusion, and to astrophysics.
For example, RTI has a profound effect in the descent of heavier, saltier waters in the North Atlantic Ocean which relates to the thermohaline circulation. Without this, the earth climate would change dramatically . In astrophysics, RTI plays a crucial role in flame propagation and the development of ignition bubbles for detonation in type Ia supernovae . These are the largest explosions in the modern universe and have been used as calibratable standard candles by cosmologists (including the 2011 recipients of the Nobel prize). RTI also determines the time-scale of mixing and burning and hence the rise time of burst in helium-ignited type I X-ray bursts. Thus, RTI is the dominant acceleration mechanism for the flame. RTI can be found in such diverse areas as sound generation by snapping shrimp owing to the collapse of cavitation bubbles , sonoluminescence  or industrial coating with thin liquid films .
Although this instability has been subjected to intense research over the last 50 years (e.g. over 100 RTI-related papers are published per year in peer-reviewed journals ), until recently numerical studies have been restricted to coarse mesh calculations of the Euler equations. On the other hand, it is notoriously difficult, in laboratory experiments, to accurately characterize the initial conditions and provide the detailed measurements required for turbulence model development and validation. Thus, a large number of open questions remain unanswered about this instability and even first-order global quantities, such as the layer growth, are not completely understood and still give rise to intense debate . Although some interesting results have recently emerged on the flow properties [11–20] much more needs to be done, and the material mixing problem remains largely unexplored, especially at high density ratios (see [15,20–25] for reviews and open problems). Nevertheless, today’s petascale computers allow fully resolved simulations of RTI at parameter ranges comparable with those attained in laboratory experiments, but providing, in carefully controlled initial and boundary conditions studies, much more information than physical experiments. These extremely high-resolution simulations are enabling a look at the physics of turbulence and testing turbulence models in unprecedented detail, hopefully contributing to a significant advance in our understanding of RTI.
A review of the simulations of material mixing in the context of RTI would take many pages, owing to the vastness of the subject. For example, materials participating in the RTI in practically relevant situations may have elasto-plastic [26,27] or thixotropic  non-Newtonian behaviour, may have complicated equations of state and react thermonuclearly , have sharp material interfaces, undergo phase changes, etc. While all these problems are fascinating, here, the attention is restricted to a case that can be successfully approached both numerically and experimentally: buoyancy-driven mixing between two inert gases at large molecular weight ratios. In this case, there are no sharp interfaces. These pose serious problems even neglecting molecular interfacial aspects, at least in a macroscopic treatment, owing to complicated topologies, with interface break-up and reconnection. Diffuse treatments of interfaces, such as those based on Korteweg stresses, are not considered here, since they do not address the miscible case (note that surface tension does not appear in inert gaseous mixtures).
The primary non-dimensional parameter characterizing the differential acceleration effects is the Atwood number: 1.1where ρ1 and ρ2 are the densities of the two fluids, with ρ2>ρ1. The Atwood number ranges from 0 to 1. For air interpenetrating helium A≈0.75; for air and hydrogen A≈0.87. Thus, the high A problem is important, in practice, yet most studies to date address the Boussinesq approximation, when .
However, even after narrowing the field to inert gaseous mixtures, the range of RTI applications still spans an immense extent of scales, from small atomistic systems, where molecular dynamics (MD) and direct simulation Monte Carlo (DSMC) become the relevant numerical tools , to macroscopic systems, where the continuum approximation applies. These areas are covered here for two reasons. First is to provide a comprehensive view of the equations considered in various studies, with their limitations and range of applicability. Second, on today’s computers, accurate simulations of RTI at relevant parameter ranges are becoming feasible for all these systems. However, the next generation of supercomputers may pose significant challenges for the present numerical treatments and cross-breeding ideas can be proved useful. It is hoped that this review will provide some cross-fertilization among several fields where RTI is important and provide an increased motivation for furthering the understanding of the physics of high A RTI mixing. Part of this material has evolved from discussions and insights gained during the Second International Conference Turbulence Mixing and Beyond .
The paper is organized as follows. Section 2 presents the governing equations used in various studies of two-fluid mixing in gaseous mixtures with applications to the RTI, starting from MD and ending with the Navier–Stokes (NS) equations. Both compressible and incompressible NS equations are considered. Because the compressible to incompressible limit is not unique and it is important to define the correct limit as the basis of comparison for compressibility effects studies, a discussion on two relevant limits for RTI has been added in §2c(ii). Section 3 surveys the numerical simulations for the equations presented. There are myriad simulations of macroscopic RTI on coarse meshes. Because high-resolution, fully resolved simulations are becoming increasingly more common on today’s computers, these coarse mesh simulations are covered only when new or unexpected physics are revealed and relevant previous reviews are mentioned for the simulations omitted. Finally, a summary and concluding remarks are provided in §4.
2. Equations describing material mixing
Here, various frameworks of approximation for describing material mixing in fluid flows, as they have been used in the study of the RTI, are discussed. Each of these frameworks fits into a hierarchy of theoretical models describing the material behaviour at a certain level. Often, in passing to the next (coarser) level, the description becomes unclosed, and closure relations are required. Thus, for each theoretical model presented, its assumptions and level of description will be stated.
The most fundamental level of description of a material consisting in N interacting particles having positions xj and velocities vj, with , is by the time-dependent Schrödinger equation (TDSE): 2.1where is the imaginary unit, is the reduced Plank constant, ψ(x1,…,xN,t) is the wave function for the N particles (the probability amplitude for different configurations of the system at different times), Δxj is the Laplace operator at position xj, mj is the mass of the jth particle and V (x1,…,xN) is the interaction potential. The integration of equation (2.1) for more than a few particles is computationally prohibitive and its application so far to systems of particles has been limited, in an approximate form, to the description of very small collections of atoms or molecules .
(a) Atomistic descriptions
The size of the quantum mechanical wave packet representing each particle compared with the average distance between particles, R=n1/3 with n the number density of particles, can be characterized by the ratio : 2.2where Λ is the thermal de Broglie wavelength, kB is the Boltzmann constant and T is the temperature. When this ratio is reasonably small (e.g. for molecular nitrogen at atmospheric conditions but not for molecular hydrogen), the quantum effects can be neglected in describing the atomic nuclei and the classical limit, in which the random variables xj,vj can be replaced by their means, becomes valid . There are additional requirements for a full classical treatment with respect to the rotational temperatures in the case of molecules and electron interactions at high temperatures and/or pressures; however, these are not addressed here. When the classical treatment can be used, the equations describing the particle motions become 2.3and 2.4where and V (x1,…,xN) is the interaction potential, also known as the force field. These are the equations solved in MD simulations . In general, the potentials used in MD simulations are empirical and neglect the effects associated with the atomic structure. In addition, to decrease computational requirements, large-scale simulations use simple two-body potentials and also neglect polarization effects.
Even though the recent increase in computer power has led to substantial progress in the size of the MD simulations, the number of molecules in the largest MD simulations to date  is still 13 orders of magnitude smaller than the number of molecules in 1 cubic metre of air at atmospheric conditions, which is approximately 3×1025. Thus, MD RTI simulations on today’s computers describe very small systems, in which departures from the continuum approximation are expected. The results of such simulations are discussed in §3a(i).
(b) Boltzmann equation
Instead of seeking to describe the position and velocity of each particle, as is the case with MD, the Boltzmann equation describes the statistical distribution of one particle in a fluid : 2.5where x, p and m are the position, momentum and mass of the particle, respectively, F is an external force field and f(x,p,t) is the probability density function (PDF) of one-particle phase space (the number of particles having the position and momentum within a volume, d3x, and momentum–space, d3x, elements about x and p, respectively). The term on the right-hand side of the equation models the collisions between particles and needs to be specified.
Equation (2.5) is one of the most important equations of non-equilibrium statistical mechanics. One powerful method for solving this equation for finite Knudsen number, Kn, fluid flows uses a random Monte Carlo approach. The Knudsen number is the ratio of the molecular mean free path to a characteristic length scale of the flow. This method, known as DSMC , uses simulation molecules that represent a large number of real molecules in a probabilistic simulation to solve the Boltzmann equation. Intermolecular collisions and molecule–surface collisions are calculated using probabilistic, phenomenological models. The fundamental assumption of the DSMC method is that the molecular movement and collision phases can be decoupled over time periods that are smaller than the mean collision time. Compared with MD, where the time step is restricted by the collision time, the decoupling between the particle movement and collision steps allows much larger time steps in DSMC. Classical DSMC assumes a collision process that leads to an ideal equation of state. DSMC has traditionally been applied to rarefied gas flows, although recently it has also been used to model near continuum flows (Kn<1) and even canonical fluid dynamics problems. Applications of DSMC to the RTI are discussed in §3a(i). A review of the method with applications and limitations can be found in reference  and a discussion on the statistical bias due to finite sampling in the presence of thermal fluctuations can be found in reference .
(i) Lattice Boltzmann equation
Direct solving of the Boltzmann equation requires a large number of particles and has severe computational limitations in terms of applicability to fluid flows in regimes where the continuum approximation is valid. A different direction for solving this equation, in a minimal form, has seen a tremendous growth over the past 20 years and constitutes an intense area of research. In this approach, the particle movement is restricted to a discrete time–space lattice, so that the particles can move only alongside a number of finite directions [41,42]. With this assumption, in the absence of external forces, the Boltzmann equation becomes the so-called Lattice Boltzmann equation [43,44]: 2.6and, the associated solution method is known as the lattice Boltzmann method (LBM). In equation (2.6), fi(x,t) is the probability of finding a particle at time t and position x with velocity vi=ci. Here, defines a set of discrete speeds, ci, connecting the nodes of a regular lattice. A typical choice in three dimensions is a lattice formed by the midpoints of the edges and faces of a cube, resulting in 19 points. The collision term in the Boltzmann equation is modelled in equation (2.6) through a relaxation to a local equilibrium which is controlled by the scattering matrix Ωij, whose eigenvalues define the time-scale for local equilibration of the various kinetic moments. In practical applications, the scattering matrix is often taken in its simplest diagonal form, resulting in the so-called lattice Bhatnagar–Gross–Krook (LBGK) scheme, after the continuum model Boltzmann equation introduced by Bhatnagar et al. . The local equilibrium is the lattice analogue of a local Maxwellian with density and flow velocity . Usually, values are calculated from a second-order expansion in Mach number of the local Maxwellian, which can be represented on a 19-point lattice in three dimensions . The second-order expansion is necessary to ensure the correct behaviour for low-Mach numbers, quasi-incompressible flows. Higher-order expansions are required to capture stronger compressibility effects. Another requirement on the expression is that it satisfies a discrete form of Boltzmann’s H-theorem  which ensures unconditional stability of the method. In LBM, the fluid properties are simply found by summing the moments of fi over the lattice points. It should be noted that the thermal fluctuations, accounted for in the DSMC, are no longer represented in the LBE.
The earlier forms of LBE were restricted to isothermal low-Mach numbers flows with ideal fluids. The BGK model is also restrictive as it implies the same diffusivity for mass, momentum and energy. Nevertheless, this is an active area of research and models with multiple relaxation times and non-ideal fluid behaviour have been proposed or under development to address complex material behaviour. In §3a(ii), recent applications to the RTI problem and developments for two-fluid mixing are discussed.
(c) Continuum descriptions
(i) Navier–Stokes equations
The multi-component equations of motion for gaseous mixtures in the continuum limit can be derived from the Boltzmann equation following the Chapman–Enskog theory . The derivation assumes near equilibrium systems, such that f(x,p,t) can be expanded about the Maxwellian distribution. To first order in the Knudsen number ξ≡cr/(Lνr), where cr is the reference flow speed, L is the reference length required for a significant change in f and νr is the reference collision frequency, the expansion gives rise to the NS equations. These equations have been found to be applicable for considerably large deviations from equilibrium.
In vector form, the governing equations describing the conservation of mass, momentum, energy and species mass fractions, for inert, miscible materials are : 2.7 2.8 2.9 and 2.10The primary dependent variables in equations (2.7)–(2.10) are the density, ρ, mass averaged velocity, u, specific internal energy, e and species mass fractions, Y α, where and . The external body force, Fα, acting on species α is considered specified (not derived). Note that ρ and u are mixture values. Each species has its own velocity, which differ from the mixture velocity by the diffusional velocity, Vα.
The molecular transport terms on the right-hand sides of equations (2.7)–(2.10) represent transport of mass (the diffusional velocities, Vα), momentum (the stress tensor, σ) and energy (the heat flux, q). These quantities cannot, in general, be related to the primary variables, because the transport properties involve higher moments of the velocity distribution function. Neglecting radiative effects and using second-order expansions in Sonine polynomials to solve the linear integral equations for the components of the first-order term in the f expansion, the Chapman–Enskog theory leads to the following expressions for the transport terms [49,50]: 2.11 2.12 and 2.13where Xα is the mole fraction of species α (), Wα is the molar mass of species α, p is the pressure of the mixture, hα is the enthalpy of species α, is the universal gas constant, T is the temperature of the mixture, I is the unit second-order tensor and (∇u)T denotes the matrix transpose of the velocity gradient.
The coefficients appearing in the expressions for the molecular transport terms are the diffusion coefficient, Dαβ, of species α and β, thermal diffusion coefficient, DTα, of species α in the mixture, the coefficients of viscosity, μ, and bulk viscosity, μB and conduction coefficient, λ. These coefficients depend on the collision term in the Boltzmann equation and they have been calculated for various expressions of this term . The coefficients of viscosity and diffusion, which are associated with the transfer of mass and momentum by the translational motions of the molecules, are successfully predicted for many poly-atomic gases by the mono-atomic theory. Semi-empirical formulae (including for mixtures) have also been published . The conduction coefficient and bulk viscosity calculations from the kinetic theory are more problematic. For example, for mono-atomic gases, the kinetic theory predicts μb=0, which is the value used in most simulations discussed here. Nevertheless, in nitrogen, for example, experimental measurements found values of μB/μ≈0.8 at moderate temperatures, so that the effects of μB cannot be neglected. In addition, recent simulations  indicate that even the classical RTI with compressible fluids develops shock waves, so that the divergence of velocity is not, in general, small. Similarly, for the conduction coefficient, so far, experimental data are still preferred over kinetic theory calculations.
Equation (2.11) for the diffusion velocities is, in general, implicit and requires the solution of a linear system of equations. For binary systems, the formula becomes explicit: 2.14The first term on the right-hand side of equation (2.14) defines the familiar Fick’s law of diffusion: 2.15if all other effects are neglected. This is the most used formula for practical binary mixing calculations. However, for the compressible RTI, the pressure gradients may not be small. Thus, if the two fluids have large molecular weight ratios, the baro-diffusion term should be retained in equation (2.14). The last term in equation (2.14) defines the thermal-diffusion or Soret effect. It was not known before the Chapman–Enskog theory and, since its experimental confirmation, represents one of the significant accomplishments of the theory. The Soret effect is usually neglected in RTI calculations; however, this term may be non-negligible in compressible two-fluid mixing problems at large molecular weight ratios . Therefore, the importance of the Soret effect, especially in the presence of the shock waves generated by the compressible RTI, needs to be assessed in future simulations.
The first term on the right-hand side of equation (2.13) defines the familiar Fourier’s law of conduction: 2.16For multi-component mixing, this is clearly inadequate, as the enthalpy diffusion term is, in general, non-zero and important for entropy conservation . In addition, the derivation described in Landau & Lifshitz  can be extended to the multi-component case to show that the second law also requires the inclusion of the baro-diffusion term in the diffusion velocities. This is important, as multi-component calculations still use Fickian diffusion even when the enthralpy diffusion is accounted for and, thus, do not satisfy the second law. The thermal diffusion term, also known as Dufour effect, is small, even when the Soret effect in the diffusion expression is important [48,52].
Several non-dimensional numbers can be introduced to define the relative importance of the molecular transport terms in equations (2.7)–(2.10). Thus, the Prandtl number, Pr=μcp/λ, where cp is the specific heat at constant pressure, is a measure of the relative importance of momentum and heat transport; the Schmidt number, Sc=μ/ρD12, is a measure of the relative importance of momentum and mass transport. The relative importance of the advective and molecular transport terms in momentum, energy and species equations can be inferred from the values of Re, RePr and ReSc, respectively. The Reynolds number, Re=ρUL/μ, where U is the characteristic velocity at scale L, is a measure of the relative importance of advective and molecular transport of momentum at scale L. In turbulence calculations, assuming that the continuum approximation holds, an appropriate Reynolds number can also be defined to characterize the range of dynamically relevant scales of motion . In most RTI applications, such a Reynolds number is large; however, Pr and Sc may vary widely from very small to very large values.
Because the theory for the collision term and the calculation of the transport coefficients beyond mono-atomic gases is still under development, a practical approach is to calibrate these coefficients to match experimental measurements, using various empirical formulae . Phenomenological models can also be developed to include the effects beyond those of gaseous mixtures (e.g. strength  or thixotropy ). This shows the versatility of the conservation equations and underscores their widespread use in scientific and engineering applications. If the transport coefficients are set to zero in equations (2.11)–(2.13), then one obtains the Euler equations. These equations can be derived in a more general context than from the collisionless Boltzmann equation, for example, from many-body quantum mechanics .
The importance of the thermal fluctuations and subsequent validity of the NS equations at scale Δx can be inferred from the ratio : 2.17where S is the magnitude of the shear rate in the flow under consideration. For most laboratory flows at atmospheric conditions, including RTI, this ratio is small. Indeed, assuming a turbulent flow, let Δx=ηK, where ηK is the Kolmogorov microscale, the scale at which most viscous dissipation occurs . The Reynolds number at the Kolmogorov microscale is 1, such that S∼μ/(ρη2K). Then, the condition C≈1 can be rewritten as 2.18The same formula for ηK at which the viscous and thermal fluctuations effects are commensurate is also obtained for the Landau–Lifshitz NS equations (LLNS). These are the usual NS equations enhanced with stochastic terms in the momentum and energy equations accounting for the random thermal fluctuations [58,59]. For air at atmospheric conditions, formula (2.18) yields ηK≈3 e−11 m, a value at least nine orders of magnitude smaller than those encountered in most laboratory experiments of the RTI. Nevertheless, in §3a(i), RTI simulations on domains small enough to lead to values of ηK close to those predicted by formula (2.18) are discussed.
The thermodynamic variables are related through relations known as equations of state. Thus, by choosing the temperature and density as the primary variables, the internal energy and pressure can be calculated as the partial derivatives of the partition function for the system . If the particles in the fluid are non-interacting, then one obtains the ideal-gas equation of state: 2.19whereas the caloric equation of state is 2.20where the specific heat at constant pressure for the mixture is calculated as the mass fraction-weighted sum of the individual specific heats. Gases obeying equation (2.20) are usually called ‘perfect’; in this case, the specific heats vary with temperature owing to energy contributions beyond those contained in translational motions. The ratio between the specific heats at constant pressure, cp, and constant volume, cv, is denoted by γ. If the mixture follows equation (2.19), then , so that cp=γR/(γ−1) and cv=R/(γ−1). In general, for a perfect gas mixture, γ=γ(T,Y α). Note that , with the maximum value obtained for mono-atomic gases obeying the ideal-gas equation of state. Equation (2.19) also leads to a simple formula for the speed of sound, c2≡(∂p/∂ρ)s=γRT=γ(p/ρ). It should be noted that the closures for the diffusion velocities and heat flux given in (2.11) and (2.13) imply ideal-gas equation of state. More general formulas can be written using derivatives of the chemical potential instead of the gradient of the mole fraction and the species enthalpy; however, these are beyond the scope of this paper.
In §§3a(i) and 3a(ii), RTI simulations with compressible materials are discussed.
(ii) The incompressible limits for two-fluid mixtures
In many practical situations, RTI occurs in regimes where gases can be considered as incompressible. For single fluid flows, the incompressible limit is not unique, for example constant density flows versus flows with background stratification. Consequently, the limiting process is also not unique. Constant density flows can be obtained as (however, mathematically, this limit is not unique; see the discussion below), whereas incompressible flows with background stratification can be obtained through the anelastic approximation . The latter case assumes that the thermodynamic variables have small variations around the background state. In the limit of small stratification, the anelastic approximation reduces to the Boussinesq approximation. For the RTI, such anelastic assumption is not justified, because the background stratification varies considerably as the instability develops: initially, the two fluids are segregated in an unstable configuration, whereas the end state is fully mixed.
For the inviscid immiscible compressible RTI in the linear regime, the condition ∇⋅u=0, can be obtained mathematically as either or ; both of them also define [61,62]. The two limits are different: leads to uniform density in each of the fluid regions, whereas allows for non-constant background density. This non-uniqueness has been studied for the barotropic version (p=p(ρ)) of the NS equations in reference . If non-ideal effects are considered (e.g. viscosity, heat conduction), then does not lead to ∇⋅u=0; however, this case is still incompressible, as the hyperbolic part of the NS equations is replaced by an elliptic equation given through the velocity divergence condition. Physically, γ cannot exceed 5/3; nevertheless, for the inviscid immiscible RTI in the linear regime, the limit is the same as that obtained if , where β is the exponent in a politropic transformation the flow is assumed to undergo . In this case, the system is not closed and the energy equation should contain an appropriate source/sink such that it can be replaced by the politropic transformation. Even though, physically, γ is bounded, the mathematical considerations show that p and γ (γ1 and γ2 for the two-fluid case) are independent compressibility parameters and need to be studied separately. Moreover, an appropriate incompressible limit needs to be defined for each case, in order to be able to make meaningful comparisons. These are constructed below.
The ‘first’ incompressible limit. This limit is taken to be when the microdensities, ρ1 and ρ2, of the two materials are constant. Let ρ2>ρ1 for the sake of clarity. The case of interest here is ρ2≫ρ1. Because the mass within an arbitrary control volume is equal to the sum of the masses of the two components, it yields 2.21
For two-fluid gas mixtures obeying equation (2.19), the microscopic densities, defined by and , approach constant values as and , with the additional constraint that . In this case, the ideal-gas equation of state reduces, at the limit, to equation (2.21). Note that the constraint needs not be satisfied by the unperturbed flow in the corresponding compressible case; any unperturbed flow asymptotes to the same incompressible limit.
For the compressible RTI with the unperturbed state in thermal equilibrium, the unperturbed pressure and density are in hydrostatic equilibrium. In this case, the unperturbed pressure and density profiles in each fluid (outside the mixing layer) are stably stratified, with the stratification determined by the reference pressure [61,65]. At larger reference pressures, the initial stratification is reduced, so that the pressure can also be interpreted as a stratification parameter [65,66]. This could lead to some analogies to the rich literature on stably stratified flows; however, within the mixing layer, finite pressure values have clear compressibility effects implications, as they can lead to the formation of shock waves . In addition, the interpretation of pressure as a stratification parameter makes sense only for some particular cases of unperturbed flows (e.g. in thermal equilibrium); for unperturbed flows out of equilibrium, with arbitrary density profiles, this interpretation is lost.
As Y 1+Y 2=1, equation (2.21) gives 2.22
When , the pressure gradient and Soret effects vanish from formula (2.14), so that the diffusion velocities are given by Fick’s law (2.15). Because the transport equations for the mass fractions should sum up to the continuity equation (2.7), in this case D12=D21=D. Then, replacing (2.22) in the mass fraction transport equations and using (2.7), it yields that the divergence of velocity is 2.23Because, as , the specific heats lose their temperature dependency, it follows that and . Replacing these relations in the energy equation, dividing the equation by p, and taking the limit , with p/T=constant yields 2.24and the energy equation reduces to the divergence formula (2.23). The , incompressible limit of the NS equations for binary gas mixtures following the ideal-gas equation of state is, thus 2.25 2.26 2.27supplemented by the divergence formula (2.23). Numerical solutions of RTI in the first incompressible case are discussed in §3b(iii).
Equations (2.25)–(2.26) also admit a Boussinesq approximation, at small Atwood numbers. A derivation of the governing equations in this approximation is given in reference . The physical meaning and equations are different than those obtained through the anelastic approximation.
The ‘second’ incompressible limit. In order to provide a basis of comparison for the case when compressibility is studied through changes in the specific heats of the two fluids, an incompressible limit can be defined mathematically as , although, strictly, such limit does not have physical sense. Then, assuming the mixture obeys the ideal-gas equation of state, the limit can be obtained as so that the energy equation reduces to 2.28where τ is the deviatoric part of the viscous stress tensor. The enthalpy diffusion term simplifies to ∇⋅(ρTD∇R); however, the rest of the molecular terms and equations (continuity, momentum and scalar transport as well as the equation of state) keep the same form as in the fully compressible equations. Note that the initial state in the ‘second’ compressible limit can be constructed with any density profile, including one that matches a corresponding compressible case. Several simplified forms of the equations can be constructed by neglecting various molecular effects. In some special cases, some physical meaning can be constructed through an analogy to a politropic transformation, as in reference ; however, these are beyond the scope of this paper. No numerical simulations have been performed for the equations defining the ‘second’ incompressible limit.
3. Numerical solutions to the governing equations
Here, numerical simulations of the miscible RTI at high A, with atomistic methods (MD and DSMC), LBM and with compressible and incompressible NS equations are discussed and results concerning the turbulent mixing problem are highlighted. Even though atomistic simulations are restricted to very small domains and the LBE for two fluid mixing at high A is still under development, the reason to include these types of simulations is twofold. First, these areas are rapidly growing both in simulation sizes and theoretical results. In addition, particle methods are expected to scale better than existent Eulerian approaches for partial differential equations (PDEs), which is especially important to take advantage of the likely future supercomputer architectures. Recent MD simulations remarkably recover, when appropriately averaged over many realizations, the linear immiscible RTI growth calculated based on the NS equations. If this result holds for the miscible case and at long times, in the turbulent stage, and the convergence rate is fast enough, then atomistic simulations may offer an alternative to the compressible NS solvers for investigating the physics of turbulent mixing on the next generation of supercomputers.
Compressible NS equations solvers, on the other hand, pose their own challenges owing to the very large range of dynamically relevant spatio-temporal scales present in the compressible RTI layer. As recent results suggest , besides the intrinsic broad range of turbulent scales, such flow also generates shock waves that further extend the dynamically relevant range of scales. One promising approach for such problems, at least at early to intermediate times, when the turbulence is highly localized, is to use an adaptive mesh refinement technique and such approaches are discussed below. Unfortunately, however, many such techniques pose significant challenges in terms of accuracy, such that the numerical diffusion remains negligible compared with the physical molecular transport terms, have directional bias, and/or do not scale well on massively parallel computers. Even though some adaptive mesh approaches are promising in terms of controlled errors and accuracy, it is not surprising that most results to date concerning the physics of mixing at high A were obtained in the incompressible case, when the acoustic fluctuations are filtered out and the range of dynamically relevant scales is significantly reduced. For the incompressible RTI problem, a number of fully resolved, NS simulations exist, in both the classical and idealized triply periodic configurations, and an overview of these simulations is presented in §3b(iii).
Owing to the significant challenges in performing fully resolved NS simulations for RTI, coarse resolution simulations have also been performed either using numerical diffusion or explicit subgrid models to regularize the equations and model the subgrid behaviour. Such simulations for both the compressible and incompressible RTI are discussed in §§3a(ii) and 3b(iii), respectively.
(a) Atomistic and particle methods
(i) Molecular dynamics and direct simulation Monte Carlo
Recent atomistic simulations have shown that even small-scale systems can exhibit RTI, with many similarities to the continuum limit. References [30,67–69] report both MD and DSMC results, with up to approximately 108 particles for MD, approximately 5.7×108 particles for quasi two-dimensional DSMC and approximately 7.1×109 particles for three-dimensional DSMC. The simulations used a splined version of the Lennard–Jones potential V (r)=ϵ[(r0/r)12−2(r0/r)6] for r less than an inflection point and cubic spline with finite range for larger r. Here, r=∥x∥, where the constants r0 and ϵ were chosen to represent methanol as the light species. Both the miscible and immiscible cases were simulated at Atwood numbers up to 0.867. The gravity was large enough in most simulations to generate non-negligible stratification (see also §3a(i) for a discussion on compressibility effects).
The largest simulations showed good agreement with the linear theory for inviscid compressible fluids and also compared well in terms of several characteristics of the mixing layer (e.g. growth rate, fractal dimension) with incompressible NS simulations and quasi-two-dimensional, controlled initial conditions experiments .
Using the physical properties of the fluids in the simulations, the Kolmogorov microscale needs to be ηK≈7×10−9 m (see (2.18)), in order for the viscous and thermal fluctuations to be commensurate. Because the most unstable modes for the simulations ranged from 40×10−9 m to 240×10−9 m, one should expect an overlap between the viscous scales and the scales where thermal fluctuations are important. Indeed, all individual MD or DSMC results presented showed significant departures from the continuum results at late times, with the edges of the mixing layer exhibiting more roughness than continuum simulations at the same A (comparisons at the same A are important, as the asymmetry of the layer and break-up at the spike side are expected to significantly increase with A, for A>0.5). Significant departures even from the linear theory were observed in smaller immiscible simulations . However, interestingly, when these smaller simulations were averaged over independent realizations, the thermal fluctuations effects vanished and the linear theory growth was recovered. This is remarkable, as it indicates, quantitatively, how the mean of the microscale fluctuations convergences to the macroscale, continuum limit. It would be interesting to see if the late time results also converge and if, at what rate, to the continuum limit. While MD and DSMC descriptions are still very far from being feasible for macroscopic fluid flows, if the convergence rate is fast enough, then these type of simulations may offer an alternative way of investigating the physics of compressible RTI to the usual NS equations, which are themselves a challenge to solve on useful sizes for compressible turbulent flows (see below).
(ii) Lattice Boltzmann method
LBM simulations have become a promising alternative approach for the immiscible RTI problem, because continuum simulations of turbulent interfacial problems are very difficult due to the complex interface topology, interface break-up and reconnection, surface tension effects, etc. Most earlier models were restricted to small A and presented various theoretical challenges ; however, newer models progressively overcome these difficulties . For the miscible case, recent progress has been made for the compressible thermal case, when the background density changes through variations in temperature  and LBM simulations with A up to 0.4 and strong background stratification are reported in reference . These models use higher-order discretizations of the distribution function, but still rely on a BGK-type single relaxation time.
Such models do not allow non-unity Pr or Sc numbers. Luo & Girimaji  examined a two-fluid isothermal model with different relaxation times for momentum and species and also accounting for mutual and self-collisions, leading to Fickian diffusion in the continuum limit. The model was further extended in reference  for species with different molecular weights, in which case the macroscopic diffusion velocities also exhibited the pressure gradient term (2.14). Recently, Shan  developed a BGK-type multiple relaxation mixture model starting from the continuum kinetic theory which also contains the thermal diffusion term in the corresponding diffusion velocity. Even though such models have not been applied to the RTI, their development gives hope to an alternative for the compressible miscible RTI simulations using the NS equations. In particular, further development is required to improve the numerical stability of the LBM models and more theoretical work on modelling species molecular properties.
(b) Continuum methods
Direct numerical simulations (DNS)  are becoming a powerful tool for solving the equations of motion for fluid flows, especially with the current increase in supercomputer power and exascale computing on the horizon. In this technique, the equations are solved using adequate algorithms to obtain a (nearly) grid independent solution with appropriately small numerical errors. In practice, only statistics up to a certain order, depending on the problem of interest, are required to be little affected by numerical errors. DNS are conducted without resort to either turbulence or subgrid modelling or the introduction of ‘artificial’ numerical dissipation or other algorithm stabilizing schemes. This offers a wealth of information pertaining to the physics of turbulent mixing, albeit at relatively low Reynolds numbers and in the presence of at most weak shocks. Note that DNS are not possible for the Euler equations that are known to diverge and require either explicit or implicit stabilization techniques. In addition, low order schemes (e.g. first- or second-order) are not generally used for DNS owing to the excessively large grids required for acceptable accuracy.
In this section, a review of the compressible and incompressible simulations for the miscible RTI is provided. For the purpose of this discussion, fully resolved simulations of the NS equations (DNS), which can provide reliable data for understanding the physics of mixing, are distinguished from Euler equations or under-resolved NS equations solutions. The last two types of numerical simulations can still be useful as tools to probe regimes outside the reach of DNS or to scope out parameters faster than DNS; however, the results need to be appropriately (and critically) considered.
Reviews of compressible and incompressible RTI simulations are also provided in references [66,79]. Here, the focus is on the high Atwood number RTI and mixing problem. Thus, some of the simulations mentioned in these reviews are not discussed here, whereas idealized triply-periodic RTI simulations, which do not appear in previous RTI reviews, are included.
(i) Direct numerical simulations of Rayleigh–Taylor instability with compressible materials
Owing to the range of scales of compressible flows, where the fast acoustic motions need to be solved, as well as increased number of variables compared with the incompressible case, presently, no DNS studies exist for the multi-mode RTI. One solution to reduce the computational effort, especially for the RTI where the flow is highly localized at early and intermediate times, is to use an adaptive mesh refinement techniques. Nevertheless, most such techniques in use today are low order and introduce directional bias, so that special care must be paid in using such techniques for DNS. Two mesh refinement approaches have been used to the study of single mode compressible RTI [79,80] so far at A>0, albeit still at small A values.
Le Creurer & Gauthier  use the AMÈNOPHIS code to study the long time relaxation towards equilibrium of a two-dimensional single-mode RTI in finite box at A=0.2. The code uses a mixed Fourier–Chebyshev polynomial expansion, and the mesh adaptation is based on a rigorous definition of the upper bound of the error. In a separate paper , the same code is used to examine the compressibility effects through changes in the specific heats, γ1 and γ2, as well as initial background stratification. The results are consistent to the linear theory [61,65]. Two-dimensional RTI single-mode simulations are also reported in reference  using the adaptive wavelet collocation method (AWCM) . AWCM uses wavelets for dynamic grid adaptation and has direct error control. The simulations, at A=0.1, investigate the potential flow theory regime. No such simulations exist for the three-dimensional case, at high A values, or multi-mode.
(ii) Large eddy simulations of RTI with compressible materials
The large eddy simulation (LES) technique is based on the paradigm that the large scales, which are influenced by external forcing, boundaries, etc., are simulated, whereas the small scales, presumably more universal, can be modelled. The subgrid modelling can be performed through explicit physical models , or handled, implicitly, by the features of the numerical algorithm , or a combination of both. Implicit LES (ILES) is attractive owing to its simplicity; however, the flows solved have nominally the same type of subgrid transfer for mass, momentum and energy, so that the Sc and Pr numbers are unity (different than unity values require different numerical schemes in the momentum, energy and species equations). In addition, it is very difficult to assess the nature of the transport processes. Nevertheless, such simulations are increasingly more common, including for RTI. Most of the LES of RTI at A>0 have actually been performed with compressible codes at low values of Mach number/stratification parameters. A complete review of such codes can certainly exceed the size limit of this article. Two points are necessary. First, all LES of compressible RTI, even at very low compressibility levels, are considered as belonging to this section if they were performed with compressible codes. Second, only simulations which discuss the physics of the mixing, beyond the value of the growth parameter, α, have been included here. For a critical assessment of α values in numerical simulations, the reader is directed to reference . Additional references concerning compressibility effects can be found in reference .
In reference , a dynamic mixed model is used to perform simulations at A=0.5 with the MIRANDA code (see below) and study the compressibility effects as a function of the turbulent Mach number. For the duration of the simulation, the turbulent Mach number remained below 0.6, so they concluded that compressibility effects, through turbulent fluctuations, are small. However, Olson & Cook  performed longer simulations, also with the MIRANDA code, but with a different subgrid model  and found an unexpected result: the pressure waves generated by the piston-like motion of the bubbles can coalesce into strong shock waves. Aside from underlying the important role compressibility can play on the RTI development, the results may also have interesting implications on the detonation initiation in type Ia supernovae or X-ray bursts, especially if a favourable temperature gradient is present. The MIRANDA code uses a hybrid Fourier-10th order compact finite differencing and also has a ‘DNS’ mode. Such simulations are discussed below. George & Glimm  report multi-mode simulations, at A=0.5, with the finite volume front tracking code Frontier, but in the untracked mode. They observe that the effective Atwood number changes owing to the mixing and defined an effective, time-varying Atwood number that collapses previous numerical and experimental data on the mixing layer growth rate. In reference , both single- and multi-mode RTI 1203 simulations are reported with the MAH-3 code. The development of the kinetic energy spectrum is discussed in connection to the onset of self-similarity in the layer growth. Nevertheless, the simulations are too coarse to draw definite conclusions. Thus, for the compressible case, results on the structure of turbulence, spectral behaviour or mixing characteristics remain scarce. The physics of turbulent mixing still represents, largely, an open problem in the compressible case.
(iii) Direct numerical simulations and large eddy simulation of Rayleigh–Taylor instability with incompressible materials
Owing to the reduced number of variables compared with the compressible case and time steps not constrained by acoustic phenomena, very high resolution simulations have been performed for the incompressible RTI (all of them using the ‘first’ incompressible limit defined above), in both the classical RTI configuration and an idealized triply periodic configuration, named homogeneous Rayleigh–Taylor (HRT). HRT allows the study of turbulence and mixing without the complications owing to inhomogeneities or mixing layer edges, but retaining the basic nonlinearities and energy production mechanism as the classical RTI case.
DNS results for the incompressible RTI problem are reported in references [15,25,87] with the MIRANDA code on up to 30723 meshes and references [20,88,89] with the CFDNS code on up to 40962×4032 meshes. All MIRANDA runs have been performed at A=0.5, whereas CFDNS has been used for A=0.04–0.9. The CFDNS code also uses a mixed spectral-compact finite differences scheme. DNS of HRT are reported in references [1,90–93] (A up to 0.5) for the unsteady case and reference  for the stationary case (A up to 0.9). LES of incompressible RTI using explicit subgrid modelling were performed with the MIRANDA code in reference  at A=0.5 and of stationary HRT using the stretched vortex model of Pullin  in reference  at A up to 0.9.
The above simulations examine in great detail the turbulence and mixing characteristics in the RT mixing layer. Below, a brief summary of some of the findings is provided.
One of the more remarkable issues explored in incompressible RTI, first reported in reference , is the marked difference in the mixing between different density fluids as opposed to mixing that occurs between fluids of commensurate densities, corresponding to the Boussinesq approximation. Thus, at large density ratios, the mixing becomes asymmetric, with the pure heavy fluid mixing more slowly than the pure light fluid. The existence of this asymmetry is significant in practical applications. For example, it is important to know both how long it takes for a pollutant to mix with the surrounding fluid and also how much is likely left at a certain time. Predicting this based on Boussinesq analogies, as it is usually done, can severely misrepresent the mixing of the pollutant. One consequence of the mixing asymmetry for the RT layer is that the penetration distance of the pure heavy fluid is larger than that of the pure light fluid . Figure 1 shows the PDF of the density at several locations across the mixing layer at A=0.75 . The density PDF varies considerably across the RT mixing layer. At the top of the RT layer, the PDF is spiked at the heavy fluid end and includes some mixed fluid. At the bottom of the layer, the PDF is spiked at the light fluid end. The PDF is clearly asymmetrical with respect to the centreline and the amount of pure heavy fluid reaching the centreline is larger than that of pure light fluid. At the edges of the layer, the consequence of the mixing asymmetry is even more obvious: the two edges develop into spikes on the light fluid side and bubbles on the heavy fluid side. This can be clearly seen in figure 2, which shows the density field from a 40962×4032 RT simulation at A=0.75 , the largest fully resolved instability simulation to date.
In general, all mix metrics in use today for RT turbulence are constructed from lower-order moments of the density PDF. A comprehensive discussion of these metrics, including rigorous bounds, can be found in reference . Thus, the usual mix measure θ=〈f1f2〉/F1F2 depends on the mean and variance of the density. Here, f1=(ρ−ρ1)/(ρ2−ρ1), f2=(ρ2−ρ)/(ρ2−ρ1) and Fl=〈fl〉, . While the density variance does appear in the dynamical equations in the Boussinesq limit , the moments equations at higher density ratios do not contain any term depending on 〈ρ2〉. From the point of view of the dynamical equations, a more appropriate mix metric at high A would be constructed from the density-specific volume covariance, θρv=1−b/bnm, where b=〈ρ′v′〉 and the no-mix value of b is . Here, overbars represent statistical means and primes fluctuations with respect to these means. Figure 3 shows that the two metrics are close at A=0.04, as expected, but they are different at high A and differences increase with A. Thus, there is a qualitatively different behaviour across the RT layer, as θ and θρv predict more mixing at the opposite sides of the layer. Nevertheless, both metrics have relatively large values across the whole layer, which misrepresents the density PDF. Nor they can capture any asymmetry in the underlying PDF. Other mix measures, for example, using a fast reaction analogy, 〈XP〉 , or the time-dependent Atwood number used in reference , whereas useful in certain instances, are still only low order representations of the underlying density PDF.
The spectral behaviour of the kinetic energy, density and mass flux has also been examined in the DNS references given above. In general, it seems that the energy develops an inertial range with a scaling close to −5/3, whereas the mass flux follows a −7/3 scaling . Note that even at the high resolutions used, the exact spectral scaling is difficult to infer. An interesting result, first discussed in reference  for HRT and  for RTI and further confirmed in reference [16,20] is that there is an anomalous anisotropy at small scales (even though the flow develops an isotropic inertial range). This anisotropy may have implications for the subgrid modelling of buoyancy driven flows.
Another important aspect explored in reference  is the existence of a mixing transition in incompressible RTI. Such a transition could have very significant applications in modelling RTI-driven flows, as it would establish a minimum Reynolds number which needs to be attained in the simulations such that no qualitative changes in the physics of mixing occur by further increasing the Reynolds number . Whether such a transition is unique over many classes of flows and whether higher-order mixing statistics have separate transitions are important open questions.
4. Summary and concluding remarks
Numerical simulations of two fluid gaseous mixtures, at high density ratios, as relevant to the RTI have been reviewed. Because such simulations cover a large range of scales, from small atomistic systems to astrophysical situations, a section of the paper has been devoted to presenting the various frameworks of approximation for describing material mixing in fluid flows, as they have been used in the study of the RTI. Each of these frameworks fits into a hierarchy of theoretical models describing the material behaviour at a certain level. Often, in passing to the next (coarser) level, the description becomes unclosed and closure relations are required. Thus, for each theoretical model presented, its assumptions and level of description have been stated. It is hoped that such presentation can introduce to various communities the equations solved at different levels of material behaviour description, to better appreciate the overall assumptions and limitations.
The theoretical models considered start from the time-dependent Schrödinger equation, continue with atomistic descriptions, Boltzmann equation, lattice Boltzmann equation and finally the NS equations. The molecular transport terms in these equations require knowledge of the collision term or potential in kinetic theory descriptions and are often found after complicated derivations. In many situations, experimental calibrations are preferred and these are also stated. Some of the effects present in these terms, not examined in previous studies, may nevertheless be important under certain conditions.
The compressible to incompressible limit of the NS equations is not unique. However, for the RTI, it is especially important to formulate the appropriate basis of comparison when compressibility effects are discussed, because there is more than one compressibility parameter and the initial configuration (e.g. if stratification is present) may obscure or exacerbate these effects. Thus, two incompressible limits, encompassing a large range of realizable initial configurations, are derived. In particular, it is shown that the ‘first’ incompressible limit corresponds to the equations used in some recent large incompressible RTI simulations and that the molecular transport terms in those equations do not require additional assumptions compared with the compressible NS equations (aside from the incompressibility). The ‘second’ incompressible limit allows more general initial configurations, but also retains almost the same formulation of the transport terms as the compressible equations. This was derived here for the first time; no simulations of these equations have been performed.
The second part of the review discusses numerical solutions to the governing equations for the RTI: MD, DSMC for the Boltzmann equation, the lattice Boltzmann method (LBM) for the lattice Boltzmann equation, and finally methods for the NS equations. Here, a distinction is made between accurate, fully resolved simulations of the NS equations, usually called DNS and coarse mesh and/or solutions where the physical transport terms are dominated by subgrid models or numerical errors. The latter are usually called LES when the subgrid model is explicit and the numerical errors negligible or ILES when the subgrid model is mostly provided implicitly, by the numerical algorithm. Fully resolved simulations of the compressible RTI are very challenging, due the broad range of dynamically relevant scales present, especially considering the unexpected recent finding that even pure RTI starting from rest may generate shock waves. Therefore, presently, such simulations do not exist for the three-dimensional multi-mode case. However, two promising mesh refinement techniques are highlighted. Most previous RTI simulations, even at very low compressibility levels, were actually performed with codes solving the compressible equations. These simulations were included in the compressible NS simulations section. Nevertheless, LES of the RTI have been discussed only when addressing issues related to turbulent mixing, aside from the mixing growth rate, which has been extensively presented in previous reviews.
Owing to significant reduction in computational requirements, most in-depth computational results for material mixing in the context of the RTI have been obtained in the incompressible case. The last part of the paper surveys these results, including some unexpected findings related to the asymmetry of the mixing and small scale anomalous anisotropy, as well as the possibility of a mixing transition in RTI turbulence. These results may provide a basis for future studies, especially in regimes outside the continuum, incompressible limit.
Finally, it is reiterated the hope that putting together results and approaches from different levels of description of RTI-driven material mixing may introduce to different communities ideas and provide increased motivation for furthering the study of such flows.
This publication was made possible in part by funding from the LDRD programme at Los Alamos National Laboratory through project number 20090058DR.
Los Alamos National Laboratory is operated by the Los Alamos National Security, LLC for the US Department of Energy NNSA under contract no. DE-AC52-06NA25396.
One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales II’.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.