## Abstract

Previous research on self-similar mixing caused by Rayleigh–Taylor (RT) instability is summarized and a recent series of high resolution large eddy simulations is described. Mesh sizes of approximately 2000 ×1000 × 1000 are used to investigate the properties of high Reynolds number self-similar RT mixing at a range of density ratios from 1.5 : 1 to 20 : 1. In some cases, mixing evolves from ‘small random perturbations’. In other cases, random long wavelength perturbations (*k*^{−3} spectrum) are added to give self-similar mixing at an enhanced rate, more typical of that observed in experiments. The properties of the turbulent mixing zone (volume fraction distributions, turbulence kinetic energy, molecular mixing parameter, etc.) are related to the RT growth rate parameter, *α*. Comparisons are made with experimental data on the internal structure and the asymmetry of the mixing zone (spike distance/bubble distance). The main purpose of this series of simulations is to provide data for calibration of engineering models (e.g. Reynolds-averaged Navier–Stokes models). It is argued that the influence of initial conditions is likely to be significant in most applications and the implications of this for engineering modelling are discussed.

## 1. Introduction

Rayleigh–Taylor (RT) instability is of current concern for inertial confinement fusion (ICF); see for example, Amendt *et al.* [1]. Such instabilities can degrade the performance of ICF implosions, where high density shells are decelerated by lower density thermonuclear fluid. In these applications and in simpler laboratory experiments, the Reynolds number is likely to be high. Turbulent mixing will then occur. The simplest situation in which RT instability occurs consists of fluid with density *ρ*_{1} resting initially above fluid with density *ρ*_{2}<*ρ*_{1} in a gravitational field *g*. Experiments with relatively low viscosity, surface tension and initial perturbations have shown that the dominant length scale increases as mixing evolves. If mixing is self-similar then dimensional reasoning suggests the length scale should be proportional to *gt*^{2}. The experiments indicate that the depth to which the mixing zone penetrates the denser fluid, generally known as the bubble distance, is given by
1.1where *α* varies slightly with Atwood number, *A*, and is typically in the range 0.04–0.08. The distance to which the mixing zone penetrates the lighter fluid is referred to as the spike distance, *h*_{s}. The ratio *h*_{s}/*h*_{b} is a slowly increasing function of the density ratio *ρ*_{1}/*ρ*_{2}. Three-dimensional simulation of this self-similar mixing problem is the subject of this paper.

For complex applications, such as ICF implosions, some form of engineering modelling, e.g. Reynolds-averaged Navier–Stokes, is considered necessary at present and such modelling requires calibration to determine model constants. Self-similar RT mixing of incompressible fluids at a plane boundary is one of the key test cases for model calibration. Some data for model calibration are available from experiments. However, there are few detailed measurements at high Atwood number (more than 1/2). Hence, it is proposed that the results of three-dimensional simulations, with experimental validation where possible, should be used to provide the data needed for model calibration.

Direct numerical simulation (DNS) is feasible at moderate Reynolds number and is essential for assessing the influence of Reynolds number and Schmidt number. However, as many of the applications are in the high-Reynolds number regime, model coefficients need to be set for this case. It is argued here that high-resolution large eddy simulation (LES) is the most computationally efficient way a calculating the high-Reynolds number behaviour and that this approximation was needed to make the present series of simulations feasible. As the problem of interest has an initial discontinuity, monotone integrated LES (MILES) is strongly favoured. Mixing of miscible fluids is considered and the Reynolds number is assumed to be high enough for the effect of Schmidt number to be unimportant, that is, the flow is beyond the ‘mixing transition’ as defined by Dimotakis [2].

Section 2 reviews some of the previous research on RT mixing, with the emphasis on the self-similar behaviour. More extensive reviews are given by Abarzhi [3] and Abarzhi & Rosner [4]. Sections 3 summarizes the numerical method used here and give details of the calculations performed. The results needed for engineering model calibration are discussed in §4 and concluding remarks are given in §5.

## 2. Summary of previous work on Rayleigh–Taylor mixing at a plane boundary

Much of the early work in the USA and UK, after G.I. Taylor's 1950 paper, focused on the growth of a single sinusoidal mode by RT instability, following its development from the linear (small amplitude) to nonlinear (large-amplitude) regimes. Early research is summarized in the review paper by Sharp [5] and an up-to-date summary of single-mode theory is included in Abarzhi [6]. However, in the early papers, there were some references to self-similar growth as described by equation (1.1). Birkhoff [7] applied a model for growth of a single mode to the case of random perturbations by making the assumption that perturbations of order 0.001*λ* are present at all wavelengths *λ* and then using a model for the growth of a single mode to estimate the growth of the dominant wavelength, *λ*_{d}, as a function of time. For Atwood number=1, this gave
2.1Fermi & von Neumann [8] did consider the case *A*≠1 and, by the use of a simple energy conservation argument, obtained *h*_{s}∼((*ρ*_{1}−*ρ*_{2})/*ρ*_{1})*gt*^{2} for single-mode spike growth. By taking the single-mode theory at *A*=1 as a starting point, Sharp & Wheeler [9], see also Sharp [5], argued that if initial perturbations were low then at late time, a bubble merger process would dominate and give self-similar growth with
2.2A very different approach was used by Belen'kii & Fradkin [10], in work published in 1965, who modelled turbulent mixing caused by RT instability by use of an analogy to shear flow mixing. According to V. B. Rozanov (2010, private communication), this paper was based on research carried out much earlier, *ca* 1950. In the instability process, potential energy is converted to turbulence kinetic energy, *k*. Mixing then arises from turbulent diffusion with diffusion coefficient , where ℓ_{t}=*αL* and *L* is the mixing zone width. The result obtained is
2.3It can be argued that current work of self-similar RT mixing is based on a combination of these two different approaches.

Experiments performed to verify equation (2.3) were described by Anuchina *et al.* [11]. However, it was not until the experiments of Read [12] that consistent results for the self-similar growth rate were obtained, indicating *α*∼0.06 in equation (1.1). Further experimental confirmation of equation (1.1) at a wider range of density ratios was given by Youngs [13] and Dimonte & Schnieder [14,15]. The latter paper indicated *α*∼0.05. A summary of observed growth rates is given in Dimonte *et al.* [16].

It is evident now that self-similar mixing can be achieved in two different ways. Firstly, case A, if the initial perturbation consists of a random combination of short wavelength perturbations, larger and larger scales will eventually evolve from nonlinear interactions between smaller scales. It seems reasonable to assume that loss of memory of the initial conditions will then occur, implying self-similar mixing. Alternatively, case B, the large scales present at late time may evolve directly from initial perturbations at these scales. The distinction between these two cases was also recognized by Sharp & Wheeler [9]. If the initial random perturbation spectrum is such that ‘amplitude ∝ wavelength’ at all scales up to the size of the experiment, then as pointed out by Inogamov [17], true self-similar mixing may occur at an enhanced rate. Note that this assumption was also on the basis of Birkhoff's [7] result, equation (2.1). In more precise terms, the power spectrum, *P*(*k*), for the initial perturbation height and the standard deviation of the surface, *σ*, in case B are given by
2.4After the experiments of Read [12], Youngs [13] and Dimonte & Schneider [14] were performed, it was thought that case A corresponded to the lower end of the observed range of values of *α*. However, when high-resolution three-dimensional simulation became feasible, calculations for this case indicated very low values of *α*∼0.03 [18]. This conclusion was supported by Dimonte *et al.* [16] in which results for seven different MILES techniques (including the TURMOIL method summarized in §3), using a mesh size of 512×256×256 all indicated similar low values, *α*=0.025±0.003. There was a high degree of consistency between the different methods used, for both the overall growth rate and the internal structure of the mixing zone. The definition used for *α* was the late-time value of d*h*_{b}/d*X* where *X*=*Agt*^{2}. The results suggested that the observed growth rates were all influenced by initial conditions. Most of the experiments described so far used immiscible fluids, whereas the simulations assume miscible mixing. It is possible that some measured values of *α* may have been influenced by surface tension. It could be argued that low values of surface tension may suppress fine-scale mixing and thereby increase *α* by increasing the effective Atwood number. However, it will be shown in §4 that very low levels of long wavelength perturbations significantly increase *α*. Hence, this is a likely explanation of the higher observed values of *α*. By use of a simple model, Dimonte [19] argued that *α* should have a weak (logarithmic) dependence on initial perturbation amplitudes. This explains why *α* does not vary greatly within a series of experiments using the same apparatus.

Very high-resolution DNS results, 3072^{3} meshes, were given by Cabot & Cook [20], corresponding to case A. In this paper, an alternative way of calculating *α* as a function of time, proposed by Ristorcelli & Clark [21], was used: . This showed that *α* was high early in the simulation and then dropped to a low value after the ‘mixing transition’, as defined by Dimotakis [2], occurred. The value of *α* then slowly increased to approximately 0.025 by the end of the simulation, close to the values obtained from the MILES, when the Reynolds number, , reached 3×10^{4}. The recent DNS of Livescu *et al.* [22] at Atwood numbers of 0.04 and 0.5 also indicates *α*∼0.025 and there seems little doubt now that the value of *α* for case A is much less than the observed values.

The available results for spike/bubble asymmetry at high Atwood number are reviewed by Dimonte & Schneider [15]. At very high Atwood number, *A*>3/4, the results obtained by Dimonte & Schneider showed greater asymmetry than the earlier results of Youngs [13] and Kucherenko *et al.* [23]. Volume fraction profiles at moderate Atwood number are given by Dimonte & Schneider [15] and Kucherenko *et al.* [23]. However, very limited measurement of the detailed structure of the mixing zone was feasible for these experiments which mainly used immiscible fluids. For experiments using miscible liquids at very low Atwood number, much more detailed measurements are available [24–27]. Of particular use for validation of the current series of three-dimensional simulations are the measurements of the molecular parameter, defined in §4.

The results discussed in this paper are for the case of *g* being constant. However, for the ICF problem *g* will vary with time. As noted by Llor [28], self-similar RT mixing can also occur for the case *g*∝*t*^{p} where positive and some negative values are permissible. Poujade & Peybernes [29] gave estimates of the coefficient *α*_{p} in the equation *h*_{b}=*α*_{p}*Ag*(*t*)*t*^{p}, for *p* in the range 0–3. These additional self-similar cases should significantly add to the understanding of more complex RT problems.

## 3. The current series of high-resolution Rayleigh–Taylor simulations

The results shown here have been obtained by using the TURMOIL Lagrange-remap hydrocode which calculates the mixing of compressible fluids. The numerical method solves the Euler equations plus advection equations for fluid mass fractions, *m*_{r}
3.1
3.2
and
3.3where *u*_{i}, *ρ*, *e*, *p* denote fluid velocity, density, specific internal energy and pressure, respectively.

A perfect gas equation of state is used. The ratio of specific heats, *γ*=5/3, is the same for each fluid in the calculations shown here. This avoids some spurious numerical oscillations which may occur if different values of *γ* were used. As low Mach number calculations are used here to model incompressible mixing, the value of *γ* does not affect the results. In the simulations, ‘sub-grid’ dissipation of both velocity and density fluctuations occurs at a scale of the order of the mesh size. Throughout the rest of this paper, it is assumed that the problem of interest is the mixing of miscible fluids at very high Reynolds number, that is, a regime is reached where the dissipation rates are determined by the Kolmgorov cascade process and insensitive to the values of viscosity and diffusivity.

The Lagrange-remap method was first used for RT mixing by Youngs [30]. More details of the method, its application to RT and Richtmyer–Meshkov mixing and further discussion of MILES and other implicit LES techniques are given in Grinstein *et al.* [31].

The computational set-up is similar to that used by Youngs in Grinstein *et al.* [31] but with increased resolution. The high-density fluid, *ρ*_{1}=1.5, 3, 7 or 20, occupies the region *x*<0 and the low-density fluid, *ρ*_{2}=1, occupies the region *x*>0. Gravity acts the *x*-direction and is chosen to give *Ag*=1. Within each initial fluid, hydrostatic equilibrium with an adiabatic variation is used. The densities at the initial interface are chosen to give the required values of *ρ*_{1} and *ρ*_{2} and the initial pressure at the interface is chosen high enough to give a low Mach number flow (*M*<0.2) which then gives a good approximation to the mixing of incompressible fluids. Low-resolution simulations were performed at even lower Mach number to check that the results quoted here are not significantly affected by compressibility. The list of calculations performed is given in table 1.

At *ρ*_{1}/*ρ*_{2}=1.5 uniform zoning, 1000×1000×1000, is used in the domain −0.5<*x*<0.5: 0<*y*, *z*<1. Coarser meshes are added at the *x*-boundaries to reduce inhibition of mixing zone growth, thus giving a total mesh size of 1050×1000×1000. At higher density ratios, the number of meshes used in the region *x*>0 is increased to allow for increased spike distances. Periodic boundary conditions are used in the *y* and *z* directions. Calculation 4 was a repeat of calculation 5 using the highest resolution feasible at the time of performing this series, 1904×1600×1600, to see whether increasing the resolution has significant effect on the self-similar properties.

Random initial amplitude perturbations, *ζ*(*y*,*z*), are imposed at the initial interface as described by Youngs in Grinstein *et al.* [31]
3.4The coefficients *a*_{mn}, *b*_{mn}, *c*_{mn}, *d*_{mn} are set to Gaussian random variables if
The standard deviation of Gaussian distribution is proportional to , where the power spectrum function is defined as in equation (2.3), and the scaling factor *S* is chosen to give the required overall standard deviation, *σ*, of the perturbation *ζ*. For self-similar growth by mode coupling (case A of §2) *P*(*k*)∼*k*, , and *σ*=0.05Δ*x*. For enhanced self-similar mixing (case B), *P*(*k*)∼*k*^{−3}, and , half the box width, and the perturbation standard deviation is . The values chosen for *ε* are given in table 1. For very low values of *ε*, the power level in the range 4Δ*x*–8Δ*x* is not allowed to fall below the value needed to trigger growth by mode coupling.

For two calculations, a different form was used for the power spectrum. Calculation 12 at *ρ*_{1}/*ρ*_{2} = 3 was a repeat of calculation 8 with the power spectrum changed to *P*(*k*)∼*k*^{−2} while keeping *σ* the same. A *k*^{−2} spectrum is more typical of some practical applications, for example Barnes *et al.* [32], and was chosen to see whether the change in spectrum altered the self-similar properties. The change in spectrum gave a slight reduction in *α* and a slightly greater reduction in *α* towards the end of simulation. Enhanced self-similar mixing requires the initial perturbation to be scale-invariant but not necessarily isotropic. The spectrum *P*(*k*,*θ*)∼*k*^{−3}*f*(*θ*), where , , could be used instead of *P*(*k*) for calculating the random amplitudes in equation (3.4). Calculation 16 at *ρ*_{1}/*ρ*_{2}=7 is a repeat of calculation 14 with the same standard deviation *σ* but with *f*(*θ*)={1+(16*θ*/*π*)^{2}}^{−1}, again to see whether the change in spectrum alters the self-similar properties.

