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 [13] 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 ρ·p<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.

Rayleigh–Taylor unstable flows are inhomogeneous and anisotropic. In contrast to single-fluid flows that are typical of turbulent flows that have been most extensively studied, RT flows involve multiple fluids with different molecular transport coefficients, molecular weights etc. In the self-similar state (corresponding to large Reynolds numbers), it is generally well-established that the average displacement of the bubble and spike fronts grow self-similarly with widths hb,s(t)=αb,sAtgt2 [58], where αb,s are the dimensionless mixing layer growth parameters for bubbles and spikes, and the Atwood number At=(ρ1ρ2)/(ρ1+ρ2)>0 characterizes the degree of initial stratification of the layer (the subscripts 1 and 2 correspond to spikes and bubbles, respectively). Thus, the late-time total mixing layer width for constant acceleration is [generalized here to include a time offset t0 and a virtual origin h(0)=αAtgt02]  
h(t)=hb(t)+hs(t)=(αb+αs)Atg(t+t0)2=αAtg(t+t0)2
(1)

(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 α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 αs (and therefore, of α) 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 cs 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 [1214]. RT instability and mixing in inertial confinement fusion [15,16] capsules decreases the efficiency of thermonuclear burning [1719] 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 [2123]. 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.24.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.25.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[0,Lz]. Also, let the subscript r label the fluid (or species). Numerical simulations of RT flows have investigated (not an exhaustive list):

  1. mixing layer widths (hb,s or h) and growth parameters (αb,s or α) based on cutoffs in some mean field on the bubble and spike front sides;

  2. integral mixing layer widths (the overbars indicate averages over the directions perpendicular to the acceleration, and nW is a numerical factor) [24]  
    W(t)=nW0Lzϕ¯1ϕ¯2dz
    (2)
    for scalars ϕr={mr,fr,Xr,cr} (mass fraction, volume fraction, mole fraction, or concentration, which are dimensionless and satisfy rϕr=1 and rϕr=0);
  3. mean mass fraction, volume fraction, and concentration profiles m¯r(z,t),f¯r(z,t), and c¯(z,t);

  4. local and global molecular mixing profiles and parameters (ϕ12¯ is the scalar variance in fluid 1) [25,26]  
    Θ(z,t)=ϕ1ϕ2¯ϕ¯1ϕ¯2=1ϕ12¯ϕ¯1ϕ¯2
    (3)
     
    Θ̂(t)=0Lzϕ1ϕ2¯dz0Lzϕ¯1ϕ¯2dz=10Lzϕ12¯dz0Lzϕ¯1ϕ¯2dz
    (4)
    respectively, where Θ̂ and 1Θ̂ are also referred to as “atomic” and “chunk” mixing fractions, respectively;
  5. ratios of kinetic energy produced to released potential energy and of dissipation to released potential energy K/ΔΦ and D/ΔΦ, 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),  
    ΔΦ=Φinitial-Φfinal=gLz[(0Lz/2ρ2zdz+Lz/2Lzρ1zdz)-0Lzρ¯zdz]
    (5)
    with g >0 the constant acceleration, the kinetic energy produced is  
    K=1Lz0Lzρv2¯2dz
    (6)
    and the dissipation is the difference between the released potential energy and kinetic energy produced, D=ΔΦK;
  6. directional components of the kinetic energy spectrum and the scalar variance spectrum, and their spatial integrals [vi(k,z,t) and ϕ(k,z,t) are Fourier-transformed velocity and scalar field fluctuations, and k=kx2+ky2 is the two-dimensional wavenumber]  
    Ei(k,t)=1h(t)-hs(t)hb(t)|vi(k,z,t)|22dz
    (7)
     
    Eϕ(k,t)=1h(t)-hs(t)hb(t)|ϕ(k,z,t)|22dz
    (8)
     
    Ei(t)=0Ei(k,t)dk,Eϕ(t)=0Eϕ(k,t)dk
    (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 Θ(z,t)=0 (corresponding to ϕ1ϕ2¯=0), the two fluids are completely segregated, and where Θ(z,t)=1 (corresponding to ϕ1ϕ2¯=ϕ¯1ϕ¯2), 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 [2830] (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 αs. 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.

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 [4046] or a gas channel [4751], 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=vct (where vc 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 αs with Atwood number.

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 α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 α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:

  1. the self-similar algebraic expression that follows directly from Eq. (1) (assuming t0=0)  
    αb,sss(t)=hb,s(t)Atgt2
    (10)
    which was used in several simulation and modeling studies [25,38,5659];
  2. the Youngs [8,25,27,31,60,61] expression  
    αb,sY(t)=dhb,s/dt2Atgt=dhb,sdX
    (11)
    (where the scaled acceleration distance is X=At[0tg(t)dt]2=Atgt2 for a constant acceleration), and is an estimate from the late-time slope of hb,s plotted versus X, which allows for an offset;
  3. the Ristorcelli–Clark [62] expression  
    αb,sRC(t)=(dhb,s/dt)24Atghb,s(t)
    (12)

    and;

  4. the Livescu et al. [63] expression  
    αb,sL(t)=[hb,s(t)hb,s(t0)gAt(tt0)]2
    (13)
    which avoids computing time derivatives and allows for an offset in time t0.

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 αs/α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.

The range of scales grows with time, so the achievable Reynolds numbers are limited by experimental facility or computational domain size. The scaling of the mixing layer (or outer-scale) Reynolds number (ν is the mixture kinematic viscosity and neglecting t0)  
Reh(t)=h(t)vh(t)ν=2(αAtg)2t3ν
(14)
[vh(t)=dh/dt and the last expression corresponds to a self-similar state] with α, At, g, and ν provides the values of Reh attainable for various values of these parameters at a given time t. With the mixture Schmidt number Sc=ν/D (D is the mixture mass diffusivity) the scalar analog of Eq. (14) is the mixing layer Péclet number Peh(t)=ScReh(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, (t)/d(t)h(t)/d(t)Reh(t)3/4, with Eq. (14) shows that the Kolmogorov scale d(t)[ν3/(8t)]1/4/αAtg decreases very slowly with time [62,64,65] (also see Ref. [54]). The analogous diffusive lengthscale is the Corrsin–Obukhov scale dϕ(t)=d(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 d and dϕ. 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 Sc0.7 [48], making them more attractive for DNS.

Another grid resolution criterion follows from the requirement to suppress the growth of spurious grid-generated modes. Viscous dissipation and mass diffusion stabilize such growth if the grid Grashof number [56] satisfies (Δρ=ρ1-ρ2)  
GrΔ=gΔρΔxΔyΔzρ(ν+D)2<1
(15)

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

 
Δz<[D2(1+Sc)2ρgΔρ]1/3

3.2.3 Atwood Number Dependence.

Another limitation arises from the increasing asymmetry between the bubble and spike mixing layer fronts with increasing Atwood number. For turbulent flows the bubble and spike fronts both grow quadratically in time according to Eq. (1), but while αb is a weak function of At, αs strongly depends on At [33,58,61,68,69] (also see Ref. [9]). Immiscible mixing experiments [33] indicated that  
αs(At)αbc1(ρ1ρ2)c2c1(1+At1At)c2
(16)

with c20.33±0.05 and c11; 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 c20.185 and c10.958 [58] and c20.18 and c11.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 αs 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 At1, 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 At0.5, although simulations at larger values have been increasingly performed in the last decade.

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

Multimode interfacial (density) perturbations are typically expressed by a linear superposition of modes (assuming that the z-direction is parallel to the acceleration, and periodic boundary conditions are taken in the x- and y-directions) [26,72,73]:  
zint(x,y)=m=mminmmaxn=nminnmaxamncos(kmx)cos(kny)+bmncos(kmx)sin(kny)+cmnsin(kmx)cos(kny)+dmnsin(kmx)sin(kny)
(17)
where the wavenumbers are {km,kn}={2πm/Lx,2πn/Ly} (m and n are integers) with wavelengths {λx,λ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=kminkmaxP(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)]  
ϕ1(x,0)={zint(x,y)/Δzifzint>01+zint(x,y)/Δzifzint<0
(18)

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

Initial velocity perturbations can be obtained using linear theory (i.e., exponential growth of modes eωkt with growth rate ωk) and the formal time derivative of Eq. (17) [26]  
vz(x,0)=m=mminmmaxn=nminnmaxωk[amncos(kmx)cos(kny)+bmncos(kmx)sin(kny)+cmnsin(kmx)cos(kny)+dmnsin(kmx)sin(kny)]
(19)
which has the same properties as Eq. (17), and vx(x,0)=vy(x,0)=0. Alternatively, initial velocity perturbations can be computed from the gradient of a potential such as [67,74]  
φ(x)=m=mminmmaxn=nminnmaxamnkmnsin(kmx)sin(kny)ekmn|z|
(20)
where the amplitude coefficients {amn} are sampled from a random distribution and kmn=km2+kn2, i.e., [67]  
v(x,0)=φ(x)Dlnρ(x,0)
(21)

with the initial density ρ(x,0) [see Eq. (25)]. The potential (20) satisfies periodic boundary conditions in the x- and y-directions, and could be generalized to a four-term expression similar to Eq. (17).

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.

Following a summary of the governing equations, a brief overview of implicit large-eddy, direct, and large-eddy simulation approaches (see Refs. [66] and [87]) is given below.

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.
The equations governing the mixing of fluids are the incompressible variable-density continuity, Navier–Stokes, and mass fraction equations (assuming Fickian diffusion) [64]  
ρt+xj(ρvj)=0
(22)
 
t(ρvi)+xj(ρvivj+pδij)=ρgi+σijxj
(23)
 
t(ρmr)+xj(ρmrvj)=Jr,jxj
(24)
respectively, where the subscripts i, j, and k are vector indices corresponding to Cartesian coordinates, σij=μ[vi/xj+vj/xi(2/3)δijvk/xk] is the viscous stress tensor (ignoring bulk viscosity), D is the binary diffusion coefficient, and Jr,j=ρDmr/xj 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ρ/ρr. It can be shown that the velocity divergence is [89,90]  
vjxj=xj(Dρρxj)0
(25)

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

4.1.1.2 Compressible approximation.
The compressible equations expressing mass, momentum, internal energy, and species conservation solved numerically for ideal gases are [88,91,92] Eqs. (22)(24) together with  
t(ρU)+xj(ρUvj)=σijvixjpvjxjxj(Φj+r=12hrJr,j)
(26)
where the internal energy is U=r=12mrUr, the kinetic energy dissipation rate is σijvi/xj, the radiative flux is Φj=κT/xj, the radiative conductivity is κ=ρcpχ, the thermal conductivity is χ, the enthalpy of fluid r is hr=Ur+pr/ρr, the mixture molecular transport coefficients (ϕ=μ,χ,D and m0.7) are  
ϕ=ϕ1m1/M1+ϕ2m2/M2m1/M1+m2/M2,ϕr=ϕr0(TT0)m
(27)

the mixture ratio of specific heats is γ=cp/cv (cp=cp1m1+cp2m2 and cv=cpRu/M=cv1m1+cv2m2) and the mixture pressure is p=ρRuT/M=(γ1)ρU with mixture molecular weight M=(m1/M1+m2/M2)1 and universal gas constant Ru. 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 [9597]. 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, Φ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] ϕ¯(x,t)=GΔ(x,x)ϕ(x,t)d3x, where GΔ(x,x) is the filter kernel satisfying a number of consistency requirements and Δ is the filter width (related to the grid spacing). For a compressible flow, ϕ̃(x,t)=[1/ρ¯(x,t)]GΔ(x,x)ρ(x,t)ϕ(x,t)d3x. 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 ρ1/ρ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/SF6 [28]. The grid resolutions for the At=1/5 case were 803 and 1603, and for the At=1/2 and 19/21 cases were 120×802 and 240×1602. 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 αbY0.04 and αbss0.054 from the finest grid simulations, which underestimated the experimental value αbY0.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×1602 simulation with At=0.5 gave an even smaller value αbY0.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 αbY values persisted for much higher resolution 720×6002 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 [106108]. 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×5762 [106]. The simulation gave αY=0.055 and the experimental value was α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=ρ(1/ρ)¯ (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×256 and 2562×512) multimode RT instability evolution with At0.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αbYAtgt2 at late times, with the codes solving the dissipative equations predicting αbY=0.025±0.003; the small value of α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 αbY that considerably underestimated the experimental values αb0.057±0.008.

Numerical simulations (2562×512 and 5123) 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 50% and the molecular mixing increases 100%. 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 (903 and 1803) 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 ρ1/ρ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×1045×102). 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 αbY=0.044±0.005 and Θ(z,t)0.75 near the center of the layer. The simulation gave αbY0.035 and Θ̂(t)0.81 at resolution 1682×230 in good agreement with the experimental values.

The simulations gave mixing efficiencies η=1ΔΦ/Φmax=0.48 (random initial perturbations) and 0.47 (with long-wavelength initial perturbation), where ΔΦ is the final potential energy loss and Φ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 pH 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 At2×103, 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 ρ1/ρ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 1 while the molecular Schmidt number was Sc103; 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×80×200: “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 k5/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)t2/5 similarity solution to the diffusion equation in an infinite domain was obtained (differing from the classical t2 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×2560 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)k3. 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×256×128 MILES using the RTI-3D code and initialized using experimental water channel velocity and density data (At=7.5×104) [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 amn, bmn, cmn, and dmn.

The azimuthally averaged kinetic energy spectra used in the simulations had the same energy for the included wavenumbers as in the experiment: kminkmaxEsim(k)θdk=kminkmaxEexp(k)dk. The density surface perturbations were converted to volume fraction fluctuations using Eq. (18) (with ϕf), and initial velocity perturbations were obtained from the gradient of the potential (20). PIV data were used to obtain the velocity data, and the amn 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 αbcl0.06 [based on the centerline value of v/(2Atgt), where v 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×760×1280 DNS using initial conditions, geometry, and physical parameters approximating those of a transitional water channel RT mixing experiment with At=7.5×104 [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 hb and the mixing layer growth parameter αb (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 αbY,αbRC0.07 at the latest time, in excellent agreement with the experiment. The late-time value Θ̂(t)0.55 also compared favorably with the experimental value 0.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)1710. 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].

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 [106108]. Simulations (288×180×240) 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 α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, αY=0.069 and αexpY=0.065–0.07; (2) the Sc=7 water channel experiment, αY=0.075 and αexpY=0.077±0.011, and; (3) the Sc=620 water channel experiment, αY=0.084 and αexpY=0.085±0.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.

A modal growth model [119] was developed to better understand the weakly nonlinear regime of multimode RT instability. The model was derived from application of the Bernoulli equation and potential theory truncated to second-order, and was originally used to estimate RT instability growth at unstable interfaces in inertial confinement fusion capsules [120]. Substituting the Fourier expansions of the velocity potentials φr(x,t)=kφkr(t)eik·x-kz(δr1-δr2), where k (and later p) is the two-dimensional wavevector into the interface evolution and Bernoulli equations, and using a Taylor expansion of φr/z truncated at first-order gives the Fourier amplitude evolution equation zk/t=kφk1(t)pk·pφp1(t)zkp(t)+=kφk2(t)pk·pφp2(t)zkp(t)+ Substituting φk1,2=±(1/k)zk/t+p[k·p/(kp)](zp/t)zkp+ and their time derivatives into the Fourier-transformed Bernoulli equation gives the second-order amplitude evolution equation (with k̂=k/|k|,ξ=1/2k̂·p̂(k̂p̂)·p̂/2, and ζ=1k̂·p̂)  
(d2dt2ωk2)zk(t)=Atkpζd2zpdt2zkp+ξdzpdtdzkpdt+,
(28)

where ωk=gAtk is the RT instability growth rate.

The nonlinear coupled mode equations may be solved iteratively beginning with the linear solution zklin(t)=zk(0)cosh(ωkt) of Eq. (28) neglecting the second-order terms and gives  
(d2dt2ωk2)zk(t)=Atkpzp(0)zkp(0)×{F+cosh[(ωp+ω|kp|)t]+Fcosh[(ωpω|kp|)t]}
(29)
 
F±ωp2(ωpζ±ω|kp|ξ)
(30)
The solution of Eq. (29) is  
zk(t)=zklin(t)+AtkpWk,p(t)zp(0)zkp(0)
(31)

with Wk,p(t)=G+{cosh[(ωp+ω|kp|)t]cosh(ωkt)}+G{cosh[(ωpω|kp|)t]cosh(ωkt)} and G±F±/[(ωp±ω|kp|)2ωk2]. The perturbation expansion converges provided that the mode amplitudes satisfy zk(t)<0.1λ, where λ=2π/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.

The general form of the buoyancy–drag equations is ρbdvb/dt=nFbn and ρsdvs/dt=nFsn, where the bubble/spike front velocities are vb,s=dhb,s/dt and the forces {Fb,sn} consist of a buoyancy, drag, added (or virtual) mass, and pressure gradient force (other forces, such as a shear force [123], can also be included). A single equation for 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 [122126]  
ρ3rdv3rdt=Fb3r(t)+Fvm3r(t)+Fd3r(t)
(32)
with the general buoyancy force, virtual mass (plus pressure) force, and drag force (Cvm=Cvm+1)  
Fb3r(t)=Cb[1E(t)](ρ2ρ1)g(t)
(33)
 
Fvm3r(t)={CvmE(t)ρ3r+[Cvm+E(t)]ρr}dv3rdt
(34)
 
Fd3r(t)=Cdbρr|v3r(t)|v3r(t)3r(t)×[Rrnδr1+Rr(h3rhr)βδr2]
(35)
with exponential factor E(t)=e-Cekhb(t) (with Ce=2 and 3 for three and two dimensions, respectively), lengthscale (i.e., radius) [127] for multimode and single-mode mixing 3r(t)=V3r(t)/[bS3r(t)]=h3r(t) and 3r(t)λ, respectively, and geometrical factor b =2/3 and 1 in three and two dimensions, respectively. The velocity is v3r(t)=dh3r/dt. The factor Rr=ρr/(ρr+ρ3r) introduces a density dependence in the drag terms, and [hs(t)/hb(t)]β 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 ϕ=ϕ¯+ϕ and ϕ=ϕ̃+ϕ for incompressible and compressible fields, respectively: ϕ¯ and ϕ̃=ρϕ¯/ρ¯ 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.

For example, Reynolds-averaging the incompressible variable-density Eqs. (22)(24) gives  
ρ¯t+xj(ρ¯ṽj)=0
(36)
 
t(ρ¯ṽi)+xj(ρ¯ṽiṽj+p¯δij)=ρ¯gi+xj(σ¯ijτij)
(37)
 
t(ρ¯m̃r)+xj(ρ¯m̃rṽj)=xj(J¯j,r+Γj,r)
(38)
where the mean mass fraction flux is J¯j,r=ρ¯D¯m̃r/xj, and the Reynolds stress tensor and the mass fraction flux vector are τij=ρ¯vivj̃ and Γj,r=ρ¯mrvj̃. These terms are closed using functionals of the mean fields. The simplest and most commonly used models are based on the gradient-diffusion hypothesis, in which the stress and flux are modeled by the algebraic closures [S̃ij=(1/2)(ṽi/xj+ṽj/xi) is the mean strain-rate tensor and S̃kk=ṽk/xk is the mean velocity divergence]  
τij=23ρ¯Kδij2μt(S̃ijδij3S̃kk),Γj,r=μtσmm̃rxj
(39)

with σm a dimensionless turbulent Schmidt number, K the turbulent kinetic energy, μt=ρ¯νt, and ν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 τij and Γ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.

For example, in the Kε model νt=CμK2/ε is computed using the modeled transport equations (neglecting the pressure–dilatation and dilatation dissipation) [132,133]  
t(ρ¯K)+xj(ρ¯Kṽj)=ajp¯xjτijṽixjρ¯ε+xj[(μ¯+μtσK)Kxj]
(40)
 
t(ρ¯ε)+xj(ρ¯εṽj)=Cε0εKajp¯xjCε1εKτijdṽixj23Cε3ρ¯εṽjxjCε2ρ¯ε2K+xj[(μ¯+μtσε)εxj]
(41)
where τijd=2μt(S̃ijδijS̃kk/3) is the deviatoric Reynolds stress tensor, Cε0Cε3 are dimensionless coefficients, and σK and σε are dimensionless turbulent Schmidt numbers. The first two terms on the right sides of Eqs. (40) and (41) represent buoyancy and shear production, respectively, and the last terms are turbulent diffusion terms. The mass flux velocity in the buoyancy production terms is modeled by  
aj=vj¯=ρvj¯ρ¯=νtσρρ¯(ρ¯xjρ¯γ¯p¯p¯xj)
(42)

where σρ 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=γp¯/ρ¯ 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¯ and f¯, of classical incompressible RT mixing was derived [134]. Asymptotic analysis was used to derive density-ratio-invariant hybrid profiles C(X)=Fm¯+(1F)f¯ with At=(R1)/(R+1) and F=(1+At)(1X)/2 and obtain an analytical expression for C(X), where X[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¯=m¯/[m¯+(1m¯)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.

A model based on the exponential growth of a small amplitude wavepacket of modes was used to obtain expressions for the bubble mixing parameter and bubble aspect ratio [72]  
αb(h0k)=Cπ4[ln(2Cπkbh0k)1]1
(43)
 
βb(h0k)=λbhb=2πC[ln(2Cπkbh0k)1]1
(44)
which depend on the initial perturbation amplitude h0k, and C is a free parameter. Applying self-similarity shows that αb increases proportional to the bubble Froude number (Db is the bubble diameter)  
Fr=vbρhρbρhgDb2
(45)
and logarithmically with h0k. The model has two parameters that can be determined by experimental or numerical simulation data. An alternate form of Eqs. (43) and (44) is  
αb(η0)=πFr8(1+At)[ln(2Frπ(1+At)η0)1]1
(46)
 
βb(η0)=2π(1+At)Fr[ln(2Frπ(1+At)η0)1]1
(47)

with η0=h0k/λ0=kh0k/(2π); 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 βb with scaled rms initial amplitude for Fr=0.49 and 0.82 from the models and numerical simulations.

A mode-coupling model for RT interface evolution that correctly models interface turnover as found in rocket-rig experiments was developed more recently [136] (see Ref. [79]). The model consists of equations for the modal amplitude in the interface-parallel direction ξ, modal amplitude in the interface-normal direction η, and vorticity amplitude ω (integrated across the interface)  
ξt=H(ω)2kη||(1+kξ,kη)||εak2ξ
(48)
 
ηt=H(ω)21+kξ||(1+kξ,kη)||εak2η
(49)
 
ωt=At2H(ωH(ω))||(1+kξ,kη)||2Atgηεak2ω
(50)

where k is the wavenumber, H(·) is the Hilbert transform, and ε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ξ,kη)|| is the slope of the interface and (a,b)ae1+be2.

A mean field formulation of the Haan–Ofer–Shvarts mode coupling model previously developed for RT instability [119,135]  
(d2dt2ωk2)zk(t)=AtkpSk,p
(51)
 
Sk,p=ωp(ωpζ+ω|kp|ξ)zp(t)zkp(t)
(52)

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 pSk,pSk=pSk,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λ/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)=2kzk2. 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 α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.

A buoyancy–drag model was derived from the Euler–Lagrange equations [138]. The Lagrangian is based on the kinetic energy of a linearly perturbed interface derived from a velocity potential, φ, and the resulting equation for 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  
T=1LxLyρ|φ|22d3x
(53)
The Euler–Lagrange equation corresponding to Eq. (53) evaluated to second-order in h and v is (d/dt)(T/v)T/h=Q, where Q=c(ρ1+ρ2)v|v|/2 is a dissipative (drag) force. According to the wavelength renormalization hypothesis, the single-mode perturbation wavelength λ is replaced by λ(t)=b|h(t)|, with bO(1) [27,139142]; n.b., b is not the geometrical factor introduced previously in Sec. 5.1.2. For RT instability  
hdvdt+(12+2πcb)v2=2πg(t)Atbh
(54)

with solution Eq. (1) and α=π/[2(b+2πc)]. The transition between the linear and nonlinear regimes may be achieved with λ(t)=max[λ0,b|h|+(1mb)λ0], where m is a parameter. Figure 16 shows a comparison of model predictions with experimental data for four different acceleration histories.

A buoyancy–drag model was developed for RT mixing using the premise that the bubble and spike regions behave as distinct, spanwise homogeneous fluids [126]. Mass conservation was applied across the mixing layer, together with a piecewise-linear approximation for the density, to obtain the average mixture densities dynamically, which were then used to calculate the inertia and buoyancy terms. The only model parameter is the drag coefficient C2.5±0.6, which was determined from turbulent RT mixing experiments over various Atwood numbers and acceleration histories. The bubble and spike widths obey the hr=αrAtgt2 scaling for constant 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 ρ0=(ρ1h1+ρ2h2)/(h1+h2) is obtained from mass conservation (ρ2ρ0)h2/2=(ρ0ρ1)h1/2. The average mixture densities ρM,r=(ρr+ρ0)/2 can be expressed as  
ρM,3r=(ρ1+ρ2)hr+2ρ3rh3r2(h1+h2)
(55)
Such a representation satisfies the continuity equation for incompressible flows and allows analytic solutions to the evolution equation; additional regions are required to represent variable profiles and compressible flows. The equations for the bubble and spike sides can be constructed by considering the regions as distinct buoyant parcels and can be combined into  
dv3rdt=Cb,3r(R,χ)gAtCd,3r(R,χ)v3r|v3r|h3r
(56)
with normalized buoyancy and drag coefficients  
Cb,1(R,χ)=1+R1+R+2χ,Cb,2(R,χ)=(1+R)χ(1+R)χ+2R,Cd,1(R,χ)=2C(1+χ)1+R+2χ,Cd,2(R,χ)=2C(1+χ)R(1+R)χ+2R
(57)

where R=ρ1/ρ2 and χ=h1/h2. Figure 17 shows αs/αb (i.e., α1/α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 KL 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 SF6/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 KL 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 αbss and α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×480 and 5122×768 ILES using the RTI-3D code. The quantities compared were visualizations of 4f1f2, contours of f1, contours of K, az, 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 KL, KLai, KLaib, 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 KLai model is more accurate than the KL 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 KLai model, particularly in the prediction of ai and K in the tilted-rig case; (4) with respect to the tilted-rig case, the BHR-2 model performs better than the KLaib 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.

A simplified form of a previously developed two-fluid, two-velocity model [7,31] motivated by the fact that the mixing layer grows due to a pressure gradient in the mixed region, which induces different accelerations in the two fluids, was developed and evaluated for RT flows [37]. The velocity difference of the fluids is suppressed by drag. The continuity and Navier–Stokes equations are  
t(frρr)+xj(frρrvjr)=0
(58)
 
t(frρrvir)+xj(frρrvirvjr)=fr(ρrgipxi)+Kd(vi3rvir)
(59)

where fr, ρr, and vir are the volume fraction, density, and velocity of fluid r, gi is the constant gravitational acceleration, Kd is the interphase drag, and f1+f2=1. Models with various algebraic and differential lengthscales (mixing lengths) entering the expression for Kd 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.

5.3.2 Turbulent Transport Modeling of Tank Experiments With Energy Deposition.

RT mixing with volumetric energy deposition was investigated experimentally using a benchtop microwave facility in which a light, weakly polar fluid initially rested on top of a heavier, higher polarity fluid [149,150]. During the application of microwave radiation, less energy is deposited into the weakly polar fluid than into the higher polarity fluid, resulting in preferential heating and density decrease of the bottom fluid, inducing a buoyancy instability. The fluids used were toluene and tetrahydrofuran. It was demonstrated that volumetric energy deposition increases the rate of mixing compared to classical RT mixing. A zero-dimensional (ordinary differential equation) Kε model in which the gradients of turbulent kinetic energy and its dissipation rate are zero was proposed:  
dKdt=Cμσtgρ0ρρ0hK2εε,dεdt=Cε0Cμσtgρ0ρρ0hKCε2ε2K
(60)

where h is the mixing layer width. Experimental data indicated that ht3 at late time, and self-similarity gives K(dh/dt)2t4 and εK/tt3 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×1024 DNS using the TRICLADE code at At=0.1 using an initial annular spectrum. Figure 25 shows the relationship between α and Θ̂ 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 αs values. The model leads to a simple correction to the formulation of the mixing layer growth, which is consistent with α 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 αb and αs as a function of At, together with the theoretical curves predicted by the model.

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.

One- and two-dimensional incompressible Boussinesq Kε models were applied to small Atwood number RT mixing [153155]. 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  
αb=CμσT(Cε0Cε2)2(4Cε03)(4Cε23)
(61)

Comparison of the model predictions with the experimental data was used to estimate Cε00.9 and σT=0.6 (with the other coefficients having their standard values); an equivalent formulation using the mean density gives the same result with σTσρ [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 νt(x,t)=02τ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 KL model [57] was further refined by introducing a transport equation for the mass flux ai [59,158]. The model in Ref. [158] was referred to as the BHR-1 (or BHR) model. The KLa 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 τij [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 ai equations. Figure 29 shows a comparison of the volume fraction and b from the model to DNS and gas channel experimental data [48] for RT flow with At0.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 ṽz,m̃1, K, ε, S=m12̃, 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 εS 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 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 χ=ε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 L2 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 α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:

  1. 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);

  2. 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);

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

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

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

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

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

  8. 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;

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

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

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

  12. applying uncertainty quantification techniques to simulation and modeling studies;

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

  14. 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 αs and perhaps for further differentiating between such miscible mixing cases from the immiscible cases that in recent years give comparable values of αb and αs. 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.

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.

References

References
1.
Rayleigh
,
Lord
,
1882
, “
Investigation of the Character of the Equilibrium of an Incompressible Heavy Fluid of Variable Density
,”
Proc. London Math. Soc.
,
14
(
1
), pp.
170
177
.10.1112/plms/s1-s14.1.170
2.
Taylor
,
G.
,
1950
, “
The Instability of Liquid Surfaces When Accelerated in a Direction Perpendicular to Their Planes. I
,”
Proc. R. Soc. London A
,
201
(
1065
), pp.
192
196
.10.1098/rspa.1950.0052
3.
Chandrasekhar
,
S.
,
1961
,
Hydrodynamic and Hydromagnetic Stability
,
Dover
,
New York
.
4.
Sharp
,
D. H.
,
1984
, “
An Overview of Rayleigh-Taylor Instability
,”
Phys. D
,
12
(
1–3
), pp.
3
18
.10.1016/0167-2789(84)90510-4
5.
Belen'kii
,
S. Z.
, and
Fradkin
,
E. S.
,
1965
, “
Theory of Turbulent Intermixing
,”
Proceedings of the P. N. Lebedev Physics Institute
Vol.
29
,
Quantum Field Theory and Hydrodynamics. P. N. Lebedev Institute
,
Nauka, Moscow, Russia
, pp.
197
228
.
6.
Anuchina
,
N. N.
,
Kucherenko
,
Y. A.
,
Neuvazhaev
,
V. E.
,
Ogibina
,
V. N.
,
Shibarshov
,
L. I.
, and
Yakovlev
,
V. G.
,
1979
, “
Turbulent Mixing at an Accelerating Interface Between Liquids of Different Density
,”
Fluid Dyn.
,
13
(
6
), pp.
916
920
.10.1007/BF01050969
7.
Youngs
,
D. L.
,
1984
, “
Numerical Simulation of Turbulent Mixing by Rayleigh–Taylor Instability
,”
Phys. D
,
12
(
1–3
), pp.
32
44
.10.1016/0167-2789(84)90512-8
8.
Youngs
,
D. L.
,
1992
, “
Rayleigh-Taylor Instability: Numerical Simulation and Experiment
,”
Plasma Phys. Controlled Fusion
,
34
(
13
), pp.
2071
2076
.10.1088/0741-3335/34/13/042
9.
Banerjee
,
A.
,
2020
, “
Rayleigh Taylor Instability: A Status Review of Experimental Designs and Measurement Diagnostics
,”
ASME J. Fluids Eng.
,
142
(
12
), p.
120801
.10.1115/1.4048349
10.
Zhou
,
Y.
,
2017
, “
Rayleigh–Taylor and Richtmyer–Meshkov Instability Induced Flow, Turbulence, and Mixing. I
,”
Phys. Rep.
,
720–722
, pp.
1
136
.10.1016/j.physrep.2017.07.005
11.
Zhou
,
Y.
,
2017
, “
Rayleigh–Taylor and Richtmyer–Meshkov Instability Induced Flow, Turbulence, and Mixing. II
,”
Phys. Rep.
,
723–725
, pp.
1
160
.10.1016/j.physrep.2017.07.008
12.
Khokhlov
,
A. M.
,
Oran
,
E. S.
, and
Wheeler
,
J. C.
,
1997
, “
Deflagration-to-Detonation Transition in Thermonuclear Supernovae
,”
Astrophys. J.
,
478
(
2
), pp.
678
688
.10.1086/303815
13.
Swisher
,
N. C.
,
Kuranz
,
C. C.
,
Arnett
,
D.
,
Hurricane
,
O.
,
Remington
,
B. A.
,
Robey
,
H. F.
, and
Abarzhi
,
S. I.
,
2015
, “
Rayleigh-Taylor Mixing in Supernova Experiments
,”
Phys. Plasmas
,
22
(
10
), p.
102707
.10.1063/1.4931927
14.
Nouri
,
A. G.
,
Givi
,
P.
, and
Livescu
,
D.
,
2019
, “
Modeling and Simulation of Turbulent Nuclear Flames in Type Ia Supernovae
,”
Prog. Aerosp. Sci.
,
108
, pp.
156
179
.10.1016/j.paerosci.2019.04.004
15.
Atzeni
,
S.
, and
Meyer-ter-Vehn
,
J.
,
2004
,
The Physics of Inertial Fusion: Beam Plasma Interaction, Hydrodynamics, Hot Dense Matter
(International Series of Monographs on Physics), Vol.
125
,
Oxford University Press
,
Oxford, UK
.
16.
Pfalzner
,
S.
,
2006
,
An Introduction to Inertial Confinement Fusion
(Series in Plasma Physics),
CRC Press
,
New York
.
17.
Freeman
,
J. R.
,
Clauser
,
M. J.
, and
Thompson
,
S. L.
,
1977
, “
Rayleigh-Taylor Instabilities in Inertial-Confinement Fusion Targets
,”
Nucl. Fusion
,
17
(
2
), pp.
223
230
.10.1088/0029-5515/17/2/005
18.
Gupta
,
N. K.
, and
Lawande
,
S. V.
,
1989
, “
Rayleigh–Taylor Instability in Multi-Structured Inertial Confinement Fusion Targets
,”
Laser Particle Beams
,
7
(
1
), pp.
27
54
.10.1017/S0263034600005826
19.
Lindl
,
J. D.
,
1998
,
Inertial Confinement Fusion: The Quest for Ignition and Energy Gain Using Indirect Drive
,
AIP Press
,
New York
.
20.
Smalyuk
,
V. A.
,
Robey
,
H. F.
,
Casey
,
D. T.
,
Clark
,
D. S.
,
Döppner
,
T.
,
Haan
,
S. W.
,
Hammel
,
B. A.
,
MacPhee
,
A. G.
,
Martinez
,
D.
,
Milovich
,
J. L.
,
Peterson
,
J. L.
,
Pickworth
,
L.
,
Pino
,
J. E.
,
Raman
,
K.
,
Tipton
,
R.
,
Weber
,
C. R.
,
Baker
,
K. L.
,
Bachmann
,
B.
,
Hopkins
,
L. F. B.
,
Bond
,
E.
,
Caggiano
,
J. A.
,
Callahan
,
D. A.
,
Celliers
,
P. M.
,
Cerjan
,
C.
,
Dixit
,
S. N.
,
Edwards
,
M. J.
,
Felker
,
S.
,
Field
,
J. E.
,
Fittinghoff
,
D. H.
,
Gharibyan
,
N.
,
Grim
,
G. P.
,
Hamza
,
A. V.
,
Hatarik
,
R.
,
Hohenberger
,
M.
,
Hsing
,
W. W.
,
Hurricane
,
O. A.
,
Jancaitis
,
K. S.
,
Jones
,
O. S.
,
Khan
,
S.
,
Kroll
,
J. J.
,
Lafortune
,
K. N.
,
Landen
,
O. L.
,
Ma
,
T.
,
MacGowan
,
B. J.
,
Masse
,
L.
,
Moore
,
A. S.
,
Nagel
,
S. R.
,
Nikroo
,
A.
,
Pak
,
A.
,
Patel
,
P. K.
,
Remington
,
B. A.
,
Sayre
,
D. B.
,
Spears
,
B. K.
,
Stadermann
,
M.
,
Tommasini
,
R.
,
Widmayer
,
C. C.
,
Yeamans
,
C. B.
,
Crippen
,
J.
,
Farrell
,
M.
,
Giraldez
,
E.
,
Rice
,
N.
,
Wilde
,
C. H.
,
Volegov
,
P. L.
, and
Johnson
,
M. G.
,
2017
, “
Mix and Hydrodynamic Instabilities on NIF
,”
J. Instrum.
,
12
(
6
), p.
C06001
.10.1088/1748-0221/12/06/C06001
21.
Houseman
,
G. A.
, and
Molnar
,
P.
,
1997
, “
Gravitational (Rayleigh–Taylor) Instability of a Layer With Non-Linear Viscosity and Convective Thinning of Continental Lithosphere
,”
Geophys. J. Int.
,
128
(
1
), pp.
125
150
.10.1111/j.1365-246X.1997.tb04075.x
22.
Gorczyk
,
W.
,
Hobbs
,
B.
, and
Gerya
,
T.
,
2012
, “
Initiation of Rayleigh–Taylor Instabilities in Intra-Cratonic Settings
,”
Tectonophysics
,
514–517
, pp.
146
155
.10.1016/j.tecto.2011.10.016
23.
Molnar
,
P.
, and
Houseman
,
G. A.
,
2013
, “
Rayleigh-Taylor Instability, Lithospheric Dynamics, Surface Topography at Convergent Mountain Melts, and Gravity Anomalies
,”
J. Geophys. Research-Solid Earth
,
118
(
5
), pp.
2544
2557
.10.1002/jgrb.50203
24.
Andrews
,
M. J.
, and
Spalding
,
D. B.
,
1990
, “
A Simple Experiment to Investigate Two-Dimensional Mixing by Rayleigh-Taylor Instability
,”
Phys. Fluids A
,
2
(
6
), pp.
922
927
.10.1063/1.857652
25.
Youngs
,
D. L.
,
1991
, “
Three-Dimensional Numerical Simulation of Turbulent Mixing by Rayleigh-Taylor Instability
,”
Phys. Fluids A
,
3
(
5
), pp.
1312
1320
.10.1063/1.858059
26.
Dimonte
,
G.
,
Youngs
,
D. L.
,
Dimits
,
A.
,
Weber
,
S.
,
Marinak
,
M.
,
Wunsch
,
S.
,
Garasi
,
C.
,
Robinson
,
A.
,
Andrews
,
M. J.
,
Ramaprabhu
,
P.
,
Calder
,
A. C.
,
Fryxell
,
B.
,
Biello
,
J.
,
Dursi
,
L.
,
MacNeice
,
P.
,
Olson
,
K.
,
Ricker
,
P.
,
Rosner
,
R.
,
Timmes
,
F.
,
Tufo
,
H.
,
Young
,
Y.-N.
, and
Zingale
,
M.
,
2004
, “
A Comparative Study of the Turbulent Rayleigh-Taylor Instability Using High-Resolution Three-Dimensional Numerical Simulations: The Alpha-Group Collaboration
,”
Phys. Fluids
,
16
(
5
), pp.
1668
1693
.10.1063/1.1688328
27.
Youngs
,
D. L.
,
1994
, “
Numerical Simulation of Mixing by Rayleigh–Taylor and Richtmyer–Meshkov Instabilities
,”
Laser Particle Beams
,
12
(
4
), pp.
725
750
.10.1017/S0263034600008557
28.
Read
,
K. I.
, and
Youngs
,
D. L.
,
1983
, “
Experimental Investigation of Turbulent Mixing by Rayleigh–Taylor Instability
,” Atomic Weapons Research Establishment, Aldermaston, UK, Report No. 011/83.
29.
Read
,
K. I.
,
1984
, “
Experimental Investigation of Turbulent Mixing by Rayleigh–Taylor Instability
,”
Phys. D
,
12
(
1–3
), pp.
45
58
.10.1016/0167-2789(84)90513-X
30.
Smeeton
,
V. S.
, and
Youngs
,
D. L.
,
1987
, “
Experimental Investigation of Turbulent Mixing by Rayleigh-Taylor Instability—Part 3
,” Atomic Weapons Establishment, Aldermaston, UK, Report No. 35/87.
31.
Youngs
,
D. L.
,
1989
, “
Modelling Turbulent Mixing by Rayleigh–Taylor Instability
,”
Phys. D
,
37
(
1–3
), pp.
270
287
.10.1016/0167-2789(89)90135-8
32.
Dimonte
,
G.
, and
Schneider
,
M.
,
1996
, “
Turbulent Rayleigh-Taylor Instability Experiments With Variable Acceleration
,”
Phys. Rev. E
,
54
(
4
), pp.
3740
3743
.10.1103/PhysRevE.54.3740
33.
Dimonte
,
G.
, and
Schneider
,
M.
,
2000
, “
Density Ratio Dependence of Rayleigh-Taylor Mixing for Sustained and Impulsive Acceleration Histories
,”
Phys. Fluids
,
12
(
2
), pp.
304
321
.10.1063/1.870309
34.
Duff
,
R. E.
,
Harlow
,
F. H.
, and
Hirt
,
C. W.
,
1962
, “
Effects of Diffusion on Interface Instability Between Gases
,”
Phys. Fluids
,
5
(
4
), pp.
417
425
.10.1063/1.1706634
35.
Olson
,
D. H.
, and
Jacobs
,
J. W.
,
2009
, “
Experimental Study of Rayleigh–Taylor Instability With a Complex Initial Perturbation
,”
Phys. Fluids
,
21
(
3
), p.
034103
.10.1063/1.3085811
36.
Horne
,
J. T.
, and
Lawrie
,
A. G. W.
,
2020
, “
Aspect-Ratio-Constrained Rayleigh–Taylor Instability
,”
Phys. D
,
406
, p.
132442
.10.1016/j.physd.2020.132442
37.
Andrews
,
M. J.
,
1986
, “
Turbulent Mixing by Rayleigh-Taylor Instability
,” Ph.D. Thesis,
Imperial College of Science and Technology
, London, UK.
38.
Linden
,
P. F.
,
Redondo
,
J. M.
, and
Youngs
,
D. L.
,
1994
, “
Molecular Mixing in Rayleigh–Taylor Instability
,”
J. Fluid Mech.
,
265
, pp.
97
124
.10.1017/S0022112094000777
39.
Dalziel
,
S. B.
,
Linden
,
P. F.
, and
Youngs
,
D. L.
,
1999
, “
Self-Similarity and Internal Structure of Turbulence Induced by Rayleigh–Taylor Instability
,”
J. Fluid Mech.
,
399
, pp.
1
48
.10.1017/S002211209900614X
40.
Snider, D. M., and Andrews
,
M. J.
, 1994, “
Rayleigh-Taylor and Shear Driven Mixing With an Unstable Thermal Stratification
,”
Phys. Fluids
6
(
10
), pp.
3324
3334
.10.1063/1.868065
41.
Snider
,
D. M.
, and
Andrews
,
M. J.
,
1996
, “
The Structure of Shear Driven Mixing With an Unstable Thermal Stratification
,”
ASME J. Fluids Eng.
,
118
(
2
), pp.
370
376
.10.1115/1.2817511
42.
Wilson
,
P. N.
, and
Andrews
,
M. J.
,
2002
, “
Spectral Measurements of Rayleigh-Taylor Mixing at Small Atwood Number
,”
Phys. Fluids
,
14
(
3
), pp.
938
945
.10.1063/1.1445418
43.
Ramaprabhu
,
P.
, and
Andrews
,
M. J.
,
2003
, “
Simultaneous Measurements of Velocity and Density in Buoyancy-Driven Mixing
,”
Exp. Fluids
,
34
(
1
), pp.
98
106
.10.1007/s00348-002-0538-0
44.
Ramaprabhu
,
P.
, and
Andrews
,
M. J.
,
2004
, “
Experimental Investigation of Rayleigh–Taylor Mixing at Small Atwood Numbers
,”
J. Fluid Mech.
,
502
, pp.
233
271
.10.1017/S0022112003007419
45.
Mueschke
,
N. J.
,
Andrews
,
M. J.
, and
Schilling
,
O.
,
2006
, “
Experimental Characterization of Initial Conditions and Spatio-Temporal Evolution of a small-Atwood-Number Rayleigh–Taylor Mixing Layer
,”
J. Fluid Mech.
,
567
, pp.
27
63
.10.1017/S0022112006001959
46.
Mueschke
,
N. J.
,
Schilling
,
O.
,
Youngs
,
D. L.
, and
Andrews
,
M. J.
,
2009
, “
Measurements of Molecular Mixing in a High Schmidt Number Rayleigh–Taylor Mixing Layer
,”
J. Fluid Mech.
,
632
, pp.
17
48
.10.1017/S0022112009006132
47.
Banerjee
,
A.
, and
Andrews
,
M. J.
,
2006
, “
Statistically Steady Measurements of Rayleigh-Taylor Mixing in a Gas Channel
,”
Phys. Fluids
,
18
(
3
), p.
035107
.10.1063/1.2185687
48.
Banerjee
,
A.
,
Kraft
,
W. N.
, and
Andrews
,
M. J.
,
2010
, “
Detailed Measurements of a Statistically Steady Rayleigh–Taylor Mixing Layer From Small to High Atwood Numbers
,”
J. Fluid Mech.
,
659
, pp.
127
190
.10.1017/S0022112010002351
49.
Akula
,
B.
,
Andrews
,
M. J.
, and
Ranjan
,
D.
,
2013
, “
Effect of Shear on Rayleigh-Taylor Mixing at Small Atwood Number
,”
Phys. Rev. E
,
87
(
3
), p.
033013
.10.1103/PhysRevE.87.033013
50.
Akula
,
B.
, and
Ranjan
,
D.
,
2016
, “
Dynamics of Buoyancy-Driven Flows at Moderately High Atwood Numbers
,”
J. Fluid Mech.
,
795
, pp.
313
355
.10.1017/jfm.2016.199
51.
Akula
,
B.
,
Suchandra
,
P.
,
Mikhaeil
,
M.
, and
Ranjan
,
D.
,
2017
, “
Dynamics of Unstably Stratified Free Shear Flows: An Experimental Investigation of Coupled Kelvin–Helmholtz and Rayleigh–Taylor Instability
,”
J. Fluid Mech.
,
816
, pp.
619
660
.10.1017/jfm.2017.95
52.
Andrews
,
M. J.
, and
Dalziel
,
S. B.
,
2010
, “
Small Atwood Number Rayleigh– Taylor Experiments
,”
Philos. Trans. R. Soc. A
,
368
(
1916
), pp.
1663
1679
.10.1098/rsta.2010.0007
53.
Taylor
,
G. I.
,
1938
, “
The Spectrum of Turbulence
,”
Proc R Soc London A
,
164
, pp.
476
490
.10.1098/rspa.1938.0032
54.
Boffetta
,
G.
, and
Mazzino
,
A.
,
2017
, “
Incompressible Rayleigh–Taylor Turbulence
,”
Annu. Rev. Fluid Mech.
,
49
(
1
), pp.
119
143
.10.1146/annurev-fluid-010816-060111
55.
Zhou
,
Y.
, and
Cabot
,
W. H.
,
2019
, “
Time-Dependent Study of Anisotropy in Rayleigh-Taylor Instability Induced Turbulent Flows With a Variety of Density Ratios
,”
Phys. Fluids
,
31
(
8
), p.
084106
.10.1063/1.5110914
56.
Cabot
,
W. H.
, and
Cook
,
A. W.
,
2006
, “
Reynolds Number Effects on Rayleigh–Taylor Instability With Implications for Type Ia Supernovae
,”
Nat. Phys.
,
2
(
8
), pp.
562
568
.10.1038/nphys361
57.
Dimonte
,
G.
, and
Tipton
,
R.
,
2006
, “
K-L Turbulence Model for the Self-Similar Growth of the Rayleigh-Taylor and Richtmyer-Meshkov Instabilities
,”
Phys. Fluids
,
18
(
8
), p.
085101
.10.1063/1.2219768
58.
Burton
,
G. C.
,
2011
, “
Study of Ultrahigh Atwood-Number Rayleigh–Taylor Mixing Dynamics Using the Nonlinear Large-Eddy Simulation Method
,”
Phys. Fluids
,
23
(
4
), p.
045106
.10.1063/1.3549931
59.
Morgan
,
B. E.
, and
Wickett
,
M. E.
,
2015
, “
Three-Equation Model for the Self-Similar Growth of Rayleigh-Taylor and Richtmyer-Meskov [sic] Instabilities
,”
Phys. Rev. E
,
91
(
4
), p.
043002
.10.1103/PhysRevE.91.043002
60.
Youngs
,
D. L.
,
2009
, “
Application of Monotone Integrated Large Eddy Simulation to Rayleigh–Taylor Mixing
,”
Philos. Trans. R. Soc. London A
,
367
(
1899
), pp.
2971
2983
.10.1098/rsta.2008.0303
61.
Youngs
,
D. L.
,
2013
, “
The Density Ratio Dependence of Self-Similar Rayleigh–Taylor Mixing
,”
Philos. Trans. R. Soc. London A
,
371
(
2003
), p.
20120173
.10.1098/rsta.2012.0173
62.
Ristorcelli
,
J. R.
, and
Clark
,
T. T.
,
2004
, “
Rayleigh–Taylor Turbulence: Self-Similar Analysis and Direct Numerical Simulations
,”
J. Fluid Mech.
,
507
, pp.
213
253
.10.1017/S0022112004008286
63.
Livescu
,
D.
,
Ristorcelli
,
J. R.
,
Petersen
,
M. R.
, and
Gore
,
R. A.
,
2010
, “
New Phenomena in Variable-Density Rayleigh–Taylor Turbulence
,”
Phys. Scr. T
,
142
, p.
014015
.10.1088/0031-8949/2010/T142/014015
64.
Cook
,
A. W.
, and
Dimotakis
,
P. E.
,
2001
, “
Transition Stages of Rayleigh–Taylor Instability Between Miscible Fluids
,”
J. Fluid Mech.
,
443
, pp.
69
99
.10.1017/S0022112001005377
65.
Chertkov
,
M.
,
2003
, “
Phenomenology of Rayleigh-Taylor Turbulence
,”
Phys. Rev. Lett.
,
91
(
11
), p.
115001
.10.1103/PhysRevLett.91.115001
66.
Pope
,
S. B.
,
2000
,
Turbulent Flows
,
Cambridge University Press
,
Cambridge, UK
.
67.
Mueschke
,
N. J.
, and
Schilling
,
O.
,
2009
, “
Investigation of Rayleigh–Taylor Turbulence and Mixing Using Direct Numerical Simulation With Experimentally Measured Initial Conditions. I. Comparison to Experimental Data
,”
Phys. Fluids
,
21
(
1
), p.
014106
.10.1063/1.3064120
68.
Cabot
,
W.
, and
Zhou
,
Y.
,
2013
, “
Statistical Measurements of Scaling and Anisotropy of Turbulent Flows Induced by Rayleigh-Taylor Instability
,”
Phys. Fluids
,
25
(
1
), p.
015107
.10.1063/1.4774338
69.
Yilmaz
,
I.
,
2020
, “
Analysis of Rayleigh–Taylor Instability at High Atwood Numbers Using Fully Implicit, Non-Dissipative, Energy-Conserving Large Eddy Simulation Algorithm
,”
Phys. Fluids
,
32
(
5
), p.
054101
.10.1063/1.5138978
70.
Kull
,
H. J.
,
1991
, “
Theory of the Rayleigh-Taylor Instability
,”
Phys. Rep.
,
206
(
5
), pp.
197
325
.10.1016/0370-1573(91)90153-D
71.
Dimonte
,
G.
,
Ramaprabhu
,
P.
,
Youngs
,
D. L.
,
Andrews
,
M. J.
, and
Rosner
,
R.
,
2005
, “
Recent Advances in the Turbulent Rayleigh-Taylor Instability
,”
Phys. Plasm.
,
12
(
5
), p.
056301
.10.1063/1.1871952
72.
Dimonte
,
G.
,
2004
, “
Dependence of Turbulent Rayleigh-Taylor Instability on Initial Perturbations
,”
Phys. Rev. E
,
69
(
5
), p.
056305
.10.1103/PhysRevE.69.056305
73.
Ramaprabhu
,
P.
,
Dimonte
,
G.
, and
Andrews
,
M. J.
,
2005
, “
A Numerical Study of the Influence of Initial Perturbations on the Turbulent Rayleigh–Taylor Instability
,”
J. Fluid Mech.
,
536
, pp.
285
319
.10.1017/S002211200500488X
74.
Ramaprabhu
,
P.
, and
Andrews
,
M. J.
,
2004
, “
On the Initialization of Rayleigh-Taylor Simulations
,”
Phys. Fluids
,
16
(
8
), pp.
L59
L62
.10.1063/1.1765171
75.
Dimotakis
,
P. E.
,
2000
, “
The Mixing Transition in Turbulence
,”
J. Fluid Mech.
,
409
, pp.
69
98
.10.1017/S0022112099007946
76.
Zhou
,
Y.
,
Remington
,
B. A.
,
Robey
,
H. F.
,
Cook
,
A. W.
,
Glendinning
,
S. G.
,
Dimits
,
A.
,
Buckingham
,
A. C.
,
Zimmerman
,
G. B.
,
Burke
,
E. W.
,
Peyser
,
T. A.
,
Cabot
,
W.
, and
Eliason
,
D.
,
2003
, “
Progress in Understanding Turbulent Mixing Induced by Rayleigh-Taylor and Richtmyer-Meshkov Instabilities
,”
Phys. Plasm.
,
10
(
5
), pp.
1883
1896
.10.1063/1.1560923
77.
Zhou
,
Y.
,
Clark
,
T. T.
,
Clark
,
D. S.
,
Glendinning
,
S. G.
,
Skinner
,
M. A.
,
Huntington
,
C. M.
,
Hurricane
,
O. A.
,
Dimits
,
A. M.
, and
Remington
,
B. A.
,
2019
, “
Turbulent Mixing and Transition Criteria of Flows Induced by Hydrodynamic Instabilities
,”
Phys. Plasmas
,
26
(
8
), p.
080901
.10.1063/1.5088745
78.
Elbaz
,
Y. S.
, and
Shvarts
,
D.
,
2018
, “
Modal Mean Field Self-Similar Solutions to the Asymptotic Evolution of Rayleigh-Taylor and Richtmyer-Meshkov Instabilities and Its Dependence on the Initial Conditions
,”
Phys. Plasmas
,
25
(
6
), p.
062126
.10.1063/1.5031922
79.
Canfield
,
J.
,
Denissen
,
N.
,
Francois
,
M.
,
Gore
,
R.
,
Rauenzahn
,
R.
,
Reisner
,
J.
, and
Shkoller
,
S.
,
2020
, “
A Comparison of Interface Growth Models Applied to Rayleigh-Taylor and Richtmyer-Meshkov Instabilities
,”
ASME J. Fluids Eng.
,
142
(
12
), p.
121108
.10.1115/1.4048341
80.
Llor
,
A.
,
2006
,
Statistical Hydrodynamic Models for Developed Mixing Instability Flows
(Lecture Notes in Physics, Vol.
681
),
Springer-Verlag
,
New York
.
81.
Peters
,
N.
,
2000
,
Turbulent Combustion
(Cambridge Monographs on Mechanics),
Cambridge University Press
,
Cambridge, UK
.
82.
Poinsot
,
T.
, and
Veynante
,
D.
,
2012
,
Theoretical and Numerical Combustion
, 3rd ed.,
R. T. Edwards
,
Philadelphia, PA
.
83.
Wyngaard
,
J. C.
,
2010
,
Turbulence in the Atmosphere
,
Cambridge University Press
,
Cambridge, UK
.
84.
Vallis
,
G. K.
,
2017
,
Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation
,
Cambridge University Press
,
Cambridge, UK
.
85.
Schwarzkopf
,
J. D.
,
Livescu
,
D.
,
Gore
,
R. A.
,
Rauenzahn
,
R. M.
, and
Ristorcelli
,
J. R.
,
2011
, “
Application of a Second-Moment Closure Model to Mixing Processes Involving Multicomponent Miscible Fluids
,”
J. Turbul.
,
12
(
49
), pp.
1
35
.10.1080/14685248.2011.633084
86.
Schwarzkopf
,
J. D.
,
Livescu
,
D.
,
Baltzer
,
J. R.
,
Gore
,
R. A.
, and
Ristorcelli
,
J. R.
,
2016
, “
A Two-Length Scale Turbulence Model for Single-Phase Multi-Fluid Mixing
,”
Flow Turbul. Combust.
,
96
(
1
), pp.
1
43
.10.1007/s10494-015-9643-z
87.
Zikanov
,
O.
,
2019
,
Essential Computational Fluid Dynamics
, 2nd ed.,
Wiley
,
Hoboken, NJ
.
88.
Livescu
,
D.
,
2020
, “
Turbulence With Large Thermal and Compositional Density Variations
,”
Annu. Rev. Fluid Mech.
,
52
(
1
), pp.
309
341
.10.1146/annurev-fluid-010719-060114
89.
Joseph
,
D. D.
,
1990
, “
Fluid Dynamics of Two Miscible Liquids With Diffusion and Gradient Stresses
,”
Eur. J. Mech. B
,
9
(
6
), pp.
565
596
.
90.
Joseph
,
D. D.
,
Funada
,
T.
, and
Wang
,
J.
,
2007
,
Potential Flows of Viscous and Viscoelastic Liquids
(Cambridge Aerospace Series),
Cambridge University Press
,
Cambridge, UK
.
91.
Williams
,
F. A.
,
1985
,
Combustion Theory: The Fundamental Theory of Chemically Reacting Flow Systems
, 2nd ed.,
Addison-Wesley
,
Menlo Park, CA
.
92.
Livescu
,
D.
,
2013
, “
Numerical Simulations of Two-Fluid Turbulent Mixing at Large Density Ratios and Applications to the Rayleigh–Taylor Instability
,”
Philos. Trans. R. Soc. London A
,
371
(
2003
), p.
20120185
.10.1098/rsta.2012.0185
93.
Drikakis
,
D.
, and
Rider
,
W. J.
,
2005
,
High-Resolution Methods for Incompressible and Low-Speed Flows
(Computational Fluid and Solid Mechanics Series),
Springer-Verlag
,
New York
.
94.
Grinstein
,
F. F.
,
Margolin
,
L. G.
, and
Rider
,
W. J.
, (eds.),
2007
,
Implicit Large Eddy Simulation: Computing Turbulent Fluid Dynamics
,
Cambridge University Press
,
Cambridge, UK
.
95.
Grinstein
,
F. F.
, and
Fureby
,
C.
,
2002
, “
Recent Progress on MILES for High Reynolds Number Flows
,”
ASME J. Fluids Eng.
,
124
(
4
), pp.
848
861
.10.1115/1.1516576
96.
Margolin
,
L. G.
, and
Rider
,
W. J.
,
2002
, “
A Rationale for Implicit Turbulence Modelling
,”
Int. J. Numer. Methods Fluids
,
39
(
9
), pp.
821
841
.10.1002/fld.331
97.
Margolin
,
L. G.
, and
Rider
,
W. J.
,
2005
, “
The Design and Construction of Implicit LES Methods
,”
Int. J. Numer. Methods Fluids
,
47
(
10–11
), pp.
1173
1179
.10.1002/fld.862
98.
Drikakis
,
D.
,
Grinstein
,
F.
, and
Youngs
,
D.
,
2005
, “
On the Computation of Instabilities and Symmetry-Breaking in Fluid Mechanics
,”
Prog. Aerosp. Sci.
,
41
(
8
), pp.
609
641
.10.1016/j.paerosci.2005.10.001
99.
Orlandi
,
P.
,
2000
,
Fluid Flow Phenomena: A Numerical Toolkit
(Fluid Mechanics and Its Applications, Vol.
55
),
Kluwer Academic
,
Dordrecht, Holland
.
100.
Lesieur
,
M.
,
Métais
,
O.
, and
Comte
,
P.
,
2005
,
Large-Eddy Simulations of Turbulence
,
Cambridge University Press
,
Cambridge, UK
.
101.
Aldama
,
A. A.
,
1990
,
Filtering Techniques for Turbulent Flow Simulation
(Lecture Notes in Engineering, Vol.
56
),
Springer-Verlag
,
New York
.
102.
Youngs
,
D. L.
,
2004
, “
Effect of Initial Conditions on Self-Similar Turbulent Mixing
,”
In Proceedings of the Ninth International Workshop on the Physics of Compressible Turbulent Mixing
,
S. B
,
Dalziel
, ed.,
University of Cambridge
,
Cambridge, UK
, July 19–23, pp.
1514
1545
.
103.
Lee
,
H.
,
Jin
,
H.
,
Yu
,
Y.
, and
Glimm
,
J.
,
2008
, “
On Validation of Turbulent Mixing Simulations for Rayleigh–Taylor Instability
,”
Phys. Fluids
,
20
(
1
), p.
012102
.10.1063/1.2832775
104.
Lim
,
H.
,
Yu
,
Y.
,
Glimm
,
J.
,
Li
,
X.-L.
, and
Sharp
,
D. H.
,
2010
, “
Subgrid Models for Mass and Thermal Diffusion in Turbulent Mixing
,”
Phys. Scr. T
,
T142
, p.
014062
.10.1088/0031-8949/2010/T142/014062
105.
Lim
,
H.
,
Iwerks
,
J.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2010
, “
Nonideal Rayleigh–Taylor Mixing
,”
Proc Natl Acad Sci
,
107
, p.
12786
.10.1073/pnas.1002410107
106.
Lim
,
H.
,
Iwerks
,
J.
,
Yu
,
Y.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2010
, “
Verification and Validation of a Method for the Simulation of Turbulent Mixing
,”
Phys. Scr. T
,
T142
, p.
014014
.10.1088/0031-8949/2010/T142/014014
107.
Kaman
,
T.
,
Melvin
,
J.
,
Rao
,
P.
,
Kaufman
,
R.
,
Lim
,
H.
,
Yu
,
Y.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2013
, “
Recent Progress in Turbulent Mixing
,”
Phys. Scr. T
,
T155
, p.
014051
.10.1088/0031-8949/2013/T155/014051
108.
Glimm
,
J.
,
Sharp
,
D. H.
,
Kaman
,
T.
, and
Lim
,
H.
,
2013
, “
New Directions for Rayleigh–Taylor Mixing
,”
Philos. Trans. R. Soc. London A
,
371
(
2003
), p.
20120183
.10.1098/rsta.2012.0183
109.
Andrews
,
M. J.
,
Youngs
,
D. L.
,
Livescu
,
D.
, and
Wei
,
T.
,
2014
, “
Computational Studies of Two-Dimensional Rayleigh-Taylor Driven Mixing for a Tilted-Rig
,”
ASME J. Fluid Eng.
,
136
(
9
), p.
091212
.10.1115/1.4027587
110.
Dimonte
,
G.
,
Ramaprabhu
,
P.
, and
Andrews
,
M.
,
2007
, “
Rayleigh-Taylor Instability With Complex Acceleration History
,”
Phys. Rev. E
,
76
(
4
), p.
046313
.10.1103/PhysRevE.76.046313
111.
Lawrie
,
A.
, and
Dalziel
,
S.
,
2011
, “
Turbulent Diffusion in Tall Tubes—I: Models for Rayleigh-Taylor Instability
,”
Phys. Fluids
,
23
(
8
), p.
085109
.10.1063/1.3614477
112.
Dalziel
,
S. B.
,
Patterson
,
M. D.
,
Caulfield
,
C. P.
, and
Coomaraswamy
,
I. A.
,
2008
, “
Mixing Efficiency in High-Aspect-Ratio Rayleigh–Taylor Experiments
,”
Phys. Fluids
,
20
(
6
), p.
065106
.10.1063/1.2936311
113.
Williams
,
R. J. R.
,
2017
, “
Rayleigh–Taylor Mixing Between Density Stratified Layers
,”
J. Fluid Mech.
,
810
, pp.
584
602
.10.1017/jfm.2016.740
114.
Lawrie
,
A.
, and
Dalziel
,
S.
,
2011
, “
Rayleigh–Taylor Mixing in an Otherwise Stable Stratification
,”
J. Fluid Mech.
,
688
, pp.
507
527
.10.1017/jfm.2011.398
115.
Davies Wykes
,
M. S.
, and
Dalziel
,
S. B.
,
2014
, “
Efficient Mixing in Stratified Flows: Experimental Study of a Rayleigh–Taylor Unstable Interface Within an Otherwise Stable Stratification
,”
J. Fluid Mech.
,
756
, pp.
1027
1057
.10.1017/jfm.2014.308
116.
Mueschke
,
N. J.
, and
Schilling
,
O.
,
2009
, “
Investigation of Rayleigh–Taylor Turbulence and Mixing Using Direct Numerical Simulation With Experimentally Measured Initial Conditions—II: Dynamics of Transitional Flow and Mixing Statistics
,”
Phys. Fluids
,
21
(
1
), p.
014107
.10.1063/1.3064121
117.
Mueschke
,
N. J.
,
2008
, “
Experimental and Numerical Study of Molecular Mixing Dynamics in Rayleigh–Taylor Unstable Flows
,” Ph.D. Thesis,
Texas A & M University
,
College Station, TX
.
118.
Abarzhi
,
S. I.
,
2010
, “
Review of Theoretical Modelling Approaches of Rayleigh–Taylor Instabilities and Turbulent Mixing
,”
Philos. Trans. R. Soc. A
,
368
(
1916
), pp.
1809
1828
.10.1098/rsta.2010.0020
119.
Haan
,
S. W.
,
1989
, “
Onset of Nonlinear Saturation for Rayleigh-Taylor Growth in the Presence of a Full Spectrum of Modes
,”
Phys. Rev. A
,
39
(
11
), pp.
5812
5825
.10.1103/PhysRevA.39.5812
120.
Haan
,
S. W.
,
1991
, “
Weakly Nonlinear Hydrodynamic Instabilities in Inertial Fusion
,”
Phys. Fluids B
,
3
(
8
), pp.
2349
2355
.10.1063/1.859603
121.
Cheng
,
B.
,
2009
, “
Review of Turbulent Mixing Models
,”
Acta Math Sci
,
29
(
6
), pp.
1703
1720
.10.1016/S0252-9602(10)60012-4
122.
Cheng
,
B.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2020
, “
The αs and θs in Rayleigh–Taylor and Richtmyer–Meshkov Instabilities
,”
Phys. D
,
404
, p.
132356
.10.1016/j.physd.2020.132356
123.
Schilling
,
O.
,
2020
, “
A Buoyancy–Shear–Drag-Based Turbulence Model for Rayleigh–Taylor, Reshocked Richtmyer–Meshkov, and Kelvin–Helmholtz Mixing
,”
Phys. D
,
402
, p.
132238
.10.1016/j.physd.2019.132238
124.
Cheng
,
B.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2000
, “
Density Dependence of Rayleigh–Taylor and Richtmyer–Meshkov Mixing Fronts
,”
Phys. Lett. A
,
268
(
4–6
), pp.
366
374
.10.1016/S0375-9601(00)00204-8
125.
Cheng
,
B.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2002
, “
Dynamical Evolution of Rayleigh-Taylor and Richtmyer-Meshkov Mixing Fronts
,”
Phys. Rev. E
,
66
(
3
), p.
036312
.10.1103/PhysRevE.66.036312
126.
Dimonte
,
G.
,
2000
, “
Spanwise Homogeneous Buoyancy-Drag Model for Rayleigh–Taylor Mixing and Experimental Evaluation
,”
Phys. Plasm
,
7
(
6
), pp.
2255
2269
.10.1063/1.874060
127.
Davies
,
R. M.
, and
Taylor
,
G. I.
,
1950
, “
The Mechanics of Large Bubbles Rising Through Extended Liquids and Through Liquids in Tubes
,”
Proc R Soc London A
,
200
(
1062
), pp.
375
390
.10.1098/rspa.1950.0023
128.
Chen
,
C.-J.
, and
Jaw
,
S.-Y.
,
1998
,
Fundamentals of Turbulence Modeling. Combustion: An International Series
,
Taylor and Francis
,
London
.
129.
Durbin
,
P. A.
, and
Pettersson Reif
,
B. A.
,
2010
,
Statistical Theory and Modeling for Turbulent Flows
, 2nd ed.,
Wiley
,
Hoboken, NJ
.
130.
Chassaing
,
P.
,
Antonia
,
R. A.
,
Anselmet
,
F.
,
Joly
,
L.
, and
Sarkar
,
S.
,
2002
,
Variable Density Fluid Turbulence
(Fluid Mechanics and Its Applications, Vol.
69
),
Kluwer Academic
,
Dordrecht, Holland
.
131.
Gatski
,
T. B.
, and
Bonnet
,
J.-P.
,
2013
,
Compressibility, Turbulence and High Speed Flow
, 2nd ed.,
Academic Press
,
Boston, MA
.
132.
Morán-López
,
J. T.
, and Schilling, O.,
2013
, “
Multicomponent Reynolds-Averaged Navier–Stokes Simulations of Reshocked Richtmyer–Meshkov Instability-Induced Mixing
,”
High Energy Density Phys.
,
9
(
1
), pp.
112
121
.10.1016/j.hedp.2012.11.001
133.
Morán-López
,
J. T.
, and Schilling, O.,
2014
, “
Multi-Component Reynolds-Averaged Navier–Stokes Simulations of Richtmyer–Meshkov Instability and Mixing Induced by Reshock at Different Times
,”
Shock Waves
,
24
(
3
), pp.
325
343
.10.1007/s00193-013-0483-2
134.
Ruan
,
Y.-C.
,
Zhang
,
Y.-S.
,
Tian
,
B.-L.
, and
Zhang
,
X.-T.
,
2020
, “
Density-Ratio-Invariant Mean-Species Profile of Classical Rayleigh-Taylor Mixing
,”
Phys. Rev. Fluids
,
5
(
5
), p.
054501
.10.1103/PhysRevFluids.5.054501
135.
Ofer
,
D.
,
Alon
,
U.
,
Shvarts
,
D.
,
McCrory
,
R. L.
, and
Verdon
,
C. P.
,
1996
, “
Modal Model for the Nonlinear Multimode Rayleigh-Taylor Instability
,”
Phys. Plasm
,
3
(
8
), pp.
3073
3090
.10.1063/1.871655
136.
Granero-Belinchón
,
R.
, and
Shkoller
,
S.
,
2017
, “
A Model for Rayleigh–Taylor Mixing and Interface Turnover
,”
Multiscale Model. Simul.
,
15
(
1
), pp.
274
308
.10.1137/16M1083463
137.
Cheng
,
B.
,
Glimm
,
J.
, and
Sharp
,
D. H.
,
2002
, “
A Three-Dimensional Renormalization Group Bubble Merger Model for Rayleigh–Taylor Mixing
,”
Chaos
,
12
(
2
), pp.
267
274
.10.1063/1.1460942
138.
Ramshaw
,
J. D.
,
1998
, “
Simple Model for Linear and Nonlinear Mixing at Unstable Fluid Interfaces With Variable Acceleration
,”
Phys. Rev. E
,
58
(
5
), pp.
5834
5840
.10.1103/PhysRevE.58.5834
139.
Mikaelian
,
K. O.
,
1989
, “
Turbulent Mixing Generated by Rayleigh-Taylor and Richtmyer-Meshkov Instabilities
,”
Phys. D
,
36
(
3
), pp.
343
357
.10.1016/0167-2789(89)90089-4
140.
Mikaelian
,
K. O.
,
1990
, “
Turbulent Energy at Accelerating and Shocked Interfaces
,”
Phys. Fluids A
,
2
(
4
), pp.
592
598
.10.1063/1.857759
141.
Alon
,
U.
,
Hecht
,
J.
,
Ofer
,
D.
, and
Shvarts
,
D.
,
1995
, “
Power Laws and Similarity of Rayleigh-Taylor and Richtmyer-Meshkov Mixing Fronts at All Density Ratios
,”
Phys. Rev. Lett.
,
74
(
4
), pp.
534
537
.10.1103/PhysRevLett.74.534
142.
Dimonte
,
G.
, and
Schneider
,
M.
,
1997
, “
Turbulent Richtmyer-Meshkov Instability Experiments With Strong Radiatively Driven Shocks
,”
Phys. Plasm.
,
4
(
12
), pp.
4347
4357
.10.1063/1.872597
143.
Zhang
,
Y.
,
He
,
Z.
,
Gao
,
F.
,
Li
,
X.
, and
Tian
,
B.
,
2016
, “
Evolution of Mixing Width Induced by General Rayleigh-Taylor Instability
,”
Phys. Rev. E
,
93
(
6
), p.
063102
.10.1103/PhysRevE.93.063102
144.
Kokkinakis
,
I. W.
,
Drikakis
,
D.
,
Youngs
,
D. L.
, and
Williams
,
R. J. R.
,
2015
, “
Two-Equation and Multi-Fluid Turbulence Models for Rayleigh–Taylor Mixing
,”
Int. J. Heat Fluid Flow
,
56
, pp.
233
250
.10.1016/j.ijheatfluidflow.2015.07.017
145.
Zhou
,
Y.
,
Zimmerman
,
G. B.
, and
Burke
,
E. W.
,
2002
, “
Formulation of a Two-Scale Transport Scheme for the Turbulent Mix Induced by Rayleigh-Taylor and Richtmyer-Meshkov Instabilities
,”
Phys. Rev. E
,
65
(
5
), p.
056303
.10.1103/PhysRevE.65.056303
146.
Kucherenko
,
Y. A.
,
Shibarshov
,
L. I.
,
Chitaikin
,
V. I.
,
Balabin
,
S. I.
, and
Pylaev
,
A. P.
,
1991
, “
Experimental Study of the Gravitational Turbulent Mixing Self-Similar Mode
,”
In Proceedings of the Third International Workshop on the Physics of Compressible Turbulent Mixing
,
D.
Besnard
,
J.
Haas
,
P.
Holstein
, and
B.
Sitt
, eds.,
Royaumont, France
, June 17–19, pp.
427
454
.
147.
Denissen
,
N. A.
,
Rollin
,
B.
,
Reisner
,
J. M.
, and
Andrews
,
M. J.
,
2014
, “
The Tilted Rocket Rig: A Rayleigh–Taylor Test Case for RANS Models
,”
ASME J. Fluid Eng.
,
136
, p.
091301
.10.1115/1.4027776
148.
Kokkinakis
,
I. W.
,
Drikakis
,
D.
, and
Youngs
,
D. L.
,
2019
, “
Modeling of Rayleigh-Taylor Mixing Using Single-Fluid Models
,”
Phys. Rev. E
,
99
(
1
), p.
013104
.10.1103/PhysRevE.99.013104
149.
Wachtor
,
A. J.
,
Mocko
,
V.
,
Williams
,
D. J.
,
Goertz
,
M. P.
, and
Jebrail
,
F. F.
,
2013
, “
Buoyancy Driven Mixing of Miscible Fluids by Volumetric Energy Deposition of Microwaves
,”
J. Microwave Power Electromagn. Energy
,
47
(
3
), pp.
210
223
.10.1080/08327823.2013.11689859
150.
Wachtor
,
A. J.
,
Mocko
,
V.
,
Jebrail
,
F. F.
, and
Andrews
,
M. J.
,
2015
, “
On Buoyancy Driven Mixing by Volumetric Microwave Energy Deposition
,”
Int. J. Heat Mass Transfer
,
86
, pp.
443
454
.10.1016/j.ijheatmasstransfer.2015.01.112
151.
Gréa
,
B.
,
2013
, “
The Rapid Acceleration Model and the Growth Rate of a Turbulent Mixing Zone Induced by Rayleigh-Taylor Instability
,”
Phys. Fluids
,
25
(
1
), p.
015118
.10.1063/1.4775379
152.
Hillier
,
A.
,
2020
, “
Self-Similar Solutions of Asymmetric Rayleigh-Taylor Mixing
,”
Phys. Fluids
,
32
(
1
), p.
015103
.10.1063/1.5130893
153.
Andrews
,
M. J.
,
1984
, “
The k-ε Model Applied to the Development of Rayleigh–Taylor Instability
,”
Imperial College of Science and Technology
,
London
, Report No. PDR/CFDU IC/13.
154.
Spitz
,
P. B.
, and
Haas
,
J.
,
1991
, “
Numerical Calibration of Rayleigh-Taylor Induced Turbulent Flows With a k-ε Mix Model
,”
In Proceedings of the Third International Workshop on the Physics of Compressible Turbulent Mixing
,
D.
Besnard
,
J.
Haas
,
P.
Holstein
, and
B.
Sitt
, eds.,
Royaumont, France
, June 17-19, pp.
511
521
.
155.
Snider
,
D. M.
, and
Andrews
,
M. J.
,
1996
, “
The Simulation of Mixing Layers Driven by Compound Buoyancy and Shear
,”
ASME J. Fluids Eng.
,
118
(
2
), pp.
370
376
.10.1115/1.2817388
156.
Besnard
,
D.
,
Harlow
,
F. H.
,
Rauenzahn
,
R. M.
, and
Zemach
,
C.
,
1996
, “
Spectral Transport Model for Turbulence
,”
Theor. Comput. Fluid Dyn.
,
8
(
1
), pp.
1
35
.10.1007/BF00312400
157.
Wilson
,
P.
,
Andrews
,
M.
, and
Harlow
,
F.
,
1999
, “
Spectral Nonequilibrium in a Turbulent Mixing Layer
,”
Phys. Fluids
,
11
(
8
), pp.
2425
2433
.10.1063/1.870103
158.
Banerjee
,
A.
,
Gore
,
R. A.
, and
Andrews
,
M. J.
,
2010
, “
Development and Validation of a Turbulent-Mix Model for Variable-Density and Compressible Flows
,”
Phys. Rev. E
,
82
(
4
), p.
046309
.10.1103/PhysRevE.82.046309
159.
Stalsberg-Zarling
,
K.
, and
Gore
,
R.
,
2011
, “
The BHR2 Turbulence Model: Incompressible Isotropic Decay, Rayleigh-Taylor, Kelvin-Helmholtz and Homogeneous Variable Density Turbulence
,”
Los Alamos National Laboratory
,
Los Alamos, NM
, Report No. LA-UR-11-04773.
160.
Schilling
,
O.
, and
Mueschke
,
N. J.
,
2010
, “
Analysis of Turbulent Transport and Mixing in Transitional Rayleigh–Taylor Unstable Flow Using Direct Numerical Simulation Data
,”
Phys. Fluids
,
22
(
10
), p.
105102
.10.1063/1.3484247
161.
Schilling
,
O.
, and
Mueschke
,
N. J.
,
2017
, “
Turbulent Transport and Mixing in Transitional Rayleigh-Taylor Unstable Flow: A Priori Assessment of Gradient-Diffusion and Similarity Modeling
,”
Phys. Rev. E
,
96
(
6
), p.
063111
.10.1103/PhysRevE.96.063111