## Abstract

Simultaneous advances in numerical methods and computing, theoretical techniques, and experimental diagnostics have all led independently to better understanding of Rayleigh–Taylor (RT) instability, turbulence, and mixing. In particular, experiments have provided significant motivation for many simulation and modeling studies, as well as validation data. Numerical simulations have also provided data that is not currently measurable or very difficult to measure accurately in RT unstable flows. Thus, simulations have also motivated new measurements in this class of buoyancy-driven flows. This overview discusses simulation and modeling studies synergistic with experiments and examples of how experiments have motivated simulations and models of RT instability, flow, and mixing. First, a brief summary of measured experimental and calculated simulation quantities, of experimental approaches, and of issues and challenges in the simulation and modeling of RT experiments is presented. Implicit large-eddy, direct numerical, and large-eddy simulations validated using RT experimental data are then discussed. This is followed by a discussion of modeling using analytical, modal, buoyancy–drag, and turbulent transport models of RT mixing experiments. The discussion will focus on three-dimensional RT mixing arising from multimode perturbations. Finally, this focused review concludes with a perspective on future simulation, modeling, and experimental directions for further research. Research in simulation and modeling of RT unstable flows, coupled with experiments, has made significant progress over the past several decades. This overview serves as an opportunity to both discuss progress and to stimulate future research on simulation and modeling of this unique class of hydrodynamically unstable turbulent flows.

## 1 Introduction

Rayleigh–Taylor (RT) instability [1–3] and mixing arise from the growth of typically multimode perturbations with some range of amplitudes and wavelengths present at an interface initially separating liquids, gases, or solids of different densities in an unstable configuration in which the lighter, lower density fluid is accelerated toward the heavier, higher density fluid [4]. The unstable configuration corresponds to misaligned density and pressure gradients that locally satisfy $\u2207\rho \xb7\u2207p<0$. The initial perturbations can arise from very small-scale random motions in the velocity or density fields present naturally or imposed through an external mechanism (e.g., mechanical vibration). The power spectrum of the perturbations can be primarily narrowband or broadband.

If surface tension is present, the fluids are immiscible; otherwise, a diffusion layer is either initially present or forms following the removal of any barrier and the fluids are miscible. In the immiscible and miscible case, the mixing process is an interpenetration and molecular, respectively. In the linear growth phase, the modes grow independently until their amplitudes become large enough that the modes begin to couple in the weakly nonlinear regime. Provided that the Reynolds number becomes large, the flow eventually becomes turbulent, with the lighter bubble-like structures interpenetrating into the heavier spike-like structures. Shear between these structures breaks down the structures into progressively smaller scales. In miscible mixing, the interface area increases dramatically with time. The acceleration (e.g., gravity) provides a sustained driving force, with the gravitational potential energy converted into kinetic energy, expanding the extent of the mixing layer in the direction parallel to the acceleration.

*t*

_{0}and a virtual origin $h(0)=\alpha \u2009At\u2009g\u2009t02$]

(see Ref. [9] for additional discussion). RT flows with larger density contrasts (larger At) lead to a wider mixing layer. Much of the experimental, computational, and theoretical work on RT flow has focused on the determination of $\alpha b,s$ and *α*, which exhibit a wide range of experimental and simulation values (see Ref. [10] for a summary).

Over the last several decades, experimental measurements and numerical simulations with various initial conditions, fluid combinations, and accelerations indicate that the values of *α _{b}* and

*α*(and therefore, of

_{s}*α*) are not universal, but depend on the spectral characteristics of the initial conditions, Atwood number (particularly for the spikes), and on other flow parameters. The instability growth rate depends on surface tension (for immiscible fluids), viscosity, mass diffusivity, Schmidt number, compressibility, yield strength (in elastic–plastic flows), magnetic fields, rotation, convergent geometry, reactions etc. (see Refs. [10] and [11] for a survey of these topics). As the Atwood number increases, the bubble and spike fronts become increasingly asymmetric; this asymmetry is also manifested in the spatial profiles of fields and terms in transport equations across the layer.

The majority of the results discussed here pertain to incompressible mixing. However, even in compressible RT unstable flows, the turbulent Mach number, $Mat=2K/cs$ (*K* is the turbulent kinetic energy and *c _{s}* is the mixture sound speed) will be small. Therefore, incompressible experimental results are still relevant to compressible applications. It is also important to note the following paradox with respect to fluid miscibility in experiments and simulations. While the principal objective of RT simulation and modeling for scientific and engineering applications is to describe miscible mixing, many experiments which have been the subject of such simulation and modeling have been immiscible in order to minimize early time fluid mixing or to utilize fluid combinations with large Atwood numbers. Consequently, measurements in such immiscible experiments are likely influenced by the stabilizing effect of surface tension. The simulations and modeling (except for simulations using front-tracking) performed thus far have not included surface tension effects.

Rayleigh–Taylor unstable flows exist in nature and in applications. RT mixing is central to the deflagration-to-detonation transition of thermonuclear flames in type Ia supernovae [12–14]. RT instability and mixing in inertial confinement fusion [15,16] capsules decreases the efficiency of thermonuclear burning [17–19] and is one of the mechanisms preventing gain and ignition of thermonuclear fuel in the laboratory [20]. In geophysics, RT instability exists in the lithosphere [21–23]. The accurate and predictive modeling of the effects of this instability on the flows in these examples is critical, and numerical simulation and experimental data have had a significant role in improving understanding of these effects and modeling them using a hierarchy of approaches.

In the examples briefly enumerated above, RT flow and mixing are only one of many complex hydrodynamic and physical effects that must be accurately described. Despite ever-increasing computational resources, it remains impractical to construct and use highly complex simulation models for applications to problems in nature and in scientific and engineering applications. As a result, a major effort consists of constructing, validating, implementing, and using reduced-order models that can describe the most important effects of complex phenomena such as buoyancy-driven RT instability and turbulent mixing. For this purpose, simplified numerical simulation models are constructed, run, and analyzed to provide detailed insight into the physical mechanisms governing RT flows. The veracity of such simulation and reduced-order models must ultimately be established by comparison of their predictions to physical experimental data. This review summarizes such efforts and comparisons to date, and identifies gaps and areas of further research. It should be noted that simulations, models, and experiments engender many inherent idealizations, simplifications, and limitations: further progress requires recognizing and addressing these as advances in theory, numerical methods, computation, and experimental diagnostics are made.

Significant progress on simulation and modeling of Rayleigh–Taylor instability-induced turbulent mixing has been made over the last several decades. This article will briefly discuss advances in simulation and modeling synergistic with experiments that have led to an enhanced understanding of the physical mechanisms in Boussinesq and variable-density RT turbulent mixing. Although considerable simulation and modeling research has been performed to investigate single-mode, two-dimensional, and immiscible RT instability, these will not be considered here (see Refs. [10] and [11] for a review of these topics). A parallel review of advances in diagnostic techniques and experiments on RT instability and mixing complements this review [9].

This paper is organized as follows. A brief overview of quantities measured in experiments and calculated in simulations is given in Sec. 2. Factors that complicate simulation and modeling of RT flow experiments are discussed in Sec. 3; experimental approaches are briefly discussed in Sec. 3.1 and simulation and modeling issues are discussed in Sec. 3.2. The numerical simulation of RT flow and mixing experiments is discussed in Sec. 4. First, the governing equations and the numerical approaches used are briefly described in Sec. 4.1. A review of progress through applications of these numerical approaches to three-dimensional RT flow experiments is then given in Secs. 4.2–4.4. The modeling of Rayleigh–Taylor flow and mixing experiments is discussed in Sec. 5. First, modeling approaches that have been used are briefly described in Sec. 5.1. A review of progress through applications of these modeling approaches to three-dimensional RT flow experiments is then given in Secs. 5.2–5.4. Finally, a synthesis and future directions for simulation, modeling, and experiments are discussed in Sec. 6.

## 2 Measured Experimental and Calculated Simulation Quantities

Let the coordinate direction aligned with the acceleration be *z* and let the vertical domain be $z\u2208[0,Lz]$. Also, let the subscript *r* label the fluid (or species). Numerical simulations of RT flows have investigated (not an exhaustive list):

mixing layer widths ($hb,s$ or

*h*) and growth parameters ($\alpha b,s$ or*α*) based on cutoffs in some mean field on the bubble and spike front sides;- integral mixing layer widths (the overbars indicate averages over the directions perpendicular to the acceleration, and
*n*is a numerical factor) [24]_{W}for scalars $\varphi r={mr,fr,Xr,cr}$ (mass fraction, volume fraction, mole fraction, or concentration, which are dimensionless and satisfy $\u2211r\varphi r=1$ and $\u2211r\varphi r\u2032=0$);$W(t)=nW\u222b0Lz\varphi \xaf1\u2009\varphi \xaf2\u2009dz$(2) mean mass fraction, volume fraction, and concentration profiles $m\xafr(z,t),\u2009f\xafr(z,t)$, and $c\xaf(z,t)$;

- local and global molecular mixing profiles and parameters ($\varphi 1\u20322\xaf$ is the scalar variance in fluid 1) [25,26]$\Theta (z,t)=\varphi 1\u2009\varphi 2\xaf\varphi \xaf1\u2009\varphi \xaf2=1\u2212\varphi 1\u20322\xaf\varphi \xaf1\u2009\varphi \xaf2\u2009$(3)respectively, where $\Theta \u0302$ and $1\u2212\Theta \u0302$ are also referred to as “atomic” and “chunk” mixing fractions, respectively;$\Theta \u0302(t)=\u222b0Lz\varphi 1\u2009\varphi 2\xaf\u2009dz\u222b0Lz\varphi \xaf1\u2009\varphi \xaf2\u2009dz=1\u2212\u222b0Lz\varphi 1\u20322\xaf\u2009dz\u222b0Lz\varphi \xaf1\u2009\varphi \xaf2\u2009dz\u2009$(4)
- ratios of kinetic energy produced to released potential energy and of dissipation to released potential energy $K/\Delta \Phi $ and $D/\Delta \Phi $, respectively [25,27], where the released potential energy is the difference between the final and initial potential energies (assuming a step function for the initial density at the interface location $Lz/2$),with(5)$\Delta \Phi =\Phi initial-\Phi final=gLz[(\u222b0Lz/2\rho 2\u2009z\u2009dz+\u222bLz/2Lz\rho 1\u2009z\u2009dz)-\u222b0Lz\rho \xaf\u2009z\u2009dz]$
*g*>and the dissipation is the difference between the released potential energy and kinetic energy produced, $D=\Delta \Phi \u2212K$;$K=1Lz\u222b0Lz\rho \u2009v2\xaf2\u2009dz$(6) - directional components of the kinetic energy spectrum and the scalar variance spectrum, and their spatial integrals [$vi(k,z,t)\u2032$ and $\varphi (k,z,t)\u2032$ are Fourier-transformed velocity and scalar field fluctuations, and $k=kx2+ky2$ is the two-dimensional wavenumber](7)$Ei(k,t)=1h(t)\u222b-hs(t)hb(t)|vi(k,z,t)\u2032|22\u2009dz$(8)$E\varphi (k,t)=1h(t)\u222b-hs(t)hb(t)|\varphi (k,z,t)\u2032|22\u2009dz$$Ei(t)=\u222b0\u221eEi(k,t)\u2009dk,\u2009\u2009E\varphi (t)=\u222b0\u221eE\varphi (k,t)\u2009dk$(9)

