## Abstract

Rayleigh–Taylor (RT) instability occurs when a dense fluid rests on top of a light fluid in a gravitational field. It also occurs in an equivalent situation (in the absence of gravity) when an interface between fluids of different density is accelerated by a pressure gradient, e.g. in inertial confinement fusion implosions. Engineering models (Reynolds-averaged Navier–Stokes models) are needed to represent the effect of mixing in complex applications. However, large eddy simulation (LES) currently makes an essential contribution to understanding the mixing process and calibration or validation of the engineering models. In this paper, three cases are used to illustrate the current role of LES: (i) mixing at a plane boundary, (ii) break-up of a layer of dense fluid due to RT instability, and (iii) mixing in a simple spherical implosion. A monotone integrated LES approach is preferred because of the need to treat discontinuities in the flow, i.e. the initial density discontinuities or shock waves. Of particular interest is the influence of initial conditions and how this needs to be allowed for in engineering modelling. It is argued that loss of memory of the initial conditions is unlikely to occur in practical applications.

## 1. Introduction

Rayleigh–Taylor (RT) instability is of current concern for inertial confinement fusion (ICF), e.g. Amendt *et al*. (2002). Such instabilities can degrade the performance of ICF capsules, where high-density shells are decelerated by a lower density thermonuclear fluid. In these applications and in many simpler laboratory experiments, the Reynolds number (*Re*) is high. Turbulent mixing will then occur. Direct numerical simulation (DNS) is feasible at a moderate *Re*. However, for most of the experimental situations, calculation of the evolution of turbulent mixing requires some form of large eddy simulation (LES).

The flows of interest here involve shocks and density discontinuities. It is then highly desirable to use monotonic or total variation diminishing numerical methods for both calculating the mean flow and the development of instabilities. As a result, monotone integrated LES (MILES) is strongly favoured when turbulent mixing occurs. Mixing of miscible fluids is considered here, and the *Re* is assumed to be high enough for the effect of the Schmidt number to be unimportant. Wall-bounded flows will not be considered in this paper. The purpose of this paper is to illustrate how MILES can currently be applied to the mixing processes of relevance to ICF. Three simplified problems are considered. Section 4 considers RT mixing at a plane boundary. Section 5 shows the results for the break-up of a dense fluid layer due to RT instability. Finally, in §6, mixing in a spherical implosion is considered. In this case, shock waves are present and mixing due to the related process, Richtmyer–Meshkov (RM) instability, also occurs.

## 2. The role of large eddy simulation

It is not, at present, feasible to apply three-dimensional simulation to turbulent mixing to full-scale applications. Hence, three-dimensional simulation is used to understand the mixing processes from first principles in simplified problems. The results are then used for calibration or validation of the engineering models (Reynolds-averaged Navier–Stokes (RANS) models), which can be applied to real problems. For the time-dependent problems considered in this paper, the ensemble average 〈*ψ*〉 is the fundamental average used for a fluid variable *ψ*. The three test problems considered in §§4–6 are one-dimensional on average and, in these cases, 〈*ψ*〉 depends on *x* (or *r*) only. The average over *y*, *z* (or *θ*, *ϕ*), , usually varies only slightly between calculations using different random numbers to construct the initial perturbation. The values of , obtained from a single three-dimensional simulation, then give an adequate approximation to 〈*ψ*〉, and the distributions of are compared in this paper with the equivalent quantities from the engineering model. The assumption is made that the *Re* is high enough for the amount of mixing to be insensitive to viscosity and diffusivity. It is argued that LES (or MILES) currently gives the best approximation to high *Re* mixing in relatively complex situations such as that considered in §6. On the other hand, DNS for simpler problems at moderate *Re* is very useful for validation of the MILES approach. Wherever possible, experimental results are also used to validate the MILES technique. However, three-dimensional simulation is now being used to obtain results beyond what can be measured in experiments. As these results are used to make improvements to engineering models, it is important to gain confidence by demonstrating approximate mesh convergence for ensemble average behaviour. Ideally, similar results should be obtained for mesh sizes Δ*x*, (1/2)Δ*x* and (1/4)Δ*x*.

As supercomputer power increases, the ‘simplified problems’ to which the MILES technique can be applied will become more and more complex and be used to validate the engineering models in greater detail. Eventually, the MILES approach will be applicable to full-scale applications. However, it is anticipated that the engineering models will be needed for several decades for quick calculations.

## 3. Summary of the numerical method used

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}, and was first used for RT mixing by Youngs (1991). The explicit calculation for each time step is divided into two parts, a Lagrangian phase and an advection or remap phase. The Lagrangian phase solves the equations(3.1)where *ρ*, ** u**,

*p*and