## 4. Results

### (a) Quantities plotted

The key quantities needed for engineering model calibration are derived from the simulations. In each cell, it is assumed that the mixture is in pressure and temperature equilibrium. The fluid volume fractions are then calculated from
4.1The ratio of specific heats for fluid *r* is denoted by *γ*_{r} and the specific heats, *c*_{vr}, are chosen to give temperature equilibrium initially at the interface. For mixing of incompressible fluids, volume fractions correspond to scaled densities, that is *f*_{1}=(*ρ*−*ρ*_{2})/(*ρ*_{1}−*ρ*_{2}). For the RT problems considered here *γ* is 5/3 for each fluid. Fluid volume fractions averaged over *y*, *z* planes are denoted by and . Bubble and spike distances, *h*_{b} and *h*_{s}, are defined to be the distances between the initial interface and the points where and , respectively. However, if bubble distances are calculated directly using this definition, results are susceptible to statistical fluctuations. Instead, an integral width, , which is insensitive to fluctuations, is calculated first and the bubble distance is defined as *h**_{b}=*βW*. If varies linearly with *x* within the mixing zone, *β* should be equal to 3. In practice, *β* is calculated as a function of the Atwood number from mean volume fraction profiles.

The growth rate parameter, *α*, as defined in equation (1.1), is calculated as a function of time by two different methods: firstly, that used by Dimonte *et al.* [16],
and secondly, that used by Ristorcelli & Clark [21],
(note that the mixing zone half-width rather than was used in the previous reference).