Early simulation studies focused on quantities that could be measured experimentally at the time. Early experiments mainly used visualization techniques to infer the mixing layer width. Later, improved diagnostics coupled with methods to track each fluid enabled measurements of the degree of molecular mixing [9]: where $\Theta (z,t)=0$ (corresponding to $\varphi 1\varphi 2\xaf=0$), the two fluids are completely segregated, and where $\Theta (z,t)=1$ (corresponding to $\varphi 1\varphi 2\xaf=\varphi \xaf1\varphi \xaf2$), they are completely molecularly mixed. With the advent of direct numerical simulation (DNS) (see Sec. 4.1.2.3), quantities studied in other areas of turbulence, such as spectra and statistics, were considered. Many of the quantities enumerated above are sensitive to the grid resolution, and have been reported as a function of resolution in many simulation studies.

## 3 Experimental Approaches and Issues for Simulation and Modeling of Rayleigh–Taylor Flow Experiments

As discussed in the review of experiments [9], there are several classes of designs that can be generically categorized as tank and channel experiments. As experimental designs have evolved, additional flow diagnostics were implemented and calibrated in order to measure additional statistical quantities. Following an overview of experimental approaches, factors that make simulation and modeling of RT instability and turbulence experiments challenging are reviewed.

### 3.1 Overview of Experimental Approaches.

The number of multimode RT instability and mixing experiments that have been simulated or modeled is quite small. Researchers in experimental groups, as well as their collaborators and others, have performed simulations and modeling of various aspects of these experiments both to validate and improve computational and theoretical models of RT flow, turbulence, and mixing. See Ref. [9] for additional background on experimental diagnostics and their limitations.

#### 3.1.1 Accelerated Tank Experiments.

In accelerated tank experiments such as performed on the rocket-rig [28–30] (also see Ref. [31]) or on the linear electric motor (LEM) [32,33] (which can achieve accelerations much larger than Earth gravity), the acceleration must be measured sufficiently accurately, and the mechanism for initially separating the fluids becomes important for the early time and nonlinear growth of the instability. Surfactants were almost always used and the majority of experiments were immiscible. If a diffusion layer is present initially, the instability growth will be modified by the binary molecular diffusion of one fluid into the other [34]. These complications, coupled with the limited optical diagnostics that were available at the time, made it difficult to extract more than the mixing layer width using backlit photography. Thus, simulations and modeling of such experiments have focused on reproducing the measured values of *α _{b}* and

*α*. Another challenging aspect of the rocket-rig experiments is that the introduction of asymmetry in the initial condition by tilting the rig led to the development and evolution of combined Rayleigh–Taylor and Kelvin–Helmholtz instability [30], which minimally requires two-dimensional simulation or modeling, as discussed in Secs. 4.2 and 5.2, respectively.

_{s}In experiments with variable (or complex) accelerations, simulations using parameterizations of the acceleration have been performed in order to reproduce the observed dynamics. These experiments have also inspired many numerical simulation and modeling studies exploring simplified piecewise constant unstable and stable acceleration cases, which are not discussed here. The initial perturbations in the accelerated tank experiments considered here were random, with no precise quantitative characterization of their modes due to limitations in the diagnostics available at the time (although as noted in Secs. 4.2.1 and 5.2.2, initial perturbations were estimated from photographs taken at early times). The fluid combinations that could be used were limited by the need for refractive index matching (as also for stationary tank and channel experiments): a limitation of the diagnostics. Since these early experiments were performed, advances have been made in the generation and characterization of the initial perturbations present in tank experiments (e.g., see Refs. [35] and [36]).

#### 3.1.2 Stationary Tank Experiments.

In stationary tank experiments, gravity is that of the Earth, but again the method of segregating the fluids initially introduces various nonideal aspects, e.g., the removal of a barrier, which produces a viscous boundary layer on the upper and lower surfaces, generates strong vortices propagating vertically away from the barrier, and leaves a wake. In early experiments using a tank that (with a fluid configuration that was initially stable) was rapidly overturned to initiate RT instability without or with a tilted interface [24,37], the overturning motion introduced shear and hysteresis effects that were challenging to estimate and model accurately at the time. In the case of experiments with barrier removal [38,39], the effect of the wake (which imposes a large initial two-dimensional perturbation) was modeled using potential flow-based expressions.

Experiments performed in stationary tanks are typically miscible, and measurements of the molecular mixing could be performed using various techniques. In addition to investigating the large-scale dynamics (such as the mixing layer width) simulations of such experiments computed local and global molecular mixing parameters in order to compare with measured values. These experiments have very short data collection times and require a large ensemble of runs, significantly limiting the acquisition of statistical data. When the mixing layer fills the vertical extent of the tank, an overturning of the two fluids occurs, establishing a stable stratification with velocity and density fluctuations decaying to zero. Simulations and modeling of such experiments are discussed in Secs. 4.3 and 5.3, respectively.

#### 3.1.3 Channel Experiments.

In experiments using a water channel [40–46] or a gas channel [47–51], a thin splitter plate separates coflowing liquids or gases and an inflow-type initial condition is present. In many early simulations of such experiments, it was assumed that at sufficiently late times, the flow becomes self-similar and independent of the details of the initial conditions. This assumption effectively circumvented the need to specify detailed and precise initial conditions in the simulations. However, as simulations could be performed at higher spatial resolutions and for longer times, it was found that the flow structures and statistics were indeed dependent on the initial conditions. Also, the flow over the splitter plate induces boundary layer growth near the plate, with the developing interface beyond the plate perturbed such that growing bubble and spike structures are distorted until phase randomization and momentum diffusion downstream ensue: modeling this boundary layer growth remains a topic of research.

Water channel experiments (which are transitional or perhaps weakly turbulent) used a small temperature differential between the upper and lower streams, which was equivalent to a small density differential that initiated small Atwood number RT mixing downstream (see [52]). Such experiments also used fresh water and salt water streams to create a small density differential [46,52]. These experiments can be described by the Boussinesq approximation. Gas channel experiments using different gas combinations later overcame the limitation to small density ratios, allowing for measurements at larger Atwood numbers to explore non-Boussinesq effects, as well as for larger Reynolds numbers. Molecular mixing parameters and other statistics were also measured in both types of experiments.

Channel experiments rely on the Taylor hypothesis [53] to convert downstream distance to time through $x=vc\u2009t$ (where *v _{c}* is the constant convective velocity). The same correspondence is used to relate time in a simulation or model and downstream distance in an experiment. When the velocities of the streams are matched, the flow is pure RT, while a combined RT/shear flow can be created using unequal velocities (varying the Richardson number). Particle-image velocimetry (PIV) and planar laser-induced fluorescence (PLIF) allow the simultaneous measurement of velocity and density fluctuations. As compared to stationary tank experiments, this class of experiments permits much longer data acquisition times (> 10 min), produces statistically convergent measurements, and does not have long-time transients. An important distinction between stationary tank experiments and channel experiments is that the initial fluids are at rest in the former, while the fluids are intially coflowing in the latter. The simulation and modeling of such experiments are discussed in Secs. 4.4 and 5.4, respectively.

### 3.2 Issues in Simulations and Modeling.

There are a number of issues that arise in performing simulation and modeling of RT flow experiments, many of which complicate the comparison to experimental results, and are briefly outlined here. Simulations and modeling of RT experiments have considered all three types of experimental configurations. Most of the modeling of RT experiments thus far has focused on predicting the variation of *α _{b}* and

*α*with Atwood number.

_{s}#### 3.2.1 Nonuniqueness.

Many quantities arising in simulation and modeling of RT flows are not uniquely defined. For example, the mixing layer widths $hb,s(t)$–which are large-scale flow quantities determining the average boundaries of the layer–are not uniquely defined. Many possible definitions exist, and have different values (e.g., based on mass or volume fraction cutoffs, integrated products of mass or volume fractions, etc.). For example, many experiments use a 5–95% or 10–90% cutoff to define the mixing layer while many simulations use a 1–99% cutoff, and proper comparisons between simulation and experimental data require consistent definitions and conventions. Integral widths, *W*, are often used as they are smoother than widths computed using cutoffs.

As the mixing layer growth parameters $\alpha b,s(t)$ are defined through $hb,s(t)$, the nonuniqueness of the definition of the mixing widths leads to nonuniqueness in the definitions of $\alpha b,s(t)$ [10,54] that can give considerably different estimates of these key growth parameters (see Ref. [55] for a discussion based on a simulation study). Expressions for the growth parameters that have been used in simulations and modeling are:

- the Youngs [8,25,27,31,60,61] expression(where the scaled acceleration distance is $X=At[\u222b0tg(t\u2032)\u2009dt\u2032]2=At\u2009g\u2009t2$ for a constant acceleration), and is an estimate from the late-time slope of $hb,s$ plotted versus(11)$\alpha b,sY(t)=dhb,s/dt2\u2009At\u2009g\u2009t=dhb,sdX\u2009$
*X*, which allows for an offset; - the Ristorcelli–Clark [62] expression(12)$\alpha b,sRC(t)=(dhb,s/dt)24\u2009At\u2009g\u2009hb,s(t)\u2009$
and;

- the Livescu et al. [63] expressionwhich avoids computing time derivatives and allows for an offset in time $t0$.(13)$\alpha b,sL(t)=[hb,s(t)\u2212hb,s(t0)g\u2009At(t\u2212t0)]2$

In a pure self-similar state, these expressions are equivalent. However, each approaches that state at a different rate and with different magnitudes (and no experiments or simulations have attained a pure self-similar state). Other expressions based on centerline velocity fluctuations have been used in experimental investigations (see Ref. [9]). The growth parameters can also be computed using expressions for *W* instead of *h*. Subtleties associated with the offset, particularly with respect to estimates of $\alpha s/\alpha b$ and $hs/hb$, are discussed in detail elsewhere [55].

#### 3.2.2 Range of Scales.