*e*denote density, velocity, pressure and specific internal energy, respectively. The Lagrangian phase is second-order accurate in space and time and is non-dissipative, unless shocks are present. Von Neumann artificial viscous pressure is used to provide the dissipation needed for shocks, rather than a Godunov technique. The remap phase solves the equation(3.2)for all fluid variables,

*ψ*, and uses a third-order van Leer monotonic advection method. Directional splitting is used, thus giving three one-dimensional monotonic steps. The numerical technique is optimized by theoretical analysis of stability and monotonicity and by testing the method on a range of simple test problems, such as growth of RT instability for a single Fourier mode (e.g. Dimonte

*et al*. 2004). It is then applied to turbulent mixing problems

*without change*and this gives a very straightforward approach to simulating turbulence. The monotonicity constraints in the van Leer method, which are a major benefit for modelling flows with discontinuities, provide dissipation at high wavenumbers, hence the method is an example of MILES. It is argued here that this implicit dissipation is similar to that provided by an explicit subgrid-scale model. More details of the Lagrange-remap method, its application to RT and RM mixing and further discussion of MILES and other implicit LES (ILES) techniques are given in Grinstein

*et al*. (2007).

## 4. Self-similar Rayleigh–Taylor mixing at a plane boundary

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 (Read 1984; Dimonte & Schneider 2000) have shown that the dominant length scale increases as mixing evolves. If mixing is self-similar, then dimensional reasoning suggests that 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(4.1)where *α* varies slightly with the Atwood number, *A*, and is typically in the range 0.05–0.07. Dimonte *et al*. (2004) give a summary of the experimental results.

Self-similar mixing can be achieved in two different ways. First, case (A), if the initial perturbation consists of a random combination of short-wavelength perturbations, then larger and larger scales will eventually evolve from nonlinear interactions between small 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 observed at late time may evolve directly from initial perturbations at these scales. If the initial random perturbation is such that ‘amplitude (proportional to) wavelength’ at all scales up to the size of the experiment, then as pointed out by Inogamov (1999), self-similar mixing may occur at an enhanced rate. In more precise terms, the power spectrum, *P*(*k*), for the initial perturbation height and the standard deviation of the surface, *σ*, are given by(4.2)where *k* is the wavenumber and *P*(*k*)∼1/*k*^{3}. After the experiments of Read (1984) and Dimonte & Schneider (2000) 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 (Youngs 1994). This conclusion was supported by Dimonte *et al*. (2004), in which results for seven different MILES techniques (including TURMOIL), using a mesh size of 512×256×256, all indicated similar low values, *α*=0.025±0.003. The results suggested that the observed growth rates were all influenced by initial conditions. By use of a simple model, Dimonte (2004) 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.

In order to illustrate this behaviour, results are given here for higher resolution TURMOIL simulations for cases (A) and (B). The number of meshes used is 720×600×600, with periodic boundary conditions in the *y*- and *z*-directions. The domain size was 1.45×1×1. For the major part of the domain, the mesh size was the same in all three directions. However, coarser meshes were added at the *x*-boundaries to reduce inhibition of mixing zone growth. At the interface, the fluid densities were *ρ*_{1}=3 and *ρ*_{2}=1. Gravity was in the *x*-direction and chosen to give *Ag*=1. Within each initial fluid, hydrostatic equilibrium with an adiabatic variation was used. The initial pressure at the interface was chosen high enough to give a low Mach number flow (*M*<0.2), and this gives a good approximation to the mixing of incompressible fluids. For case (A), random initial amplitude perturbations were used with wavelengths in the range 4Δ*x* to 8Δ*x*, and for case (B), the perturbation spectrum (4.2) was used with wavelengths up to half the domain width and *σ*=0.00025 (a very low value).

The iso-surfaces (*f*_{1}=0.99) in figure 1 show the increase in scale length as time proceeds (*f*_{1} and *f*_{2} are the volume fractions for fluids 1 and 2, respectively). For a mixture in pressure and temperature equilibrium, the fluid volume fractions are calculated from(4.3)The 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, i.e. .

In order to demonstrate the *gt*^{2} growth of the mixing zone width, an integral measure of the width is used (which is insensitive to statistical fluctuations), , where are values of the volume fractions averaged over the *y* and *z* planes. If varies linearly with *x* within the mixing zone, the bubble penetration should be *h*_{b}=3*W*. The edges of the mixing zone are slightly diffuse and the approximation *h*_{b}=3.3*W* will be used here. Figure 2*a* shows a plot of *h*_{b} against *Agt*^{2} for cases (A) and (B). For both cases, three calculations are shown with different choices for the random numbers used to generate the initial perturbations. For case (A), *α*∼0.027, in agreement with previous results, and for case (B) with *σ*=0.00025, *α*∼0.056, typical of observed values. These results strongly suggest that a low level of long-wavelength initial perturbations is needed to account for the experimental growth rates.

Very-high-resolution DNS results, 3072^{3} meshes, were given by Cabot & Cook (2006), corresponding to case (A). In this paper, an improved way of calculating *α* as a function of time was used, , which showed that *α* was high early in the simulation and dropped to a low value after a ‘mixing transition’, as defined by Dimotakis (2000), occurred. The value of *α* then slowly increased to ∼0.025 by the end of the simulation when reached 3×10^{4}. Similar plots of *α* versus time from the MILES are given in figure 2*b*. For case (A), the trend is similar to that in the DNS. There is some doubt about the precise late-time, high *Re* limiting value of *α*, and higher resolution simulations are needed to investigate this. However, it can be concluded that relatively low resolution MILES gives similar late-time results to DNS at *Re*=3×10^{4}. The low values of *α* for case (A) from both sets of results show that initial conditions are likely to have a large effect in practical applications, and this needs to be taken into account both in comparison with experiment and in calibrating engineering models. This is an important conclusion and it is worth noting that it came from three-dimensional simulations *not* from experiments. For case (B), *α* reduces a little at the end of the simulation. This is probably due to the cut-off in the initial perturbation spectrum at low wavenumbers. Higher mesh resolution is needed to obtain a period of more nearly constant *α* for this case.

There is a strong similarity between the results given here and the behaviour of mixing in turbulent shear flows according to George (in press), who has argued for many years that the influence of initial and upstream conditions persists throughout the entire flow.

Much more detailed information can be extracted from the MILES results, and some results are given in figure 3, which shows plots of and the molecular mixing parameter, , versus scaled distance, (*x*−*x*_{0})/*W*, and probability density functions (PDFs) for *f*_{1} at *x*=*x*_{0}, where *x*_{0} is the initial interface location. Initial conditions influence the values of *Θ*; the slightly higher mean value for case (B), *Θ*∼0.7, is close to the measured values of Ramaprabhu & Andrews (2004). Values of the PDF for *f*_{1} at three different times show that the internal structure changes little when the resolution of the mixing zone increases by a factor of four.

## 5. Break-up of a dense layer due to Rayleigh–Taylor instability

In this section, a more complex RT mixing problem is considered: the break-up of a dense fluid layer of depth (1/2)*H*. The initial configuration, with *x* measured downwards in the direction of gravity, hasandThe size of the computational domain is 2*H*×*H*×*H* and three mesh resolutions are used: 256×128×128; 512×256×256; and 1024×512×512. Initial random long-wavelength perturbations are added at the lower unstable boundary using the power spectrum given by (4.2), with *σ*=0.00025*H*. This gives an RT growth parameter of *α*≈0.056. Two different density ratios have been considered, *ρ*_{1}*/ρ*_{2}=1.5 and 10, and the results are presented as a function of a scaled time parameter (Boussinesq scaling). In terms of this scaled time, the evolution of the flow at the lower density ratio is similar to that of the low Atwood number experiments of Jacobs & Dalziel (2005), which provide useful validation of the results shown here. For the higher density ratio, the break-up of the layer is shown in a sequence of two-dimensional sections (figure 4).