The molecular mixing parameter, , which lies in the range 0–1 is used to quantify the degree of molecular mixing as a function of *x*. The global molecular mixing parameter, , is also calculated.

Vertical and horizontal components of the turbulence kinetic energy, *k*_{v} and *k*_{h}, are defined as
The total turbulence kinetic energy is then *k*=*k*_{v}+2*k*_{h}.

The total loss of potential energy is *P*=*K*+*D*, where is the total kinetic energy and *D* is the dissipation. For the MILES results shown here *D* is the implicit ‘sub-grid’ dissipation. For self-similar mixing, *Θ* and *D*/*P* should both approach limiting values.

### (b) Contour plots

Figure 1 shows plane sections of the dense fluid volume fraction for case A, that is evolution of mixing from short wavelength random perturbations. In the images, the vertical extent is limited to the range −0.5<*x*<0.55. The bubble distances are similar for the four density ratios, that is the coefficient *α* varies only slightly. As the density ratio increases the spikes penetrate further. However, the effect is not large and the images at the four density ratios are surprisingly similar.

### (c) Integral properties

Bubble distances are calculated from the integral mix width by using the relation *h*_{b}=*βW*, where *β* is chosen so that *h*_{b} corresponds closely to the position where . As shown in §4*d*, for self-similar mixing at a given density ratio, the volume fraction profiles change little as the growth rate parameter, *α*, is forced to increase. Hence, *β* should depend only on the density ratio. A good estimate is given by *β*=3.433−0.387*A*^{2}.

Figure 2 shows the time variation of *α* for the calculations where mixing evolved by mode coupling from initial short wavelength perturbations. Initial values are high owing to a combination of the influence of initial conditions and the inability to resolve molecular mixing. Figure 2*a* includes results for the higher resolution simulation, 4, and compares the two ways of calculating *α* (*α*^{DY} and *α*^{RC}). At the end of the simulations *α*∼0.022–0.029. There is some uncertainty in the precise limiting value as there are insufficient large-scale structures at the end of the simulations to eliminate statistical fluctuations. The results are very similar to those given in Dimonte *et al.* [16] using much lower resolution, 512×256×256 meshes. For growth by mode coupling at *ρ*_{1}/*ρ*_{2}=3, the higher resolution simulations have given increased confidence in the results but have not significantly changed the conclusions. There are theoretical arguments in favour of using *α*^{RC}. However, for the calculations shown here *α*^{DY} appears to settle down to a limiting value faster than *α*^{RC}. Hence, *α*^{DY} is used throughout the rest of the paper and figure 3*b* shows plots of *α*^{DY} for growth by mode coupling for the four density ratios considered. This confirms the experimental result that the growth rate parameter *α* should vary little with density ratio.

Higher values of *α* are obtained if long wavelength perturbations with a 1/*k*^{3} spectrum are added. Figure 3 shows the effect of increasing the parameter *ε* which determines the magnitude of the long wavelength perturbations. For a low value of *ε*=0.0005, the value of *α* increases to about 0.06, typical of experimental results. The long wavelength random perturbation is truncated at a maximum wavelength of half the domain width as there are too few modes beyond this wavelength to give good statistical behaviour. This truncation of the spectrum appears to give some reduction of *α* towards the end of the simulation. When the long wavelength perturbations are added it proves more difficult to establish a period of accurately constant *α*.

Figure 4*a* shows the variation of *Θ*, global molecular mixing parameter, and *D*/*P*, kinetic energy dissipation fraction, with the mixing zone width *W* at density ratio *ρ*_{1}/*ρ*_{2}=3. For the high-resolution simulation 4, these quantities vary little over a wide range *W* = 0.02–0.16 and settle down to limiting values of *Θ*∼0.80 and *D*/*P*∼0.48. The value of *Θ*∼0.80 agrees very well with the recent DNS results of Livescu *et al.* [22]. Similar results are obtained from the lower resolution simulation 5. When the long wavelength perturbations are added, simulation 8, these quantities are reduced to *Θ*∼0.68 and *D*/*P*∼0.3, but are not so precisely constant at late time. As the mixing zone is forced to grow more rapidly, the process becomes less dissipative. The behaviour is further illustrated in figure 4*b* which shows plots of *Θ* and *D*/*P* versus *α*. As *α* is not precisely constant for each simulation, several points are taken from each simulation corresponding to the range *W*=0.06 to *W*=0.125. It is interesting to note that the effect of density ratio is relatively small. The points for simulation 12, which has a 1/*k*^{2} initial perturbation spectrum, lie on the same curve as for the 1/*k*^{3} initial perturbation spectrum. The change of spectrum has reduced the values of *α* slightly, as shown in figure 3, but not the correlation of internal structure with *α*. For calculation 16, it was also found that the use of a highly anisotropic perturbation spectrum had little effect on the integral properties.

Figure 5 shows plots of the spike/bubble ratio *h*_{s}/*h*_{b}. This ratio is difficult to estimate both experimentally and computationally. The curve for the rocket-rig experiments, taken from Youngs [13], gives *h*_{s}/*h*_{b}=(*ρ*_{1}/*ρ*_{2})^{0.231}. This was based on visual estimates of the bubble and spike positions. There was significant uncertainty in the spike positions owing to the diffuse nature of the boundary. The results of Kucherenko *et al.* [23] were based on X-ray measurements of the mean volume fraction profile and the criteria used for determining bubble and spike positions were and , respectively. The same criteria have been used in figure 5 for the TURMOIL simulations. The ratio *h*_{s}/*h*_{b} seems to vary little with *α* but is very susceptible to statistical fluctuations. There is approximate agreement between the TURMOIL results and those of Kucherenko. In one case, *ρ*_{1}/*ρ*_{2}=3, an observed volume fraction profile is available and a comparison will be given in the next subsection. The rocket-rig data show a little more asymmetry than the simulations but this may be due the uncertainties in interpreting the experiments. Dimonte & Schneider [15] reported greater asymmetry based on measuring *α* separately for bubbles and spikes. Their result, at high Atwood number, *α*_{s}/*α*_{b}=(*ρ*_{1}/*ρ*_{2})^{0.33}, does appear to disagree with the current simulations. Hence, confirmation of the results shown here from alternative simulation techniques would be valuable.

### (d) Mixing zone structure

In this subsection, the distributions of mean volume fraction, molecular mixing parameter and turbulence kinetic energy components , *θ*, *k*, *k*_{v}, *k*_{h}) across the mixing layer are presented. Quantities are plotted against *x*/*W* in order to show the self-similar behaviour. Figure 6 shows profiles of and *θ* for case A when mixing evolves via mode coupling, for the four density ratios considered. The variation with density ratio is surprisingly small. At the higher density ratios, spike distances are greater, in terms of *x*/*W*. At low Atwood number, the profile of *θ* is nearly flat, in agreement with experimental results of Ramaprabhu & Andrews [26]. As the density ratio increases, *θ* on the spike side becomes higher than on the bubble side. For two of the density ratios considered, *ρ*_{1}/*ρ*_{2}=3 and 20, volume fraction profiles are shown for calculations with different values of the mixing rate parameter, *α*, in figure 7. At both density ratios, there is little change in the profile as *α* increases. At *ρ*_{1}/*ρ*_{2}=3, a comparison is made with the experimental results of Kucherenko *et al.* [23] and this shows excellent agreement. Unfortunately, the equivalent experimental data are not available at *ρ*_{1}/*ρ*_{2}=20.

Profiles of turbulence kinetic energies are shown in figure 8. Again, it is found that distributions at a given density ratio vary little as *α* increases. At *ρ*_{1}/*ρ*_{2}=3, the distribution for *k* peaks at *x*/*W*∼1. As the density ratio increases, the location of the peak moves further into the spike region and occurs at *x*/*W*∼3 for *ρ*_{1}/*ρ*_{2}=20. There is a high degree of asymmetry with velocity fluctuations being much higher in the vertical (*x*) direction compared with the horizontal directions. At *ρ*_{1}/*ρ*_{2}=3, peak *k*_{v}/*k*_{h}∼3.7 and this degree of asymmetry is consistent with previous simulations, for example Dimonte *et al.* [16], Cabot & Cook [20]. It also agrees with the experimental estimate of a ratio of approximately 2 for vertical to horizontal velocity fluctuations at low Atwood number [26]. At *ρ*_{1}/*ρ*_{2}=20, the asymmetry increases and gives *k*_{v}/*k*_{h}∼6.

Figure 9 shows the variation of the molecular mixing parameter, *θ*, with *x*/*W*. As noted in §4*a*, the global mixing parameter, *Θ*, decreases as the growth rate increases. Hence, in this case, profiles are expected to change as *α* varies. At both density ratios shown, the values of *θ* reduce as *α* increases and the profile becomes slightly more peaked on the spike side. At *ρ*_{1}/*ρ*_{2}=20, the degree of molecular mixing is much higher on the spike side compared with the bubble side and this can be explained by the higher velocity fluctuations on the spike side as shown in figure 8*b*. Measurements of the molecular mixing parameter at low Atwood number are given in Ramaprabhu & Andrews [26] and these results provide very useful validation of the results presented here. The global mixing parameter for experiments with growth rate parameter *α*∼0.07 and Schmidt number approximately 1, is *Θ*∼0.7 and this ties in reasonably well with figure 4. More convincing validation is obtained by using the volume fraction probability density functions (PDFs) given in Wilson & Andrews [25]. Figure 10 shows a comparison with results from simulation 2 for *ρ*_{1}/*ρ*_{2}=1.5 which has *α*=0.068. The experimental results give volume fraction PDFs for positions above and below the centreline, *y*/*h*_{b}=0.14 and *y*/*h*_{b}=−0.10, where *y* is measured vertically upwards and *h*_{b} is measured to the point where . The computational results are for the equivalent plane sections. On the whole, there is excellent agreement but the peaks at *f*_{1}=1 are more pronounced in the simulation. However, the measured peak at *f*_{1}=1 for *y*/*h*_{b}=0.14 is less than the peak at *f*_{1}=0 for *y*/*h*_{b}=−0.10 and this seems physically unreasonable as there should be more pure fluid at the point furthest from the centreline (as in the simulation). It seems probable that the experimental measurements are uncertain if *f*_{1} is close to unity. In that case, the comparison with experiment confirms that the MILES approach gives good results for molecular mixing in high-Reynolds number flows.

## 5. Concluding remarks

The simple case of self-similar mixing described in equation (1.1) is of fundamental importance both for understanding RT mixing from first principles and for calibration of engineering models. For the ideal case, evolution of mixing from short wavelength perturbations, that is, growth by mode coupling, the current series of simulations give *α*=0.025–0.03. At *ρ*_{1}/*ρ*_{2}=3.0, the results are not significantly different from those of Dimonte *et al.* [16]. However, the significantly increased mesh resolution gives much greater confidence in the results. Experimental values of *α* are usually much higher than this and, in many cases, this is most likely to be owing to the presence of a low level of long wavelength initial perturbations. Enhanced self-similar mixing is obtained if long wavelength perturbations are added with a 1/*k*^{3} spectrum. It is unlikely that experimental initial perturbations are precisely of this form but it is suggested that this provides a better model for experimental situations.

Experiments are always finite. The real problem is not mixing at an infinite plane boundary with a perturbation having finite standard deviation (this should asymptote to *α*∼0.03) but mixing in a finite domain of size *L*, with low levels of perturbations with wavelengths up to size *L*. Then, in most cases, it is likely that the influence of initial conditions will persist throughout the duration of the experiment. Engineering models are essential for complex applications for which LES is not currently feasible. This influence of initial conditions has important implications for engineering models of RT mixing which range from simple models, such as the buoyancy-drag model of Dimonte & Schneider [15] or the kL model of Dimonte & Tipton [33], to more advanced models, such as the BHR model of Besnard *et al.* [34], the multi-fluid model of Youngs [18,35] or the 2SFK model of Llor [36]. All of these models are used to represent mixing of miscible fluids at high Reynolds number. In the BHR, multi-fluid and 2SFK models, the degree of molecular mixing is explicitly modelled and the values of *Θ* quoted here may be used for calibration. In all cases, when model coefficients are chosen a particular value of *α* is obtained (or possibly a value of *α* which varies with Atwood number) and this raises the question ‘what value of *α* should the model be chosen to match?’ A suggested answer is given below.

Firstly, it is noted that the experimental series of Dimonte & Schneider [14,15], which covered a wide range of acceleration profiles and density ratios, was well modelled by a simple buoyancy-drag model with constants chosen to give *α*=0.05, i.e. for this series *α*_{effective}=0.05. Dimonte [19] showed that *α* should vary only slowly (logarithmically) with initial conditions and this explains why it is likely that *α* will be approximately constant for a given set of experiments. Hence, it is proposed that model constants should be derived for a range of values of *α* and an appropriate set chosen for a given complex application. This requires choosing *α*_{effective} for the application. A possible first-principles approach to this might be to estimate the experimental initial conditions, perform a three-dimensional simulation for a simplified version of the complex application and use the mixing rate from this three-dimensional simulation to choose *α*_{effective}.

The MILES results presented in this paper show that the results for self-similar RT mixing at a range of density ratios and for various values of *α* fit in to a simple pattern. The simulations provide much of the data needed to calibrate an engineering model and produce coefficient sets for a range of values of *α*. High Reynolds number DNS is now becoming feasible for the relatively simple problem considered here, see for example Livescu *et al.* [22], and could provide very valuable confirmation of the results presented here. However, the use of three-dimensional simulation on more complex problems (e.g. [37]), which may additionally involve Richtmyer–Meshkov and Kelvin–Helmholtz mixing, to determine *α*_{effective} as suggested above for example, is likely to require the MILES approach.

## Acknowledgements

copyright 2013 British Crown owned copyright/AWE.

## Footnotes

One contribution of 13 to a Theme Issue ‘Turbulent mixing and beyond: non-equilibrium processes from atomistic to astrophysical scales II’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.