Rayleigh–Taylor instability is initiated at perturbed interfaces separating different-density fluids from perturbations that can have small or large amplitudes and a wide range of wavelengths. Resolving the initial perturbations is essential for accurately describing linear instability growth and subsequent nonlinear mode coupling. RT flows develop progressively larger-scale structures through bubble competition, corresponding to an inverse energy cascade, i.e., transfer of kinetic energy from small to larger scales (see Ref. [10]). As a result, a fully developed RT flow contains a very broad range of interacting scales from large-scale bubble and spike structures to small-scale eddies generated by shear.

*ν*is the mixture kinematic viscosity and neglecting

*t*

_{0})

*α*, At,

*g*, and

*ν*provides the values of Re

*attainable for various values of these parameters at a given time*

_{h}*t*. With the mixture Schmidt number $Sc=\nu /D$ (

*D*is the mixture mass diffusivity) the scalar analog of Eq. (14) is the mixing layer Péclet number $Peh(t)=Sc\u2009Reh(t)=h(t)vh(t)/D$.

Combining the estimate of the ratio of the integral scale to the Kolmogorov viscous dissipation scale from homogeneous incompressible turbulence, $\u2113(t)/\u2113d(t)\u223ch(t)/\u2113d(t)\u223cReh(t)3/4$, with Eq. (14) shows that the Kolmogorov scale $\u2113d(t)\u2248[\nu 3/(8\u2009t)]1/4/\alpha \u2009At\u2009g$ decreases very slowly with time [62,64,65] (also see Ref. [54]). The analogous diffusive lengthscale is the Corrsin–Obukhov scale $\u2113d\varphi (t)=\u2113d(t)/Sc3/4$ [66], and larger Sc favors smaller Corrsin–Obuhkov scales. The grid spacing in a DNS (see Sec. 4.1.2.3) of a turbulent flow with viscous dissipation and mass diffusion is determined by $\u2113d$ and $\u2113d\varphi $. In particular, the requirement to resolve the Corrsin–Obukhov scale sets a stringent resolution requirement for a DNS [67] of a water channel experiment using hot/cold water with Sc* *=* *7 [45], and a DNS of a Sc* *=* *620 water channel experiment using fresh/salt water [46] is currently not feasible. However, gas channel experiments have $Sc\u223c0.7$ [48], making them more attractive for DNS.

Assuming equal grid spacings in each spatial direction $\Delta x=\Delta y=\Delta z$, this implies that the grid spacing should minimally satisfy

#### 3.2.3 Atwood Number Dependence.

*α*is a weak function of At,

_{b}*α*strongly depends on At [33,58,61,68,69] (also see Ref. [9]). Immiscible mixing experiments [33] indicated that

_{s}with $c2\u22480.33\xb10.05$ and $c1\u22481$; similar estimates have more recently been obtained for miscible mixing in gas channel experiments (see Ref. [9]). However, miscible large-eddy simulations (LES) of large Atwood number RT mixing indicate a smaller value $c2\u22480.185$ and $c1\u22480.958$ [58] and $c2\u22480.18$ and $c1\u22481.02$ [69], which implies less asymmetry between the spike and bubble fronts than indicated by experiments. Very few numerical simulations have been performed over a wide range of At that achieve a state where *α _{b}* and

*α*both asymptote to provide a reliable estimate of the ratio (16) or of $hs/hb$. At small At, the bubble and spike mixing layer fronts are nearly the same. At large At, the spike front is larger than the bubble front, further exacerbating the vertical domain length and zoning needed on the spike side (and in the free-fall limit $At\u22721$, the spike structures become narrower with tips having infinite curvature, as shown by potential flow theory [70]). This is one of the principal reasons that most simulations have $At\u22640.5$, although simulations at larger values have been increasingly performed in the last decade.

_{s}#### 3.2.4 Initial Perturbations.

Experimental and numerical evidence for RT flows suggests that the details of the initial perturbations are not “forgotten” at late times, implying that the statistical properties of a given flow depend in some way on those details [71]. Most studies have used an analytically specified initial perturbation spectrum with some distribution of modes having random amplitudes and a specified standard deviation. The dynamics of RT instability growth are different between initially sharp and diffuse interfaces [34]. As the initial conditions in RT experiments are very difficult to impose, measure, and control [9], it has been historically challenging to initialize numerical simulations with perturbations and other features (e.g., boundary layers and wakes) closely matching those in physical experiments.

Realistic initial perturbations that model those present in physical experiments are quite different from the simpler initial conditions commonly used in most RT simulations. In general, experiments have both initial density and velocity perturbations: simulations typically use only one type of perturbation, although given one type of perturbation, linear theory can be used to deduce the other. Experiments also have asymmetry in the initial conditions manifested through different perturbation spectra in different directions: this means that the perturbations cannot be completely characterized using an isotropic wavenumber $k=kx2+ky2$, as used in most simulations. The initial perturbation spectra must be constructed using density and velocity fluctuation data measured using techniques such as PLIF and PIV, respectively [44,45].

*z*-direction is parallel to the acceleration, and periodic boundary conditions are taken in the

*x*- and

*y*-directions) [26,72,73]:

*m*and

*n*are integers) with wavelengths ${\lambda x,\lambda y}={Lx/m,Ly/n}$, the amplitude coefficients ${amn,bmn,cmn,dmn}$ are sampled from a random distribution, a given power spectrum

*P*(

*k*) is assumed, and a root-mean-square amplitude $zrms=\u222bkminkmaxP(k)dk$ is specified [61]. The density perturbations can also be converted to heavy mass fraction or volume fraction fluctuations [74], e.g., [with $x=(x,y,z)$]

with, for example, $zint(x,y)$ given by Eq. (17). Note that this expression should be bounded by $\varphi 1(x,0)\u2208[0,1]$.

*ω*) and the formal time derivative of Eq. (17) [26]

_{k}#### 3.2.5 Boundary Conditions and Computational Domain.

An important limitation of numerical simulations is manifested through the choice or requirements of boundary conditions and computational domain (or discretization of the equations). The recognition that transport in RT flows is statistically inhomogeneous in the direction parallel to the acceleration relates to the common use of periodic boundary conditions in the directions perpendicular to the acceleration. For example, when higher order accurate numerical methods are used in a code, spectral methods are used in the periodic directions and suitable high-order schemes are used in the nonperiodic direction. This strategy is helpful for achieving accuracy in the spatial discretization of the equations solved with higher computational efficiency, but is not representative of flows in nature, the laboratory, or in technological applications, e.g., the experimentally studied RT flows previously described are all confined to bounded (nonperiodic) finite rather than infinite domains. Topologically, the space defined by such two-dimensional periodic boundary conditions is mapped onto a torus (compactification), which is clearly different from a realistic geometry. The predictions of multidimensional applications of turbulent transport models using periodic boundary conditions are also subject to artifacts of this choice.

The high resolution requirements imposed for turbulent flows coupled with computational resource limitations continues to limit the size of the domain in RT flow simulations. Therefore, the computational domain size that is affordable in a simulation of a typical RT experiment is a subset of the full geometry of a tank or channel (and represents the “core” of the mixing layer) (e.g., see [67]). As mentioned in Sec. 3.2.2, this limitation of computational domain size currently precludes DNSs of RT mixing with $Sc>O(1)$.

#### 3.2.6 Transition to Turbulence and Self-Similarity.

The transition to turbulence in RT flow appears to be a continuous process. The “mixing transition” was proposed, both as an extension from stationary [75] turbulent flow and nonstationary flow [76] (see Ref. [77]). The latter also specifies a “minimum state” of turbulence. Only a few RT simulations have achieved the criteria for either mixing transition, but these simulations used simplified initial conditions and were not related to any experiment.

A fully turbulent RT flow is expected to be self-similar: arguments based on a multimode growth model were used to establish that asymptotic self-similar states in RT flow are achieved when the initial perturbation spectrum is dominated by short-wavelength perturbations, and at least 3–4 mode coupling generations have occurred [78]. This is also a stringent criterion that has not likely been satisfied in the vast majority of RT simulations and experiments performed to date, as the duration of layer growth is geometrically constrained by the boundaries in finite-sized experimental and computational domains.

#### 3.2.7 Initial Conditions.

As for numerical simulations of RT flows, the predictions of reduced-order models of these flows also strongly depend on the initial conditions chosen. This is further exacerbated by the fact that initial conditions are required for many more fields than needed for numerical simulations, e.g., turbulent fields in turbulent transport models. The predictions of these models are very sensitive to the initial conditions. As the initial conditions are not generally known for the turbulent fields, they are often treated as model parameters and adjusted to postdict some specifically chosen data.

Additional spatial information is needed for two- and three-dimensional modeling. Initial perturbations are usually characterized by their spectra in experiments, but most models require inputs as a function of space. As the number of turbulence model equations increases so does the number of initial conditions. Sufficiently accurately estimating the initial turbulent fields remains challenging, as detailed data is needed to validate initial conditions models. A recent study that uses experimental validation data has begun to address this issue [79].

#### 3.2.8 Closure Modeling.

The application of LES (see Sec. 4.1.2.4) using explicit subgrid-scale modeling to RT experiments has thus far been limited to a few studies. More research is needed to investigate the effects of buoyancy on the subgrid-scale dynamics of RT flows, which may suggest potential improvements to these models. Similarly, closure models for turbulent transport models of RT flows [80] are much less developed than for traditional shear-driven flows. Many of the closure models are simple adaptations or extensions of existing incompressible and compressible flow models, and do not self-consistently incorporate all of the physics peculiar to this class of buoyancy-driven flows. Most closure models are based on single-species, rather than multispecies flows, and many models omit important terms that have been included in turbulent combustion and reacting flow [81,82] and geophysical flow [83,84] modeling.

Model coefficients are usually set in two ways: (1) by application of a model to various cases and adjusting the coefficients associated with buoyancy production and other mechanisms until good agreement with data is obtained on average (e.g., [85,86]), and; (2) by using self-similarity constraints to calibrate the model coefficients (e.g., [57,59]). The former approach can be time-consuming and is generally inconsistent with self-similarity (which is not necessarily a disadvantage). The latter approach is expedient, but may limit the models to predicting large Reynolds number behavior which may not be realized in a particular application. The self-similarity approach is based on a specific form of the one-dimensional spatial profiles of the turbulent fields, and may drive the numerical solutions of the models toward such profiles (which may not agree with those computed from experimental or numerical simulation data).

## 4 Numerical Simulations of Rayleigh–Taylor Flow and Mixing Experiments

Many of the early simulations of Rayleigh–Taylor instability and mixing were motivated by observations and measurements in tank experiments and were (monotone) implicit large-eddy simulations [(M)ILES]. Later, the advent of DNS at higher resolutions allowed for more detailed simulations and comparisons to data. Several LES studies of RT experiments were subsequently performed. Relatively few simulations of experiments have been performed compared to simulations of “generic” configurations, i.e., using assumed simplified initial conditions that are not related to any particular experiment.