Figure 5 shows the effect of changing the density ratio. In terms of mean volume fractions, (rather than mass fractions), and the scaled time *τ*, this effect is surprising small. There are some differences: at *τ*=2.5, for the higher density ratio, there is asymmetry in the mean volume fraction profile and the molecular mixing parameter. At both density ratios, a high degree of molecular mixing occurs, and this prevents the dense fluid from simply accumulating at the bottom of the tank at late time. Experiments of this type are very difficult to perform at high density ratios, and the results given here illustrate how MILES may be used to further the understanding of the mixing process.

Figure 6 shows the effect of mesh resolution for the case *ρ*_{1}/*ρ*_{2}=10. The evolution of the mean volume fraction profiles is very similar for the three mesh resolutions and this gives confidence in the results. As previously noted, the current purpose of the LESs is the validation of an engineering model, and this is illustrated here. The model, described in Youngs (1994, 1995), is based on the equations of multi-phase flow, with turbulent diffusion and Reynolds stress terms added. Mass exchange between the phases is used to represent the decay of concentration fluctuations (i.e. model the time evolution of the molecular mixing parameter). The results obtained from the model, a one-dimensional calculation using 80 meshes, have been added to figure 6. As the initial conditions, in reality, do influence the mixing, they must be allowed for in the model. For the case considered here, this is achieved by choosing model coefficients to match the properties of self-similar RT mixing, case (B) with *α*=0.056. When this is done, the evolution of the distributions of the three fluids is quite well matched. Again, there is a similarity to turbulent shear flow according to George (in press), who argued that one-point closure model constants should depend on initial conditions. This section gives an illustration of how MILES can provide validation of the engineering model, which cannot be achieved by comparison with currently available experimental data.

## 6. Mixing in a simple spherical implosion

The final problem considered here is closer to the ICF applications and was previously considered by Youngs & Williams (2008). The geometry is greatly simplified. However, the initial random interface perturbations are considered to have a realistic power spectrum. An important issue is how the initial perturbations influence the mixing at the end of the implosion. As for the previous problem, comparison is made with the one-dimensional engineering model. The initial geometry for the unperturbed spherical implosion, using dimensionless units, is

Perfect gas equations of state are used with *γ*=5/3. The ratio of the specific heats for the two fluids is 20 : 1, thus giving equal initial temperatures in the two regions. The implosion is driven by applying a prescribed pressure, *p*, versus time at the outer boundary, initially at *r*=12, which moves in a Lagrangian manner. The applied pressure is obtained by linear interpolation in time between the valuesThe resulting radius–time plot obtained from a one-dimensional Lagrangian calculation is shown in figure 7. RM instability occurs when the incident shock wave passes through the dense fluid/light fluid interface at *t*∼0.5. This amplifies the initial perturbation. Turbulent mixing occurs towards the end of the implosion, when the interface decelerates, due to a combination of RT and RM instabilities. The perturbations that initiate the late-stage mixing are now a combination of three factors: (i) the initial surface finish, (ii) amplification of the perturbation by the leading shock wave, and (iii) spherical convergence. LES is essential for answering the question: does the engineering model adequately represent the influence of initial conditions? Another feature of this test case is that the time variation of the turbulence (expressed in terms of Lagrangian time derivatives) is rapid compared with many other RANS model applications. This further strengthens the need for the use of LES for model validation.

For this problem, TURMOIL uses spherical polar coordinates and three-dimensional simulation is carried out for a sector of the sphere centred at the equator,This reduces the computational resources needed and avoids the mesh singularity at *θ*=0. A further extension of the numerical method is also used: Lagrangian zoning is used in the radial direction. Moreover, regions of one-dimensional Lagrangian zoning are used near the origin, 0<*r*_{0}<2.5, and near the outer boundary, 11.3<*r*_{0}<12, where *r*_{0} denotes the initial radius. This overcomes the problem associated with the mesh singularity at the origin and also limits the three-dimensional calculation to the region near the interface where mixing occurs. The number of zones used in the *r*, *θ* and *ϕ* directions is as follows: coarse mesh, 220×120×120; standard mesh, 440×240×240; and fine mesh, 880×480×480. The fine mesh simulation ran for 40 hours using 360 processors of the Atomic Weapons Establishment's Cray XT3. Random amplitude perturbations were initially applied to the light fluid/dense fluid interface. The power spectrum used was given by*C* is chosen so that the s.d. of the pertubation is σ=0.0005, where Figure 8 shows sections from the fine mesh simulation. The shading represents the volume fraction of the denser fluid. Figure 9*a* shows plots of the turbulent mixing zone limits versus time for the three simulations. These limits are defined as the values of *r* where the dense fluid volume fraction averaged over *θ* and *ϕ*, , is equal to 0.01 and 0.99. The width of the mixing zone reduces slightly as the mesh is refined. However, the effect of mesh size is not very large, and this suggests that the fine mesh calculation should give a good indication of the amount of mixing. Figure 9*b*,*c* shows distributions of the dense fluid volume fraction and the molecular mixing parameter at the end of the simulations. Again, the influence of mesh size is small, and this gives considerable confidence in the results. Figure 9 also includes a comparison with the one-dimensional engineering model previously referred to in §5. If coefficients are chosen to match the properties of self-similar RT mixing with a higher value of *α*=0.07, the model represents the mixing process well. The need to use a higher value of *α* reflects a greater influence of initial conditions compared with the problem considered in §5. Without guidance from LES results, it would be very difficult to calibrate the engineering model.

## 7. Concluding remarks

It is argued here that a relatively simple MILES technique is able to give accurate results for high *Re* mixing in complex problems with shocks and initial density discontinuities. DNS is very useful for understanding *Re* effects in simple problems and in justifying the MILES approach. On the other hand, MILES can be applied to problems of much greater complexity. At present, many full-scale practical applications are not feasible for MILES. However, the MILES technique currently makes an essential contribution to the validation of the engineering models, which can be used for such problems and greatly adds to what can be achieved by comparison with experimental data. A key issue discussed in this paper is the influence of initial conditions. For the situations considered here, it is argued that initial conditions are likely to influence the amount of mixing throughout the entire duration of the problem. LES is essential for understanding the influence of initial conditions in relatively complex flows and for validating the treatment of initial conditions in the engineering model.

## Acknowledgments

© Crown Copyright 2009, the MoD, UK.

## Footnotes

One contribution of 16 to a Discussion Meeting Issue ‘Applied large eddy simulation’.

- © 2009 The Royal Society