### 4.1 Numerical Approaches.

#### 4.1.1 Governing Equations.

The governing equations in the variable-density incompressible and compressible approximation used in simulations of miscible RT mixing experiments are briefly summarized here (summation over repeated indices is implied) (see Ref. [88]). In binary mixing the sum of the mass fractions is $m1+m2=1$.

##### 4.1.1.1 Incompressible variable-density approximation.

*i*,

*j*, and

*k*are vector indices corresponding to Cartesian coordinates, $\sigma ij=\mu [\u2202vi/\u2202xj+\u2202vj/\u2202xi\u2212(2/3)\u2009\delta ij\u2009\u2202vk/\u2202xk]$ is the viscous stress tensor (ignoring bulk viscosity),

*D*is the binary diffusion coefficient, and $Jr,j=\u2212\rho \u2009D\u2009\u2202mr/\u2202xj$ is the molecular diffusion flux of fluid

*r*. Summing Eq. (24) over

*r*and using the fact that the mass fractions sum to unity recovers Eq. (22). The flow is incompressible in the sense that the internal energy equation is decoupled. Volume fractions may be computed from the mass fractions as $fr=mr\u2009\rho /\rho r$. It can be shown that the velocity divergence is [89,90]

and is a consequence of the diffusive mixing of two fluids with different densities.

##### 4.1.1.2 Compressible approximation.

*χ*, the enthalpy of fluid

*r*is $hr=Ur+pr/\rho r$, the mixture molecular transport coefficients ($\varphi =\mu ,\chi ,D$ and $m\u22480.7$) are

the mixture ratio of specific heats is $\gamma =cp/cv$ ($cp=cp1\u2009m1+cp2m2$ and $cv=cp\u2212Ru/M=cv1m1+cv2m2$) and the mixture pressure is $p=\rho $$Ru\u2009T/M=(\gamma \u22121)\rho \u2009U$ with mixture molecular weight $M=(m1/M1+m2/M2)\u22121$ and universal gas constant *R _{u}*. Other definitions of mixture quantities can also be used. Most compressible RT simulations have used the nonconservative internal energy formulation without the last term in Eq. (26). The initial sound speed is chosen large enough to eliminate Mach number effects.

##### 4.1.2 Simulation Approaches.

The simulation approaches used for RT flows are briefly summarized here. Implicit large-eddy simulation, direct numerical simulation, and large-eddy simulation are briefly described. These approaches are discussed in chronological order of their general use as historically applied to RT flows.

##### 4.1.2.1 Implicit Large-eddy simulation.

An implicit large-eddy simulation (ILES) utilizes the inherent numerical dissipation/diffusion in the algorithm to provide an effective “subgrid” dissipation/diffusion at the smallest scales [93,94]. This approach is based on the assumption that the inherent dissipation/diffusion well-approximates the physical dissipation/diffusion and that the Schmidt number is $O(1)$. When low-order finite volume methods are used, the numerical truncation error has a diffusive effect and it may be possible to estimate the artificial viscosity using a modified equation analysis [95–97]. This numerical diffusion enforces the monotonicity-preserving property of the numerical algorithm, suppressing unphysical overshoots and undershoots (arising, for example, from material discontinuities such as a sharp density interface). When the implicit filter is strongly biased toward damping high wavenumber modes, the method is a monotone integrated large-eddy simulation (MILES) [98].

In (M)ILES, the fluid dynamics equations are solved without physical viscosity, diffusivity, and conductivity, i.e., the terms with *σ _{ij}*, $\Phi j$, and $Jr,j$ are not included. The small-scale “mixing” is numerically generated, rather than physically generated by viscous, diffusive, and conductive terms. The applicability of (M)ILES is largely based on the assumption that the dissipation of energy is driven by a cascade of energy to smaller scales and is essentially independent of the detailed viscous and diffusive processes. (M)ILES cannot describe the transition to turbulence, and applies to flows that have exceeded the “mixing transition” [75]. Early comparisons of (M)ILES code results with available RT experimental data indicated that (M)ILES was a reasonably accurate and computationally efficient approach.

##### 4.1.2.3 Direct numerical simulation.

A DNS resolves all dynamically relevant spatial and temporal scales in the velocity and scalar fields (i.e., up to the Kolmogorov and Corrsin–Obukhov scales, respectively [99]), limiting DNS to small-to-moderate Reynolds number, transitional RT flows. The computational expense of DNS typically leads to a single realization of a high-resolution case rather than an ensemble of simulations. Velocity and density fluctuations, and therefore statistics, can be computed using DNS data. Turbulence statistics from statistically stationary channel experiments, for example, are computed using ensemble-averaged data: the statistical quantities computed from experimental data are thus smoother than those computed from DNS data. The incompressible variable-density or compressible Navier–Stokes equations summarized in Sec. 4.1.1 are solved.

DNS requires high-order accurate spatial and temporal discretization methods. Periodic boundary conditions are typically used in the directions perpendicular to the acceleration: (pseudo)spectral methods are used in these coordinate directions and a high-order finite difference (e.g., compact difference) method is used in the direction parallel to the acceleration. High order accurate boundary conditions must be implemented to maintain overall accuracy. The progressive increase in computational capabilities has enabled well-resolved simulations of turbulent RT flows for moderate Reynolds numbers in the last two decades. In addition to resolving the Kolmogorov scale and the Corrsin–Obukhov scale, the initial diffuse interface thickness and the dominant initial interface perturbation amplitude must be resolved.

##### 4.1.2.4 Large-eddy simulation.

A LES resolves the largest (filtered) scales and models the effects of the smaller (unresolved) scales on the large scales using an explicit subgrid-scale model [100]. Although this approach may require less grid resolution than DNS, it is still computationally expensive. Compared to (M)ILES and DNS, relatively few LESs of RT flows have been performed to date. The filtered equations are solved (usually, the filtering is implicit rather than explicit). Applying the filter to the equations gives the filtered (resolved-scale) equations with unresolved quantities modeled using the resolved fields.

The filtered field has general form [100,101] $\varphi \xaf(x,t)=\u222bG\Delta (x,x\u2032)\varphi (x\u2032,t)d3x\u2032$, where $G\Delta (x,x\u2032)$ is the filter kernel satisfying a number of consistency requirements and Δ is the filter width (related to the grid spacing). For a compressible flow, $\varphi \u0303(x,t)=[1/\rho \xaf(x,t)]\u222bG\Delta (x,x\u2032)\rho (x\u2032,t)\varphi (x\u2032,t)d3x\u2032$. LES typically requires the same low dissipation/diffusion numerical methods as DNS. The subgrid-scale models used for RT simulations have largely been adapted from those used for LES of homogeneous turbulence.

### 4.2 Simulations of Accelerated Tank Experiments.

Numerical simulations of accelerated tank RT experiments are discussed here.

#### 4.2.1 Simulations of Rocket-Rig Experiments.

One of the earliest applications of MILES to three-dimensional multimode RT instability growth used the Lagrange-remap TURMOIL code (also called TURMOIL3D) [25]. RT mixing was simulated over a wide range of density ratios $\rho 1/\rho 2=1.5$, 3.0, and 20.0 corresponding to At* *=* *1/5, 1/2, and 19/21: the two smaller density ratio cases were motivated by rocket-rig experiments using water/pentane and NaI solution/pentane, respectively [29], and the third case was motivated by experiments using pentane/SF_{6} [28]. The grid resolutions for the At* *=* *1/5 case were 80^{3} and 160^{3}, and for the At* *=* *1/2 and 19/21 cases were $120\xd7802$ and $240\xd71602$. Several quantities were reported including growth parameters, molecular mixing parameter, ratio of the total energy dissipation to the potential energy released etc. Figure 1 shows perspective views of the bubbles and spikes visualized using isosurfaces of the heavy volume fraction for the largest At case. At the latest times $\alpha bY\u22480.04$ and $\alpha bss\u22480.054$ from the finest grid simulations, which underestimated the experimental value $\alpha bY\u223c0.06$. It was hypothesized that insufficient grid resolution led to growth that did not become self-similar, and that the simulations were miscible while nearly all of the experiments were immiscible. A subsequent $270\xd71602$ simulation with At* *=* *0.5 gave an even smaller value $\alpha bY\u22480.03$ [27]. This small value was possibly attributed to an incomplete loss of memory of initial conditions [102] and the numerical dissipation. The discrepancy between the numerical and experimental $\alpha bY$ values persisted for much higher resolution $720\xd76002$ simulations [60], which was attributed to small amplitude, long-wavelength initial perturbations present in the experiments but not accounted for in the simulations.

LESs of RT mixing using the FronTier front-tracking code were performed using a dynamic Smagorinsky-like model [103] that includes subgrid-scale concentration variance [104]. Similar simulations were used to demonstrate that *α* is indeterminate, as it depends on regularizing physics and on nonideal initial conditions [105]. Simulations modeling RT accelerated tank experiments [30] were also performed as a validation of the numerical method [106–108]. It is noteworthy that the initial perturbations were estimated from experimental photographs at early time in Ref. [108]. Rocket-rig experiment number 112 considering the miscible mixing of NaI solution/water (At* *=* *0.308) was simulated using DNS at resolution $96\xd75762$ [106]. The simulation gave $\alpha Y=0.055$ and the experimental value was $\alpha Y=0.052$. Figure 2 shows a comparison of the mixing layer between the simulation and experiment at late time.

#### 4.2.2 Simulations of Tilted Rocket-Rig Experiments.

An early simulation of a tilted rocket-rig experiment [30] was performed using the TURMOIL3D code in three dimensions [8]. These experiments (with a tilt angle) involved a large-scale overturning motion superposed with RT mixing. The simulation reproduced the large-scale features of the flow (see Fig. 3), but underresolved the small-scale mixing near the center of the layer. The resolution and other details of the simulation were unspecified.

The predictions of two MILES codes (TURMOIL and RTI-3D) and a DNS code (CFDNS) applied to tilted rocket-rig experiments with At* *=* *0.482 were compared [109]. The details of the three-dimensional simulations were specified to make them potentially reproducible by other researchers. One- and two-dimensional mixing metrics were used to examine the development of the mixing layer and the overall flow evolution. Comparison of the simulations with experiment indicated that the large-scale overturning, central mixing widths, and spike/bubble sidewall penetrations were well-described in all of the simulations. A comparison between MILES and DNS showed overall good agreement between metrics such as the molecular mixing. The DNS indicated a dependence on Reynolds number that should be further investigated experimentally. Figure 4 shows the (negative) density–specific volume correlation, $b=\u2212\rho \u2032(1/\rho )\u2032\xaf$ (used algebraically or differentially in a number of turbulence models) at late time.

#### 4.2.3 Simulations of Linear Electric Motor Experiments.

A comparative study of three-dimensional ($1282\xd7256$ and $2562\xd7512$) multimode RT instability evolution with $At\u22480.5$ in the strong mode-coupling limit (short-wavelength initial perturbations) was performed using the TURMOIL3D, FLASH, WP/PPM, NAV/STK, RTI-3D, HYDRA, and ALEGRA codes [26]. The HYDRA and ALEGRA codes have an interface reconstruction capability (to keep the interface sharp), and were used to perform simulations with and without interface reconstruction. Simulation results were also compared to LEM experimental data. Following an initial transient, it was found that $hb\u2248\alpha bY\u2009At\u2009g\u2009t2$ at late times, with the codes solving the dissipative equations predicting $\alpha bY=0.025\xb10.003$; the small value of $\alpha bY$ compared to experimentally inferred values was attributed to a density dilution due to small-scale mixing or entrainment, and applies whether or not interface reconstruction was used. The multimode interfacial perturbations were given by Eq. (17). Volume fraction profiles, integral mixing widths, molecular mixing profiles, and molecular mixing parameters were computed. Various quantities from the codes were compared to LEM data (e.g., volume fractions) and shown to be in good agreement, although the codes gave values of $\alpha bY$ that considerably underestimated the experimental values $\alpha b\u223c0.057\xb10.008$.

Numerical simulations ($2562\xd7512$ and 512^{3}) of the hexane/salt water (At* *=* *0.48) LEM experiments with acceleration, deceleration, and acceleration were performed using the RTI-3D code [110]. Figure 5 shows a comparison between backlit images from LEM experiments and simulations with two different initial perturbation spectra, and Fig. 6 shows the acceleration, amplitude, and aspect ratio for bubbles versus displacement. It was found that the dominant bubbles and spikes produced during the initial acceleration phase-reverse their direction during the deceleration phase and “shred”: the dominant structures are reduced in size by $\u223c50%$ and the molecular mixing increases $\u223c100%$. This results in a reduced growth rate during the final acceleration phase. Simulations were performed using initial spectra with only short-wavelength perturbations and with broadband perturbations: the latter gave good agreement with the experimental data.

### 4.3 Simulations of Stationary Tank Experiments.

Numerical simulations of stationary tank RT experiments are discussed here.

#### 4.3.1 Simulation of the Linden, Redondo, and Youngs Experiments.

MILES (90^{3} and 180^{3}) using the TURMOIL code were performed to model stationary tank experiments in which RT instability was induced following the removal of a barrier initially separating salt water and fresh water [38,52]. The density ratio $\rho 1/\rho 2=1.5$ was larger than in the experiment to ensure that the initial density difference is large compared to any small density variations within each fluid arising from a finite Mach number (the experiments had $At=1\xd710\u22124$–$5\xd710\u22122$). A random initial interfacial perturbation as a sum of Fourier modes was used. The removal of the barrier induced a long-wavelength two-dimensional perturbation, which was modeled by the addition of a single-mode perturbation with amplitude and wavelength estimated from visual observations. The experiments gave $\alpha bY=0.044\xb10.005$ and $\Theta (z,t)\u22480.75$ near the center of the layer. The simulation gave $\alpha bY\u22480.035$ and $\Theta \u0302(t)\u22480.81$ at resolution $1682\xd7230$ in good agreement with the experimental values.

The simulations gave mixing efficiencies $\eta =1\u2212\Delta \Phi /\Phi max=0.48$ (random initial perturbations) and 0.47 (with long-wavelength initial perturbation), where $\Delta \Phi $ is the final potential energy loss and $\Phi max$ is the change in potential energy if the fluids interchange positions without any mixing. These values are significantly larger than the experimental values. Large-scale overturning dynamics arising from the circulation induced by the barrier removal were observed in the experiments at late times, which may have reduced the mixing efficiency in the experiments. Concentration fluctuations were measured using conductivity probes placed at various locations in the tank. The amount of mixed product was quantified using the color change of a *p*H indicator to determine the extent of a chemical reaction. There was good agreement between the experimental and simulation results.

#### 4.3.2 Simulation of the Dalziel, Linden, and Youngs Experiments.

The self-similarity and internal structure of RT mixing with $At\u22482\xd710\u22123$, in which a salt water layer was initially separated from a fresh water layer by a barrier that is removed to generate mixing, was studied experimentally and numerically using the TURMOIL3D code [39,52]. As the barrier was removed, a vortex sheet perturbation was generated. Discrepancies between previous experimental and numerical studies [38] were attributed to incomplete modeling of the experimental conditions, and it was suggested that the spatial structure of the initial conditions is more important than the amplitude of the perturbations in determining the transition to self-similar mixing. An initial streamfunction was used to model the vortical perturbation arising from the barrier removal. The density ratio $\rho 1/\rho 2=1.2$ ensured that small density fluctuations due to compressibility had little effect, and small enough that the Boussinesq approximation was valid. The numerical Schmidt number was $\u223c1$ while the molecular Schmidt number was $Sc\u223c103$; the Reynolds number in the experiment was assumed large enough that small-scale mixing was insensitive to the Schmidt number.

The simulation resolution was $160\xd780\xd7200$: “idealized” and “real” initial conditions were used. The idealized conditions had $v(x,0)=0$ and an interfacial perturbation consisting of a sum of Fourier modes with randomly chosen amplitudes. The real initial conditions had $vx(x,0)$ and $vz(x,0)$ generated by the streamfunction and used the same random interfacial perturbations as in the ideal case. The initial conditions were chosen such that the memory of the small-scale initial conditions was lost in the simulations. An ensemble of simulations was performed, and ensemble-averaged statistics were reported. The bubble and spike mixing widths were obtained by thresholding the mean concentration profiles. Other integral measures using the mean concentration profile were also used. An unambiguous value of *α _{b}* could not be determined from the experiments or simulations. The concentration variance spectrum was consistent with a Kolmogorov $k\u22125/3$ power-law. The spectra from the simulations with ideal and real initial conditions were fit with an exponent −1.79 and −1.63, respectively. Concentration probability density functions, and local and global molecular mixing parameters were also computed. Figure 7 shows a comparison at two late times between a typical experiment (LIF images) (top), simulations using idealized initial conditions (middle), and experimental initial conditions (bottom).

#### 4.3.3 Simulations of the Lawrie and Dalziel Experiments.

A density diffusion equation was developed to describe RT mixing in a tall, narrow (i.e., very large aspect ratio) tube arising from initially homogeneous layers above and below the unstable interface [111]. The mixing length hypothesis was used to obtain a variable turbulent diffusivity and a nonlinear diffusion equation. A new $h(t)\u221dt2/5$ similarity solution to the diffusion equation in an infinite domain was obtained (differing from the classical *t*^{2} RT scaling due to the presence of an additional lengthscale–the tank width–in the similarity analysis). The equation was integrated numerically to match the boundary conditions of an At* *=* *0.0015 experiment corresponding to a finite domain (a tall tube), and the results were compared with $642\xd72560$ ILES using the MOBILE code. Figure 8 shows the self-similar collapse of concentration profiles from an overturned tank experiment with large aspect ratio [112] and ILES.

The TURMOIL code was used to simulate RT mixing with At* *=* *0.091 [113] using density profiles (rather than uniform densities in the upper and lower fluids) based on stratified RT experiments [114,115]. The power spectrum of the initial interfacial perturbations was $P(k)\u221dk\u22123$. The grid resolution was not specified. It was found that the late-time mixing profiles are similar to the experimental results for similar initial conditions; a range of additional initial conditions was also considered to investigate the sensitivity of the results. Comparison of a model for the late-time mixing layer structure based on the maximization of configurational entropy with the simulation results indicated good agreement. Figure 9 shows the passive scalar 0.1% and 99.9% isosurfaces at early and late times from the simulation.

### 4.4 Simulations of Channel Experiments.

The numerical simulations of water and gas channel RT experiments are discussed here.

#### 4.4.1 Implicit Large-Eddy Simulations of Water Channel Experiment.

Three-dimensional $128\xd7256\xd7128$ MILES using the RTI-3D code and initialized using experimental water channel velocity and density data ($At=7.5\xd710\u22124$) [44] showed that simulations using nonzero initial velocity yielded an *α _{b}* in better agreement with late-time experimental data than simulations using zero initial velocity field [74]. The two-dimensional density perturbations were expressed as in Eq. (17). The amplitudes of the longest modes obtained from the Fourier amplitudes of the measured density fluctuations at

*x*=

*0.1 cm were used for*

*a*,

_{mn}*b*,

_{mn}*c*, and

_{mn}*d*.

_{mn}The azimuthally averaged kinetic energy spectra used in the simulations had the same energy for the included wavenumbers as in the experiment: $\u222bkminkmax\u27e8Esim(k)\u27e9\theta dk=\u222bkminkmaxE\u2009exp\u2009(k)dk$. The density surface perturbations were converted to volume fraction fluctuations using Eq. (18) (with $\varphi \u2192f$), and initial velocity perturbations were obtained from the gradient of the potential (20). PIV data were used to obtain the velocity data, and the *a _{mn}* were chosen so that the integral of velocity spectra from the simulations and experiment agreed at the measurement point. Three simulations using only density perturbations gave values $\alpha bcl\u223c0.06$ [based on the centerline value of $v\u2032/(2\u2009At\u2009g\u2009t)$, where $v\u2032$ is the vertical velocity fluctuation], somewhat smaller than measured experimentally, and also lagged behind the experiments in time. It was hypothesized that this was due to the presence of velocity perturbations in the experiment. A fourth simulation was performed using velocity perturbations which showed better agreement with the experimental

*α*.

_{b}#### 4.4.2 Direct Numerical Simulation of Water Channel Experiment.

A $1152\xd7760\xd71280$ DNS using initial conditions, geometry, and physical parameters approximating those of a transitional water channel RT mixing experiment with $At=7.5\xd710\u22124$ [45,52] was performed using the MIRANDA code [67]. Temperature diffusion with Prandtl number Pr* *=* *7 was modeled by mass diffusion with an equivalent $Sc=7$. The density and velocity fluctuations measured just off of the splitter plate in this experiment were parameterized to provide realistic, anisotropic initial conditions for the DNS.

Large-scale quantities such as *h _{b}* and the mixing layer growth parameter

*α*(see Figs. 10 and 11), statistics such as the centerplane molecular mixing parameter, and vertical velocity and density variance spectra from the DNS were in favorable agreement with the data. The DNS predicted $\alpha bY,\u2009\alpha bRC\u22480.07$ at the latest time, in excellent agreement with the experiment. The late-time value $\Theta \u0302(t)\u22480.55$ also compared favorably with the experimental value $\u22480.60$. Density and vertical velocity variance spectra were also in excellent agreement between the DNS and experiment. At the latest time, the mixing layer Reynolds number attained was $Reh(t)\u223c1710$. Other quantities were also obtained from the DNS such as the integral- and Taylor-scale Reynolds numbers, Reynolds stress and dissipation anisotropy, hypothetical chemical product formation measures, other local and global mixing parameters, and the statistical composition of mixed fluid [116].

_{b}#### 4.4.3 Large-Eddy Simulations of Water and Gas Channel Experiments.

LESs of RT water [45,46,117] and gas channel [47] experiments using a dynamic Smagorinsky-like model in the FronTier code were performed as a validation of the numerical method [106–108]. Simulations ($288\xd7180\xd7240$) of the water channel experiments at $Sc=7$ and 620 were performed using subgrid-scale mass diffusion and front-tracking. The computational domain and initial conditions were those in Ref. [67]. Excellent agreement with quantities such as $\alpha b$ measured in these experiments spanning a range of $At$ and Sc was reported. The simulations and experiments gave for: (1) the Sc* *=* *0.7, At* *=* *0.035 Banerjee–Andrews air/He gas channel experiment, $\alpha Y=0.069$ and $\alpha expY=0.065$–0.07; (2) the Sc* *=* *7 water channel experiment, $\alpha Y=0.075$ and $\alpha expY=0.077\xb10.011$, and; (3) the Sc* *=* *620 water channel experiment, $\alpha Y=0.084$ and $\alpha expY=0.085\xb10.005$. Comparing simulations of rocket-rig experiment 112 with the water channel simulations indicated that the former had no significant long-wavelength noise, less viscous dissipation, and larger initial mass diffusion compared to the latter.

## 5 Modeling of Rayleigh–Taylor Mixing Experiments

This section discusses progress on modeling RT flow and mixing experiments. Statistical modeling of RT flows encompass models for the mixing layer properties, modal growth models, bubble merger models, buoyancy–drag models, and turbulent transport models. Progress on various aspects of theoretical RT modeling, as well as open questions, were reviewed in Ref. [118].

### 5.1 Modeling Approaches.

Several approaches that have been used for the modeling of RT flow and mixing are described here: modal models, buoyancy–drag models, and turbulent transport models.

#### 5.1.1 Modal Models.

Modal growth models were developed to predict the transition from early time linear growth to weakly nonlinear growth of multimode perturbations. Several models were developed to predict the bubble/spike front asymmetry with increasing Atwood number.

**k**(and later

**p**) is the two-dimensional wavevector into the interface evolution and Bernoulli equations, and using a Taylor expansion of $\u2202\phi r/\u2202z$ truncated at first-order gives the Fourier amplitude evolution equation $\u2202zk/\u2202t=k\u2009\phi k1(t)\u2212\u2211pk\xb7p\u2009\phi p1(t)\u2009zk\u2212p(t)+\cdots =\u2212k\u2009\phi k2(t)\u2212\u2211pk\xb7p\u2009\phi p2(t)\u2009zk\u2212p(t)+\cdots $ Substituting $\phi k1,2=\xb1(1/k)\u2202zk/\u2202t+\u2211p[k\xb7p/(kp)](\u2202zp/\u2202t)\u2009zk\u2212p+\cdots $ and their time derivatives into the Fourier-transformed Bernoulli equation gives the second-order amplitude evolution equation (with $k\u0302=k/|k|,\u2009\xi =1/2\u2212k\u0302\xb7p\u0302\u2212(k\u0302\u2212p\u0302)\xb7p\u0302/2$, and $\zeta =1\u2212k\u0302\xb7p\u0302$)

where $\omega k=g\u2009At\u2009k$ is the RT instability growth rate.

with $Wk,p(t)=G+{cosh[(\omega p+\omega |k\u2212p|)t]\u2212cosh(\omega k\u2009t)}+G\u2212{cosh[(\omega p\u2212\omega |k\u2212p|)t]\u2212cosh(\omega k\u2009t)}$ and $G\xb1\u2261F\xb1/[(\omega p\xb1\omega |k\u2212p|)2\u2212\omega k2]$. The perturbation expansion converges provided that the mode amplitudes satisfy $zk(t)<0.1\u2009\lambda $, where $\lambda =2\pi /k$ is the perturbation wavelength. Each mode grows in two distinct stages: (1) a linear regime, in which RT instability leads to exponential growth with the final amplitude proportional to the initial amplitude and; (2) an asymptotic constant velocity regime attained when a mode reaches a nonlinear threshold.

#### 5.1.2 Buoyancy–Drag Models.

Buoyancy–drag models are phenomenological coupled ordinary differential equation models in time for the bubble and spike mixing widths (see Ref. [121] for a review). These models are based on applying Newton's second law to bubble and spike structures. The models can describe two- or three-dimensional bubble fronts, and can incorporate early time behavior consistent with linear theory and late-time behavior consistent with self-similarity. Such models are not limited to small Atwood number flows. Time-dependent accelerations and densities can be implemented. Several buoyancy–drag models were developed to describe RT experiments with time-dependent accelerations and predict the bubble/spike front asymmetry with increasing Atwood number. See Ref. [122] for a recent survey of growth parameters for RT mixing obtained from buoyancy–drag models.

*h*can also be used. With bubbles and spikes corresponding to

*r*=

*1 and*

*r*=

*2, respectively, the bubble and spike equations can be combined into the single generic equation [122–126]*

*b*=

*2/3 and 1 in three and two dimensions, respectively. The velocity is $v3\u2212r(t)=dh3\u2212r/dt$. The factor $Rr=\rho r/(\rho r+\rho 3\u2212r)$ introduces a density dependence in the drag terms, and $[hs(t)/hb(t)]\beta $ increases the drag on the spikes as the spike over bubble width increases [126].*

#### 5.1.3 Turbulent Transport Models.

Turbulent transport models utilize statistical averaging of the fluid dynamics equations to define mean flow equations that contain unclosed terms arising from nonlinearity [66]. The fields are written as a sum of a mean and fluctuating field $\varphi =\varphi \xaf+\varphi \u2032$ and $\varphi =\varphi \u0303+\varphi \u2033$ for incompressible and compressible fields, respectively: $\varphi \xaf$ and $\varphi \u0303=\rho \varphi \xaf/\rho \xaf$ are the Reynolds- and Favre-averaged fields, respectively. First-order closure models are most often used to obtain expressions for these terms and require a turbulent viscosity. Various increasingly complex approaches are available to obtain the turbulent viscosity (see Ref. [121] for a review), ranging from one-, two-, three-, and four-equation turbulence models. Second-order Reynolds stress closure models have also been developed. Several turbulent transport models were applied to RT experiments. Historically, early model applications were to water channel experiments; later applications focused on the tilted rocket-rig experiments.

with *σ _{m}* a dimensionless turbulent Schmidt number,

*K*the turbulent kinetic energy, $\mu t=\rho \xaf\u2009\nu t$, and $\nu t$ a turbulent viscosity to be determined by the solution of modeled turbulent transport equations for

*K*and another mechanical turbulent field (the turbulent kinetic energy dissipation rate

*ε*or turbulent lengthscale

*L*). Second-order closure models model

*τ*and $\Gamma j,r$ using transport equations rather than using algebraic closures. For incompressible flows see Refs. [128] and [129] and for compressible flows see Refs. [130] and [131] for additional background and details.

_{ij}*K*–

*ε*model $\nu t=C\mu \u2009K2/\epsilon $ is computed using the modeled transport equations (neglecting the pressure–dilatation and dilatation dissipation) [132,133]

where $\sigma \rho $ is a dimensionless coefficient. If the flow is incompressible, the mean pressure gradient term (arising from a neutrally stable adiabatic stratification) vanishes as the sound speed $cs=\gamma p\xaf/\rho \xaf$ is formally infinite.

### 5.2 Modeling of Accelerated Tank Experiments.

A variety of analytical, modal, bubble merger, buoyancy–drag, and turbulent transport models have been developed and applied to accelerated tank RT experiments as discussed here.

#### 5.2.1 Analytical Modeling of Linear Electric Motor Experiments.

A model for the density-ratio-invariant mean mass and volume fraction profiles, $m\xaf$ and $f\xaf$, of classical incompressible RT mixing was derived [134]. Asymptotic analysis was used to derive density-ratio-invariant hybrid profiles $C(X)=Fm\xaf+(1\u2212F)f\xaf$ with $At=(R\u22121)/(R+1)$ and $F=(1+At)(1\u2212X)/2$ and obtain an analytical expression for *C*(*X*), where $X\u2208[0,1]$ is a mixing-width-rescaled coordinate. The invariance was validated by comparing model predictions to experimental (rocket-rig and LEM) and DNS data. Combining the *R* invariance of *C*(*X*) and the relation $f\xaf=m\xaf/[m\xaf+(1\u2212m\xaf)R]$ derived from the definitions gives the mean profiles of mass fraction and volume fraction at different density ratios. Figure 12 shows a comparison of *C*(*X*) from the model and experimental and simulation data.

#### 5.2.2 Modal Modeling of Rocket-Rig and Linear Electric Motor Experiments.

A two-dimensional multimode model of incompressible RT instability that extends the Haan [119] model and describes low-order nonlinear mode coupling, mode growth saturation, and mode coupling at long times was developed [135]. The predictions of the model were compared to those of the two-dimensional, finite difference ALE interface tracking code LEEOR2D at At* *=* *0. An initially flat interface was used with velocity perturbations (a “band-type” initial condition was also evolved). Figure 13 shows the spectrum of a rocket-rig case as predicted by the model and compared to the predictions of the Haan model.

with $\eta 0=\u2329h0k\u232a/\lambda 0=k\u2329h0k\u232a/(2\pi )$; Fr =* *0.82 provided the best fit to the linear electric motor experimental results. In the application of the model to the experiments, early time photographs were used to estimate the initial perturbations. Figure 14 shows the variation of *α _{b}* and

*β*with scaled rms initial amplitude for Fr

_{b}*=*

*0.49 and 0.82 from the models and numerical simulations.*

*ξ*, modal amplitude in the interface-normal direction

*η*, and vorticity amplitude

*ω*(integrated across the interface)

where *k* is the wavenumber, $H(\xb7)$ is the Hilbert transform, and $\epsilon a$ is a linear artificial viscosity that regularizes large wavenumber instabilities associated with the pseudo-spectral method used to numerically solve the equations. The quantity $||(1+k\xi ,k\eta )||$ is the slope of the interface and $(a,b)\u2261ae1+be2$.

was used to show that self-similar states are achieved when the initial perturbation spectrum is dominated by short-wavelength perturbations, and at least 3–4 mode coupling generations have occurred [78]. The mean field approximation allows the sums to be replaced by their root-mean-squares $\u2211pSk,p\u2192\u2329S\u232ak=\u2211pSk,p2$, so that the solutions of the resulting equations represent smooth ensemble-averages for all realizations having the same initial perturbation spectrum but different phases. Saturation is treated as in the Haan model [119]. At saturation, the growth of a mode with wavenumber **k** transitions from the value given by Eq. (51) to the asymptotic values $dzk/dt\u221d\lambda /t$. Saturated modes do not participate in mode coupling and stop growing once their amplitude reaches a few $Sk$. The bubble front width is $hb(t)=2\u2211kzk2$. Figure 15 shows a comparison of model results for miscible and immiscible mixing, miscible simulations, and LEM experimental data.

#### 5.2.3 Bubble Merger Modeling of Rocket-Rig Experiments.

A three-dimensional bubble merger model was developed using renormalization group methods [125,137]. The model predicts the bubble growth parameter $\alpha b=0.0635$–0.0742, the mean bubble radius, and the bubble height separation at merger time. Comparisons to rocket-rig experimental data (experiments 104, 105, and 114) [30] were performed, and indicated good agreement. The model consists of three essential components: the velocity of a single bubble in the unstable regime, a mixing envelope velocity, and a merger model. The bubble radii variance at fixed time is a model parameter.

#### 5.2.4 Buoyancy–Drag Modeling of Rocket-Rig and Linear Electric Motor Experiments.

*h*(

*t*) is consistent with Eq. (1), where

*α*depends on the kinetic energy dissipation rate. The expression for the kinetic energy is based on a wavelength renormalization hypothesis similar to that proposed in Ref. [27]. The kinetic energy per unit area is

*h*and

*v*is $(d/dt)(\u2202T/\u2202v)\u2212\u2202T/\u2202h=Q$, where $Q=\u2212c(\rho 1+\rho 2)\u2009v|v|/2$ is a dissipative (drag) force. According to the wavelength renormalization hypothesis, the single-mode perturbation wavelength

*λ*is replaced by $\lambda (t)=b|h(t)|$, with $b\u223cO(1)$ [27,139–142]; n.b.,

*b*is not the geometrical factor introduced previously in Sec. 5.1.2. For RT instability

with solution Eq. (1) and $\alpha =\pi /[2(b+2\pi \u2009c)]$. The transition between the linear and nonlinear regimes may be achieved with $\lambda (t)=max[\lambda 0,b|h|+(1\u2212m\u2009b)\lambda 0]$, where *m* is a parameter. Figure 16 shows a comparison of model predictions with experimental data for four different acceleration histories.

*g*. For the bubbles,

*α*

_{2}is insensitive to At, and for the spikes,

*α*

_{1}increases as a power law with the density ratio. The density at the initial interface location $\rho 0=(\rho 1h1+\rho 2h2)/(h1+h2)$ is obtained from mass conservation $(\rho 2\u2212\rho 0)h2/2=(\rho 0\u2212\rho 1)h1/2$. The average mixture densities $\rho M,r=(\rho r+\rho 0)/2$ can be expressed as

where $R=\rho 1/\rho 2$ and $\chi =h1/h2$. Figure 17 shows $\alpha s/\alpha b$ (i.e., $\alpha 1/\alpha 2$) predicted by the model as a function of density ratio *R* with *C *=* *1, 2, and 4 (solid lines). The points are measurements on the LEM.

A generalized buoyancy–drag model was developed by incorporating conservation of mass and momentum and a symmetry factor to account for density variations [143]. A parabolic velocity profile was assumed, and the local rate of advance of the mixing fronts was approximated by a volume-averaged rate. Three model parameters were determined through the use of asymptotic solutions. The model was validated by comparison to data from several experiments with a large range of density ratios and different acceleration histories, including the constant acceleration rocket-rig experiments [31,61] and with variable acceleration experiments [32,33]. The predictions of this model were in better agreement with data than those of the Dimonte [126] model. Figure 18 shows a validation of the model for constant accelerations and at all density ratios.

#### 5.2.5 Turbulent Transport Modeling of Rocket-Rig and Linear Electric Motor Experiments.

A compressible two-dimensional multiphase model coupled to a *K*–*L* turbulence model was developed [31] (see Ref. [144]) and included differential acceleration induced by the pressure gradient between the fluids, drag proportional to the velocity difference squared, mass flux, and increase of bubble or spike radius proportional to the mixing layer width. Interphase drag is a source of turbulent kinetic energy, and was needed to describe late-time mixing by turbulent diffusion observed in rocket-rig acceleration/deceleration experiments. Figure 19 shows a comparison between mixing widths from an acceleration/deceleration experiment and the model.

A two-scale *K*–*ε* model coupled to a buoyancy–drag model was developed for RT mixing [145]. The buoyancy–drag equation describes the large-scale dynamics, and an energy equation corresponding to the large-scale dynamics is added to the *K* and *ε* equations. The model was applied to a LEM SF_{6}/freon experiment [33] and to an At* *=* *1/2 experiment with large acceleration [146], and the predicted mixing layer widths were in good agreement with the data.

A compressible *K*–*L* model was developed for RT turbulent mixing using gradient-diffusion closures and production and dissipation terms motivated by buoyancy–drag modeling [57]. The model coefficients were calibrated using self-similar expressions for the growth parameters and exponents with specified values. Comparisons of the model predictions using the two-dimensional CALE code to RT experimental data were favorable. Figure 20 shows $\alpha bss$ and $\alpha sss$ as a function of Atwood number computed from the model and compared to LEM data, and Fig. 21 shows a comparison of the volume fraction profiles across the mixing layer from the model and LEM experiments with At* *=* *0.2 and 0.32.

Tilted rocket-rig experiments were also used to validate the four-equation BHR-2 model in two dimensions in the arbitrary Lagrangian–Eulerian (ALE) FLAG code [147]. Model predictions at grid resolutions 152 × 240, 232 × 360, and 304 × 480 were compared to $3602\xd7480$ and $5122\xd7768$ ILES using the RTI-3D code. The quantities compared were visualizations of $4f1f2$, contours of *f*_{1}, contours of *K*, *a _{z}*, and

*b*, and integral quantities including the bubble and spike widths, mixing layer width, tilt angle, and spatially integrated turbulent kinetic energy. Figure 22 shows the contours of heavy-fluid volume fraction from the model and simulation. The energy budget was also computed as a function of time.

A comparison of the predictions of the *K*–*L*, *K*–*L*–*a _{i}*,

*K*–

*L*–

*a*–

_{i}*b*, and BHR-2 models applied to one-dimensional RT mixing and to the two-dimensional tilted rocket-rig case was presented [148]. The models were implemented in the CNS3D code, and results were also compared to those from the TURMOIL code. Figure 23 shows the integral mixing width evolution with scaled time for the models and ILES. The main conclusions were: (1) the more complex tilted-rig case was necessary to identify differences between the models; (2) the

*K*–

*L*–

*a*model is more accurate than the

_{i}*K*–

*L*model and uses a simpler form of the turbulent kinetic energy production term; (3) using a modeled

*b*transport equation somewhat improves the accuracy compared to the algebraic expression used in the

*K*–

*L*–

*a*model, particularly in the prediction of

_{i}*a*and $K$ in the tilted-rig case; (4) with respect to the tilted-rig case, the BHR-2 model performs better than the

_{i}*K*–

*L*–

*a*–

_{i}*b*model on the spike side but worse on the bubble side, and; (5) overall, the BHR-2 model provided the closest predictions to ILES in the bubble and spike regions.

### 5.3 Modeling of Stationary Tank Experiments.

Relatively little modeling of stationary (nonaccelerated) RT tank experiments has been performed and is discussed here.

#### 5.3.1 Two-Fluid Modeling of Overturning Tank Experiments.

where *f ^{r}*,

*ρ*, and $vir$ are the volume fraction, density, and velocity of fluid

^{r}*r*,

*g*is the constant gravitational acceleration,

_{i}*K*is the interphase drag, and $f1+f2=1$. Models with various algebraic and differential lengthscales (mixing lengths) entering the expression for

_{d}*K*were used to simulate overturning tank RT experiments without and with a tilted interface in one and two dimensions, respectively [37]. As the experiments were performed using a narrow tank, the behavior of the flows was largely two-dimensional for tilted interfaces. It was shown that the two-fluid model, coupled with a novel transport equation for the lengthscale, gave good predictions of the flow morphology as compared to the experimental data. In contrast to the Youngs model [31], the models evaluated by Andrews [37] did not account for turbulent diffusion effects; however, the lengthscale equation (with some modifications) was used in the Youngs [31] model.

_{d}#### 5.3.2 Turbulent Transport Modeling of Tank Experiments With Energy Deposition.

*K*–

*ε*model in which the gradients of turbulent kinetic energy and its dissipation rate are zero was proposed:

where *h* is the mixing layer width. Experimental data indicated that $h\u223ct3$ at late time, and self-similarity gives $K\u223c(dh/dt)2\u223ct4$ and $\epsilon \u223cK/t\u223ct3$ consistent with this simple model. Figure 24 shows a comparison between the experimental and theoretical mixing layer width as discussed here.

### 5.4 Modeling of Channel Experiments.

A variety of analytical and turbulent transport models have been developed and applied to water and gas channel RT experiments as discussed here.

#### 5.4.1 Analytical Modeling of Water and Gas Channel Experiments.

A rapid acceleration model was developed to describe RT turbulent mixing layer growth subject to a very strong acceleration [151]. This model encompasses linear rapid distortion theory applied to unstably stratified flows and an equation for the mixing layer width. The nonlinear mechanism arising from the interaction between the turbulent quantities and the mean concentration field leads to a self-similar regime. Comparisons were also made to a $5122\xd71024$ DNS using the TRICLADE code at At* *=* *0.1 using an initial annular spectrum. Figure 25 shows the relationship between *α* and $\Theta \u0302$ from different experiments and numerical simulations (point data) compared to the formula for *α* that was derived.

Self-similar solutions of asymmetric RT mixing were obtained using a piecewise step function approximation where the main constraint imposed is mass conservation [152]. This model was used to predict the mean density within the layer and the layer asymmetry for a given value of At. The model predictions were compared to gas channel [48] and LEM [32] experimental data and to simulation results, confirming the predictions for the asymmetry of the *α _{b}* and

*α*values. The model leads to a simple correction to the formulation of the mixing layer growth, which is consistent with

_{s}*α*for a given system that is independent of the density difference for both immiscible fluids and miscible fluids with small mass diffusion. Additionally, the model can be used to compute the energy released by the instability. Figure 26 shows experimental results for

*α*and

_{b}*α*as a function of At, together with the theoretical curves predicted by the model.

_{s}#### 5.4.2 Turbulent Transport Modeling of Water and Gas Channel Experiments.

One of the earliest applications of a turbulence model to RT experimental data considered thermally stratified flow in a water channel. This motivated the development of a *K*–*ε* and spectral–physical space model. Later, three- and four-equation models were developed, followed by second-moment closure models.

*K*–

*ε*models were applied to small Atwood number RT mixing [153–155]. A model formulation using the temperature was applied to water channel experiments with and without shear [155]. Figure 27 shows the calculated mixing layer width compared to water channel experimental data, indicating very good agreement until wall effects begin to affect the experimental measurements. A self-similar analysis, in which inverted parabolic profiles for

*K*and

*ε*and a linear profile for the mean temperature were substituted into the

*K*and

*ε*equations, was used to derive the turbulent fields and the bubble growth parameter as a function of the model coefficients

Comparison of the model predictions with the experimental data was used to estimate $C\epsilon 0\u22480.9$ and $\sigma T=0.6$ (with the other coefficients having their standard values); an equivalent formulation using the mean density gives the same result with $\sigma T\u2192\sigma \rho $ [154].

Several spectral–physical models were subsequently developed for RT flow. The BHRZ model was developed, in which transport equations for the two-point, single-time correlations describing variable-density incompressible flow are solved [156]. This approach recognized that models of turbulent RT flows should incorporate spectral nonequilibrium and allow spectral initial conditions to be used. The equations are very difficult to solve. Later, a simpler model that solves the spherically averaged (in wavevector space) density–velocity correlation and density variance spectra transport equations was proposed and applied to water channel experiments [157]. The turbulent viscosity is $\nu t(x,t)=\u222b0\u221e2\tau ii/k3dk$. The density–velocity equation contains a buoyancy source term, a form drag and a turbulent viscous drag, another source term that represents a disordered contribution to the fluid interpenetration, a forward cascade term, and a turbulent diffusion term. The density variance equation contains a source term derived from a kinematic expression of advecting mass through a control volume, a term representing the inverse cascade, and a direct cascade term combining advective and diffusive transport.

The *K*–*L* model [57] was further refined by introducing a transport equation for the mass flux *a _{i}* [59,158]. The model in Ref. [158] was referred to as the BHR-1 (or BHR) model. The

*K*–

*L*–

*a*model in Ref. [59] differs in an important way from the BHR-1 model in that the

*L*equation in the former model does not include a buoyancy production term. Both models use an algebraic model for the (negative) density–specific volume correlation

*b*, and are immiscible. Figure 28 shows a comparison between the normalized turbulent kinetic energy from the BHR-1 model with small Atwood number gas channel data [48].

The four-equation BHR-2 model was subsequently developed, replacing the algebraic model for *b* in the BHR-1 model by a modeled transport equation [159]. This formulation emphasized the modeling of molecular mixing using a transport equation that allows *b* to decrease as miscible mixing progresses. The BHR-1 and BHR-2 models (as well as all of the other previously mentioned models) used the Boussinesq model for the Reynolds stress tensor *τ _{ij}*. The BHR-3 model replaced this algebraic closure by a modeled transport equation for

*τ*[85], recognizing that buoyancy manifested by the mean pressure gradient directly affects the Reynolds stress (as also elucidated in Ref. [160]), which appears in the mean momentum and

_{ij}*a*equations. Figure 29 shows a comparison of the volume fraction and

_{i}*b*from the model to DNS and gas channel experimental data [48] for RT flow with $At\u22480.04$.

Thus far, very little theoretical work has been performed to better understand turbulent transport mechanisms or improve turbulent transport models for RT flows using simulation data validated by experimental data. Data from DNS of a transitional RT mixing layer [67] was used to investigate the spatiotemporal structure of mean and turbulent transport and mixing [160]. The budgets of the $v\u0303z,\u2009m\u03031$, *K*, *ε*, $S=m1\u20332\u0303$, and *ε _{S}* equations were constructed using Reynolds-averaged DNS data. The relative importance of turbulent production, dissipation and destruction, and transport across the mixing layer was investigated as a function of time (equivalently Reynolds number) to provide insight into the flow dynamics unavailable from experiments at the time. The budget analysis supported the assumption for small Atwood number RT flows that the principal transport mechanisms are buoyancy production, turbulent production, turbulent dissipation, and turbulent diffusion (shear production is negligible). Distinctions between momentum and scalar transport were noted:

*K*and

*ε*both grew in time and were peaked near the layer centerplane, while

*S*and

*ε*initially grew and then decreased as mixing progressed and reduced density fluctuations. All terms in the transport equations grew or decayed, with no qualitative change in their profile, except for the pressure flux contribution to the total turbulent kinetic energy flux, which changed sign early in time. The production-to-dissipation ratios corresponding to

_{S}*K*and

*S*were large and varied strongly at small times, decreased with time, and nearly asymptoted as the flow approached self-similarity.

A four-equation model based on the turbulent kinetic energy and the incompressible turbulent kinetic energy dissipation rate, *K* and *ε*, and on the heavy-fluid mass fraction variance and its dissipation rate, *S* and $\chi =\epsilon S$ [160] was proposed for small Atwood number, intermediate Reynolds number RT flows [161]. The validity of gradient-diffusion and similarity closures was investigated a priori. A methodology was used to estimate model coefficients as a function of time in the closures by minimizing the *L*_{2} norm of the difference between the model and correlations constructed using the simulation data. It was shown that gradient-diffusion and similarity closures describe the exact, unclosed profiles well. Order-of-magnitude estimates [160] for the terms in the exact transport equations and their closure models showed that several standard closures for the turbulent production and dissipation (destruction) must be modified to include Reynolds-number scalings appropriate for RT flow at small to intermediate Re. The late-time, large Reynolds number coefficients are different from those used in shear flows. It was also shown that the Boussinesq model predictions for *τ _{ij}* agree better with the data when additional buoyancy terms are included. Thus, an unsteady RANS paradigm is needed to predict transitional flow dynamics, analogous to the small Reynolds number modifications in RANS models of wall-bounded flows in which the production-to-dissipation ratio is far from equilibrium.

## 6 Future Directions in Simulation, Modeling, and Experiments

The overview of simulations and modeling of Rayleigh–Taylor unstable flows and mixing presented here has illuminated several areas of progress and trends. The early accelerated tank experiments have been the subject of most simulation and modeling efforts. Overall, there has been somewhat less effort to simulate and model experiments that were developed later, such as the stationary tank and channel experiments. The former experiments were more systematically performed over a range of Atwood numbers, but were largely limited to measurements of the large-scale mixing layer quantities $hb,s$ and $\alpha b,s$. The latter experiments began to use more sophisticated fluid mechanics diagnostics and provided more detailed statistical data–more challenging for simulations and models to predict. As a result, most simulation and modeling efforts have focused on comparisons of large-scale mixing properties to those measured in the experiments. There have been relatively few simulations that have extracted and compared statistical data to measured data. (M)ILES have dominated DNS and LES in numerical simulation studies. A large variety of modeling approaches have been applied to experiments. Most modeling has been applied to accelerated tank and channel experiments.

Ultimately, simulations and models must be validated by experimental data. Strong coupling and collaboration between experimental and numerical groups is vital to further success in these efforts. Making experimental datasets available to other researchers could also encourage applications of simulation and modeling to broader classes of experiments. Potential areas for future work include:

Encouraging measurements of initial conditions in experiments and using them in simulations (in particular, investigating possibilities for creating reproducible methods for generating variations in the spectral content, similar to what has been done numerically);

performing simulations and modeling of specific experiments, and then varying parameters away from these configurations (in contrast to performing simulation and modeling with such variations but with no relationship to any experiment);

further expanding the use of LES with subgrid-scale models appropriate for buoyancy-driven RT flows;

developing and applying simulation and modeling to predict quantities from a wider variety of different classes of experiments;

developing models that can predict additional mean and turbulent properties of mixing layers;

performing simulation and modeling of the larger Atwood and Reynolds number gas channel experiments;

performing more intercomparisons of the predictions of different numerical approaches [i.e., DNS, LES, and (M)ILES] and methods (i.e., based on different discretizations and algorithms) applied to experiments to better understand their relative advantages and disadvantages and overall accuracy;

performing more intercomparisons of the predictions of different modeling approaches and specific models applied to experiments to better understand their relative advantages and disadvantages and overall accuracy;

experimentally and numerically better understanding the differences between the liquid and gas phase instability and mixing evolution (e.g., Schmidt and Prandtl number effects);

better understanding the differences and similarities between miscible and immiscible simulations, models, and experiments;

encouraging parameterization of experimental and numerical simulation data to facilitate cross-comparisons of quantities;

applying uncertainty quantification techniques to simulation and modeling studies;

performing simulations and modeling of combined RT/KH and RT/RM instability-induced mixing experiments, and

designing an apparatus to conduct a new generation of controled acceleration/deceleration (or complex acceleration) experiments with modern diagnostics to provide new data for simulation and modeling.

A series of similarly designed gas channel experiments over a wide range of Atwood numbers (rather than separate experiments) would be beneficial for reducing the uncertainty in the measurements of *α _{b}* and

*α*and perhaps for further differentiating between such miscible mixing cases from the immiscible cases that in recent years give comparable values of

_{s}*α*and

_{b}*α*. It should be noted that, in general, the areas of future work noted above apply equally well to simulation and modeling of Richtmyer–Meshkov instability-induced mixing.

_{s}As discussed in the review of experimental work on RT flows [9], there are many opportunities for future experimental research that can complement ongoing computational and theoretical studies. Further addressing the challenges outlined in this review will open and inspire new opportunities in simulation and modeling of Rayleigh–Taylor instability, turbulence, and mixing.

## Acknowledgment

This article is dedicated to the memory of Dr. Malcolm John Andrews—friend, collaborator, colleague, and scientist extraordinaire. As apparent from this review, he was a major contributor to the scientific literature on Rayleigh–Taylor experiments, simulations, and modeling, and inspired many students and scientists in this field of fluid dynamics. He understood the importance of and promoted the synergy among experiments, simulations, and modeling. The author thanks Dr. Nicholas J. Mueschke for many years of fruitful collaboration. I am also grateful to the reviewers for suggestions that have improved this article.

This article has been authored by Lawrence Livermore National Security, LLC under Contract No. DEAC52-07NA27344 with the U.S. Department of Energy.

This document was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes.