Rayleigh–Taylor mixing: direct numerical simulation and implicit large eddy simulation

Rayleigh–Taylor (RT) instability plays an important role in many areas of research. It can degrade the performance of inertial confinement fusion (ICF) capsules, see for example Amendt et al (2002). It also occurs in many astrophysical flows, see for example Fryxell et al (1991). The simplest case consists of fluid with density ρ₁ resting initially above fluid with density ρ₂ < ρ₁ in a gravitational field g. Experiments using incompressible fluids with low viscosity, low surface tension and random initial perturbations 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². The experiments indicate that the depth to which the mixing zone penetrates the denser fluid, generally known as the bubble distance, is given by:

$$\begin{eqnarray}&&{h}_{{\rm{b}}}=\alpha \displaystyle \frac{{\rho }_{1}-{\rho }_{2}}{{\rho }_{1}+{\rho }_{2}}g{t}^{2}=\alpha Ag{t}^{2},\end{eqnarray} \tag{ 1 }$$

where α 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}_{{\rm{s}}}.$ The ratio ${h}_{{\rm{s}}}/{h}_{{\rm{b}}}$ is a slowly increasing function of the density ratio ${\rho }_{1}/{\rho }_{2}.$ Three-dimensional numerical simulation of this self-similar turbulent mixing problem, and some more complex extensions of the simple case, is the subject of the present paper.

Some experimental data is available. However, there are few detailed measurements at high-Atwood number (A > ½). Much of the understanding of the mixing process has been obtained from 3D numerical simulations. Direct numerical simulation (DNS) is feasible at moderate Reynolds number and is essential for assessing the influence of Reynolds number (Re) and Schmidt number (Sc). However, many of the experiments and applications are in the high-Re regime. It is argued here that high-resolution large eddy simulation is the most computationally efficient way of calculating the high-Re behavior and that this approximation is needed for the more complex situations in which RT mixing occurs. As the problems of interest have initial discontinuities,  implicit large eddy simulation (ILES) is strongly favored. However, a few researchers have used large eddy simulation (LES) with explicit models for the dissipative processes. Mixing of miscible fluids is considered here and for ILES or LES to be valid the Reynolds number is assumed to be high enough for the flow to be beyond the 'mixing transition' where, according to the hypothesis of Dimotakis (2000), the Schmidt number does not affect the amount of molecular mixing.

Section 2 reviews previous research on 3D numerical simulation of RT mixing, DNS, ILES and LES. The emphasis is on the simple incompressible situation (1). The choice of numerical technique is a matter of considerable controversy. However, very few papers have made direct comparisons between the various approaches. The main purpose of the current paper is to give some direct comparisons of DNS and ILES results using the same computational framework. The emphasis is on quantities which are of most importance in practical applications, such as ICF: the extent to which to which the fluids mix and the degree to which mixing occurs at a molecular level. Section 3 summarizes the numerical method used here for both DNS and ILES. Results for four test cases: (i) single-mode RT instability, (ii) self-similar RT mixing, (iii) three-layer mixing and (iv) a tilted-rig RT experiment are given section 4. Concluding remarks are given in section 5.

The present paper focuses on 3D numerical simulation of RT mixing, which has been feasible since the 1990s. Other aspects of research into RT instability have been reviewed by various authors. For example, Sharp (1984) gives a review of the early research into Rayleigh–Taylor instability. Kull (1991), Abarzhi (2010), Abarzhi and Rosner (2010) give reviews of theoretical modeling approaches. Low Atwood number experiments, for which detailed diagnostics are feasible, are reviewed by Andrews and Dalziel (2010). Wide-ranging reviews of many aspects of the instability processes are given in the theme issue of the Philosophical Transactions of the Royal Society edited by Sreenivasan et al (2013).

Much of the content shown here was presented at the 'Turbulent Mixing and Beyond' workshop held in 2014. The paper addresses key themes and topics of the TMB program, in particular non-equilibrium turbulent processes, Rayleigh–Taylor instabilities and interfacial mixing.

A gt² scaling law for the evolution of RT instability from random perturbations was proposed in a number of early publications, Birkhoff (1954), Sharp and Wheeler (1961), Belen'kii and Fradkin (1965) and Anuchina et al (1978). However, the form of the scaling law, given by equation (1) was established by the experiments of Read (1984) and Youngs (1989) which indicated α ∼ 0.06. Later experiments of Dimonte and Schneider (1996, 2004), gave α ∼ 0.05 and also considered cases where g varied with time. However, when high-resolution 3D simulation became feasible, calculations with 'ideal initial conditions' (small random short wavelength random perturbations) indicated very low values of α ∼ 0.03, Youngs (1994), Linden et al (1994). In the latter paper, a simple model was used to show that if the value of α for ideal initial conditions is so low, then mixing would be significantly enhanced by very low levels of long wavelength initial perturbations. The low values of α for ideal initial conditions were supported by Dimonte et al (2004) in which results for seven different ILES techniques (including the TURMOIL method described in section 3) were given. The mesh size used was 512 × 256 × 256 and all simulations 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. These results confirmed that the observed growth rates are likely to be 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.

Cook and Dimotakis (2001), Young et al (2001) and Cook and Zhou (2002) performed the first DNS for RT mixing and gave results for the early stages of the mixing process and the detailed structure of the mixing zone. Cook et al (2004) used LES with 1152³ meshes and an explicit dissipation model. Both of these approaches showed the development of a 'mixing transition', as defined by Dimotakis (2000), when an inertial range begins to form. The LES gave α ∼ 0.027, very similar to the ILES results of Dimonte et al (2004). Ristorcelli and Clark (2004) used DNS to understand the time variation of α in the early stages of the mixing process. Subsequently very high resolution DNS results, using 3072³ meshes, were given by Cabot and Cook (2006). This paper showed that α was high early in the simulation and the dropped to a low value after the mixing transition. The value of α then slowly increased to ∼0.025 by the end of the simulation, when $Re={h}_{{\rm{b}}}{\dot{h}}_{{\rm{b}}}/\nu ,$ reached 3 × 10⁴, close to the values obtained from the previous ILES/LES. The Atwood number used in these papers did not exceed A = 0.5. A form of LES was also used by Burton (2011) to perform simulations at high Atwood number, A = 0.5–0.96. In all cases, a low value of α ∼ 0.02 was obtained. The spike to bubble ratio, ${h}_{{\rm{s}}}/{h}_{{\rm{b}}},$ was estimated and found to be significantly less at high Atwood number than reported for the immiscible-fluid experiments of Dimonte and Schneider (2000). The review paper, Statsenko et al (2013), shows results for the total mixing zone width at density ratios 3:1 to 40:1. The results at 3:1, assuming ${h}_{{\rm{s}}}\sim {h}_{{\rm{b}}},$ indicate α ∼ 0.03. High resolution (up to 4032 × 4096²) DNS results for RT mixing at a range of Atwood numbers, A = 0.04–0.9 have given by Livescu et al (2010, 2011) and a review of DNS approaches is given by Livescu (2013). All the calculations referred to so far (DNS, ILES, LES) have indicated low values of α ∼ 0.02–0.03 for mixing of miscible fluids arising from ideal initial conditions.

It appears that α only has a unique value, corresponding to loss of memory of the initial conditions, if small random short wavelength perturbations are assumed. The best estimate is then α ∼ 0.025 with little dependence on Atwood number. However, some doubt still remains. Cabot and Cook (2006) have speculated that there might be a difference in the behavior at ultra-high Reynolds number, beyond present computational capabilities.

Several researchers have tried to improve agreement between simulations and experiments by taking the initial conditions into account. In Dalziel et al (1999) measurements of the perturbations generated when the barrier separating the two initial fluids was withdrawn were used to give a more realistic ILES model of the experiments. Ramaprabhu and Andrews (2004) used measurements of initial perturbations for channel flow RT experiments to initialize an ILES model. This enabled the observed value of α ∼ 0.07 to be matched. Mueschke and Schilling (2009a, 2009b) used a DNS model for the same experiments to obtain detailed mixing statistics and further analysis of the DNS data was given in Schilling and Mueschke (2010). Ramapraphu et al (2005) and Banerjee and Andrews (2009) used an incompressible ILES technique to show that RT mixing has a strong dependence on initial conditions. In particular, it was found that if long wavelength initial perturbations were used with a k⁻³ power spectrum (k = wavenumber), values of α up to 0.08 could be achieved. Many of the experiments have used immiscible fluids whereas the simulations assume miscible mixing. It is possible that measured values of α have been influenced by surface tension. It could be argued that low values of surface tension suppress fine-scale mixing and thereby increase α by increasing the effective Atwood number. Glimm et al (2001) advocated the use of front-tracking simulations in such cases and expressed concern about the treatment of molecular mixing in the ILES approach. The paper by Glimm et al (2013) uses front tracking and LES to model a range of experiments, many of which have used immiscible liquids. Initial conditions were estimated, from the early-time experimental photographs in some cases, and good agreement was obtained with the observed growth rates which had values of α in the range 0.05–0.08. For turbulent mixing experiments with random initial perturbations, initial conditions are never known precisely. However, it is clear that several authors have achieved improved agreement with experiment by estimating the initial conditions as accurately as possible.

Since Youngs (1994), TURMOIL ILES have been performed with increasing mesh resolution. In Youngs (2003) the mesh resolution was increased to 720 × 480 × 480. For random short wavelength initial perturbations the value of α was 0.027, in agreement with the earlier simulations. Inogamov (1978) pointed out that self-similar mixing at an enhanced rate could be obtained if random long wavelength perturbations with amplitude ∝ wavelength, which corresponds to a k⁻³ power spectrum. It was found that if such a perturbation spectrum was used with a very small amplitude standard deviation of only 0.00025× the computational domain width, then α was increased to about 0.06, a typical experimental value. Similar results using increased resolution were given in Youngs (2009). In Youngs (2013) ILES results were shown for mesh resolutions ∼2000 × 1000 × 1000. The increase in mesh resolution had little effect on the low value of α obtained when short wave length initial perturbations are used. A range of calculations was carried at different density ratios and with different levels of the k⁻³—spectrum perturbations, giving a range of values of α ∼0.025–0.1. Properties of the mixing zone were derived as functions of α.

In this section emphasis has been given to the values of α calculated by various researchers and to the agreement with experiment. An important conclusion is that there is a large difference between the mixing rates for simulations using ideal conditions and those observed experimentally. The papers referred to here have also given many results for the detailed structure of the mixing zone. There has been some experimental validation of this. However, as regards the detailed structure, simulations have provided much more information than can be obtained experimentally.

Several papers have considered simple extensions to the self-similar mixing case (1). ILES results for a three-fluid mixing problem were given Youngs (2009) and further results for this case are shown here in section 4. In Andrews et al (2014) the initial interface was tilted relative to the direction of gravity; DNS, incompressible and compressible ILES results were given. In section 5 here further comparisons of DNS and ILES are shown for this case. RT mixing in a tall tube was considered by Lawrie and Dalziel (2011) and ILES was compared with a simple model which predicted ${h}_{{\rm{b}}}\sim {t}^{2/5}.$ Williams (2017) has used ILES (TURMOIL) for RT mixing with non-uniform stratification in each fluid. Variable acceleration, with g alternating in sign was investigated by incompressible ILES, Dimonte et al (2007) and Ramaprabhu et al (2013). Alternating sign acceleration and evolution of RT mixing with an additional localized perturbation was considered in Statsenko et al (2013). In the paper by Olson et al (2011) initial shear is added to investigate the combined effect of Rayleigh–Taylor and Kelvin–Helmholtz instabilities. This paper also compares results for DNS and LES with an explicit dissipation model. It is shown that the LES results converge to the DNS result at high resolution.

The choice of the most suitable numerical methods for 3D simulation of Rayleigh–Taylor and Richtmyer–Meshkov (RM) mixing remains an area of significant uncertainty. High-order or spectral methods are often favored for DNS, in particular for the accurate calculation of fine-scale properties, and have been used in previous the RT DNS of Cabot and Cook (2006) and Livescu (2013). Use of a simple approach is considered here, a 2nd/3rd order Lagrange-remap method (TURMOIL). It is argued that this method is suitable for relatively complex problems with initial density discontinuities and shocks. Rehagen et al (2017) gives a comparison of a similar Lagrange-remap technique with the high-order method of Cabot and Cook (2006), for RT mixing at moderate mesh resolution. Agreement between the two numerical techniques was found for global quantities such as mixing zone widths and the degree of molecular mixing. In these simulations the dissipative scales were well resolved at the start of the calculation (near DNS) but at the end the simulations were effectively LES. A more advanced LES approach, the stretched-vortex subgrid-scale model, has been applied to RM mixing, Lombardini et al (2011), but has not been used for the RT mixing problems considered here. A wide range of numerical techniques is available. However, the most appropriate choices for simulating the variety of RT and RM mixing cases is very complex problem and additional studies are required to resolve this challenging issue.

The results shown here have been obtained by using the TURMOIL Lagrange-remap hydrocode which calculates the mixing of compressible fluids. This was initial used for ILES of RT mixing, Youngs (1991) and solved the Euler equations plus advection equations for fluid mass fractions. More details of the Lagrange-remap method, its application to RT and RM mixing and further discussion of ILES techniques are given by Grinstein et al (2007). The remap phase uses the monotonic advection method of Van Leer (1977) to prevent spurious oscillations. In turbulent flows, the monotonicity constraints provide the required dissipation at high-wavenumbers and this is the basis of the ILES approach. Perfect gas equations of state are used for each fluid species with constant values of the specific heat ratio, ${\gamma }_{r},$ and specific heat, ${c}_{{\rm{v}}r}.$ The mixture is assumed to be in pressure and temperature equilibrium. The pressure of the mixture is then given by

$$\begin{eqnarray*}&&p=(\gamma -1)\rho e\,{\rm{with}}\,\gamma =\displaystyle \sum _{r}\,{\gamma }_{r}{c}_{{\rm{v}}r}{m}_{r}/\displaystyle \sum _{r}{c}_{{\rm{v}}r}{m}_{r}.\end{eqnarray*}$$

For DNS the Navier–Stokes equations are used and include the effects of viscosity (ν = kinematic viscosity), species diffusion and heat conduction. Simplified material properties are used. The same diffusivity, D, is used for all fluid species and for heat (i.e. the conductivity coefficient is $\rho {c}_{{\rm{p}}}D$). This leads to a simplified form of the equations, previously used by Kokkinakis et al (2015), for density, ρ, velocity, ${u}_{i},$ mass fractions, ${m}_{r}$ and specific internal energy, e:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{\partial (\rho {u}_{j})}{\partial {x}_{j}}=0\\ &&\displaystyle \frac{\partial (\rho {u}_{i})}{\partial t}+\displaystyle \frac{\partial (\rho {u}_{i}{u}_{j})}{\partial {x}_{j}}=\,-\displaystyle \frac{\partial p}{\partial {x}_{i}}-\displaystyle \frac{\partial {S}_{ij}}{\partial {x}_{j}}+\rho {g}_{i}\\ &&{\rm{with}}\,{S}_{ij}=-2\rho \nu {e}_{ij}^{\prime },{\rm{where}}\,{e}_{ij}^{\prime }\,{\rm{is}}\,{\rm{the}}\,{\rm{deviatoric}}\,{\rm{strain}}\,{\rm{rate}}\\ &&\displaystyle \frac{\partial (\rho {m}_{r})}{\partial t}+\displaystyle \frac{\partial (\rho {m}_{r}{u}_{j})}{\partial {x}_{j}}\,\,=\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\rho D\displaystyle \frac{\partial {m}_{r}}{\partial {x}_{j}}\right),\,\,\,\\ &&\displaystyle \frac{\partial (\rho e)}{\partial t}\,+\displaystyle \frac{\partial (\rho e{u}_{j})}{\partial {x}_{j}}\,=-p\displaystyle \frac{\partial {u}_{j}}{\partial {x}_{j}}+\displaystyle \frac{\partial }{\partial {x}_{i}}\left(\rho D\displaystyle \frac{\partial h}{\partial {x}_{j}}\right),\\ &&{\rm{where}}\,h=\gamma e\,=\,{\rm{enthalpy}}.\end{eqnarray} \tag{ 2 }$$

For ILES the Euler equations are used; ν and D are set to zero. For the calculations shown here, the adiabatic constants are ${\gamma }_{1}={\gamma }_{2}=5/3$ and the specific heats are chosen to give temperature equilibrium initially at the unstable interface: ${\rho }_{1}{c}_{{\rm{v}}1}={\rho }_{2}{c}_{{\rm{v}}2}$ (specific heats are only needed for the calculation of volume fractions—see below). Heat conduction does not have a significant effect. For DNS the Schmidt number is unity, i.e. ν = D, as in many previous RT DNS. The initial interface pressure is set to a sufficiently high value to give a peak Mach number less than 0.2 and this gives near incompressible flow for which mean volume fraction distributions are unaltered by increase in pressure. Hence the results shown are representative of incompressible gaseous mixing.

In the numerical method, the Lagrange and remap phases are unchanged for DNS. In particular, monotonic advection is retained in the remap phase. Extra steps in the calculation are included for viscosity and diffusion of mass fractions and enthalpy. For accurate DNS, the viscous and diffusive scales need to be well enough resolved for dissipation in the remap phase to be negligible.

A convenient way of quantifying the extent to which species r mixes with the other fluids is by use of the mixed mass, Zhou et al (2016)

$$\begin{eqnarray}&&{ {\mathcal M} }_{r}=\displaystyle \int \rho {m}_{r}(1-{m}_{r}){\rm{d}}V={M}_{r}-\displaystyle \int \rho {m}_{r}^{2}{\rm{d}}V.\end{eqnarray} \tag{ 3 }$$

M_r denotes the total mass of fluid r which remains constant. Hence the production of ${ {\mathcal M} }_{r}$ corresponds to the dissipation of ${m}_{r}^{2}.$

The amount of mixing which occurs in the diffusion steps (the resolved mixing) can be calculated from

$$\begin{eqnarray*}&&{ {\mathcal M} }_{{\rm{r}}}^{{\rm{res}}}=\displaystyle \iint 2\rho D\displaystyle \frac{\partial {m}_{r}}{\partial {x}_{j}}\displaystyle \frac{\partial {m}_{r}}{\partial {x}_{j}}{\rm{d}}V{\rm{d}}t.\end{eqnarray*}$$

For accurate DNS, ${ {\mathcal M} }_{{\rm{r}}}\approx { {\mathcal M} }_{{\rm{r}}}^{{\rm{r}}{\rm{e}}{\rm{s}}}.$ However, as some mixing will occur in the remap steps, ${ {\mathcal M} }_{{\rm{r}}}\gt { {\mathcal M} }_{{\rm{r}}}^{{\rm{r}}{\rm{e}}{\rm{s}}}.$ A measure of the accuracy of the DNS is given by the ratio

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{\rm{r}}}=\displaystyle \frac{{ {\mathcal M} }_{{\rm{r}}}-{ {\mathcal M} }_{{\rm{r}}}^{{\rm{r}}{\rm{e}}{\rm{s}}}}{{ {\mathcal M} }_{{\rm{r}}}}.\end{eqnarray} \tag{ 4 }$$

An alternative approach would be to use a non-dissipative remap phase. The resolution would need to be sufficient to make any deleterious effects due to spurious oscillations negligible and it is thought that this would be difficult to quantify. The advantages of the approach chosen here is that a robust computational method is retained and a method of checking the accuracy of the DNS is available. A similar check on viscous dissipation is also made; for accurate DNS, the loss of kinetic energy in the remap phase should be small compared to the dissipation directly due to viscosity.

For quantifying the amount of mixing that takes place, it is useful to construct fluid volume fractions

$$\begin{eqnarray*}&&{f}_{r}=\displaystyle \frac{({\gamma }_{r}-1){c}_{{\rm{v}}r}{m}_{r}}{\displaystyle {\sum }_{s}({\gamma }_{s}-1){c}_{{\rm{v}}s}{m}_{s}}.\end{eqnarray*}$$

Many previous researchers e.g. Youngs (1994), Andrews et al (2014), have used fluid volume fractions to define local (θ_r) and global (Θ_r) molecular mixing parameters

$$\begin{eqnarray}&&{\theta }_{r}=\displaystyle \frac{\overline{{f}_{r}(1-{f}_{r})}}{{\bar{f}}_{r}(1-{\bar{f}}_{r})},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{r}=\displaystyle \frac{\displaystyle \int {f}_{r}(1-{f}_{r}){\rm{d}}V}{\displaystyle \int {\bar{f}}_{r}(1-{\bar{f}}_{r}){\rm{d}}V}.\end{eqnarray} \tag{ 6 }$$

In general, the overbar denotes the ensemble average. For the cases shown here the ensemble average is approximated by using plane averaging (sections 4.2 and 4.3) or line averaging (section 4.4). If only two fluids are present, ${\theta }_{1}={\theta }_{2}=\theta \,{\rm{and}}\,{{\rm{\Theta }}}_{1}={{\rm{\Theta }}}_{2}={\rm{\Theta }}.$ The values of the molecular mixing parameters vary from 0 (no molecular mixing) to 1 (fully mixed at a molecular level).

As explained by Youngs (2007) the spatial variation of the dissipation in TURMOIL ILES can be accurately quantified. The Lagrange phase is non-dissipative in the absence of shocks as in the test cases considered here. The remap phase is split into separate steps for x, y and z advection. Consider the x-remap of a quantity ϕ per unit mass. The equations solved are

$$\begin{eqnarray*}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{\partial (\rho u)}{\partial x}=0,\\ &&\displaystyle \frac{\partial (\rho \phi )}{\partial t}+\displaystyle \frac{\partial (\rho \phi u)}{\partial x}=0.\end{eqnarray*}$$

These equations imply

$$\begin{eqnarray*}&&\displaystyle \frac{\partial (\rho {\phi }^{2})}{\partial t}+\displaystyle \frac{\partial (\rho {\phi }^{2}u)}{\partial x}=0.\end{eqnarray*}$$

Hence in the absence of numerical dissipation ${\textstyle \int \rho {\phi }^{2}{\rm{d}}V}$ should be conserved in the remap step. However, if monotonic advection is used for ϕ, dissipation of this quantity will occur. The spatial distribution of the dissipation can be derived as follows. Suppose the difference equations for conservation of mass and ϕ are written as

$$\begin{eqnarray*}&&{M}_{j}^{1}={M}_{j}^{0}+\delta {M}_{j-\tfrac{1}{2}}-\delta {M}_{j+\tfrac{1}{2}},\\ &&{\phi }_{j}^{1}{M}_{j}^{1}={\phi }_{j}^{0}{M}_{j}^{0}+{\phi }_{j-\tfrac{1}{2}}^{VL}\delta {M}_{j-\tfrac{1}{2}}-{\phi }_{j+\tfrac{1}{2}}^{VL}\delta {M}_{j+\tfrac{1}{2},}\end{eqnarray*}$$

where ${M}_{j}^{0}$ and ${M}_{j}^{1}$ are the cell masses at the start and end of the step, $\delta {M}_{j+\tfrac{1}{2}}$ is the mass flux accross the cell boundary $x={x}_{j+\tfrac{1}{2}}\,{\rm{and}}\,{\phi }_{j+\tfrac{1}{2}}^{VL}$ is a van Leer limited value of $\phi .$

These difference equations are used to calculate the reduction in ${\textstyle \int \rho {\phi }^{2}{\rm{d}}V.}$ The summation is re-arranged to give

$$\begin{eqnarray}&&\displaystyle \sum _{j}\{{M}_{j}^{0}{({\phi }_{j}^{0})}^{2}-{M}_{j}^{1}{({\phi }_{j}^{1})}^{2}\}=\displaystyle \sum _{j}\,{{\rm{\Phi }}}_{j+\tfrac{1}{2}}\delta {M}_{j+\tfrac{1}{2}}.\end{eqnarray} \tag{ 7 }$$

The contribution to the dissipation, ${{\rm{\Phi }}}_{j+\tfrac{1}{2}}\delta {M}_{j+\tfrac{1}{2}},$ is then divided between meshes $j\,{\rm{and}}\,j+1.$ If ϕ is set to each of the velocity components in turn, this technique gives the spatial distribution of the dissipation of kinetic energy in the remap phase. For the ILES this gives the 'sub-grid-scale' energy dissipation. Further details are given in Youngs (2007). This idea was first used by DeBar (1974) to calculate kinetic energy loss, but not in the context of turbulence simulation. If ϕ is set to m_r, the dissipation calculated corresponds to the production of ${ {\mathcal M} }_{r}$ (3).

The TURMOIL method used second-order accurate numerics in the Lagrange phase and third-order monotonic advection in the remap phase. Second-order, central differencing is used for the viscous and diffusion terms. As noted in section 2, high order methods have, as a rule, been used in previous DNS studies of RT mixing. Comparisons between TURMOIL and higher order methods, typically of order 5, for a simple viscid vortical flow are given in Shanmuganathan et al (2014) and Tsoutsanis et al (2015). TURMOIL compares very favorably with the higher order methods. It does require a finer mesh (by a factor ∼1.5) to obtain converged results but is computationally much faster. Hence it is considered a feasible approach for DNS.

4.1. Single mode Rayleigh–Taylor instability

The first application is single mode Rayleigh–Taylor instability at density ratio ${\rho }_{1}/{\rho }_{2}=20$ and with Ag = 1. The width of the domain, set to unity here, is twice the perturbation wavelength, λ, and the kinematic viscosity (the same for each fluid) is chosen so that the most unstable wavelength (if diffusivity is neglected) is

$$\begin{eqnarray}&&{\lambda }_{{m}}=4\pi {\left(\displaystyle \frac{{v}^{2}}{Ag}\right)}^{\mathrm{1/3}}=\tfrac{1}{2}\lambda .\end{eqnarray} \tag{ 8 }$$

Diffusivity, with $D=\nu ,$ also reduces the growth rate of short wavelength perturbations. However, the parameter ${\lambda }_{m}$ is nevertheless considered to be a useful estimate of the viscous/diffusive scale. The Mach number of the simulations is sufficiently low to give near-incompressible flow. There is initially a sharp interface but this quickly becomes diffuse. Numerical results are shown in figure 1.

Zoom In Zoom Out Reset image size

Volume fraction contours $({f}_{1}=0.1,\,0.3,\,0.7,\,0.9)$ for three different mesh sizes are compared at t = 2.5. Figure 1(a) shows results at an earlier time with viscosity and diffusivity omitted. This clearly shows that viscosity and diffusivity have a large effect. The effect of viscosity and diffusivity is accurately calculated with ${\lambda }_{m}=32\times \,{\rm{the}}\,\,\text{mesh}\,\text{size},\,{\rm{\Delta }}x,$ figure 1(d). However, a good approximation is obtained with ${\lambda }_{m}=8{\rm{\Delta }}x,$ figure 1(b). In that case Δ, as defined by (4) is 0.07 and this implies that some of the mixing occurs in the remap phase. It is concluded that the influence of diffusivity on the amount of mixing is adequately calculated with less than fully resolved DNS.

4.2. Self-similar mixing at a plane boundary

One of the ILES calculations described in Youngs (2013) is repeated with viscosity and diffusivity included. For the case chosen, high-density fluid, ${\rho }_{1}=20,$ occupies the region x < 0 and the low density fluid, ${\rho }_{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. Short wavelength random perturbations are used with wavelengths in the range 4Δx to 8Δx. The resolution used for the DNS is same as that used for the ILES, 1550 × 1000 × 1000 meshes. The width of the domain in the y and z directions is unity. For viscosity and diffusivity $\nu =D={\rm{a}}\,{\rm{constant}}$ and the value of the constant is chosen so that the viscous scale, (8), is ${\lambda }_{m}=17\,{\rm{meshes}}.$ Then, according to the results of the previous sub-section, the early time effect of viscosity and diffusivity should be well resolved. However, the dissipation scales also need to be resolved at later times when turbulence develops. For homogeneous turbulence dissipation occurs at the Kolmogorov scale, ${\eta }_{{\rm{K}}},$ which is related to ε, the kinetic energy dissipation rate per unit mass by

$$\begin{eqnarray*}&&{\eta }_{{\rm{K}}}={\left(\displaystyle \frac{{\nu }^{3}}{\varepsilon }\right)}^{\mathrm{1/4}}.\end{eqnarray*}$$

For RT mixing, experiments and simulations, Dalziel et al (1999), Cabot and Cook (2006) show high wavenumber spectra similar to k^−5/3. Hence a similar scaling law for the dissipation scales (at least at Sc ∼ 1) is expected. For self-similar RT mixing, ε is proportional to potential energy loss rate per unit mass, $\sim Ag{\dot{h}}_{{\rm{b}}}$ and this implies that the dissipation scales should decrease slowly with time, ${\eta }_{{\rm{K}}}\sim {h}_{{\rm{b}}}^{\mathrm{1/8}}.$ Hence the need to resolve the turbulent dissipation scales will eventually determine the mesh resolution required.

The ILES and DNS give very similar values of α ∼ 0.030. Viscosity and diffusivity reduce the growth rate of the short wavelength initial perturbations used here. This leads to a small time delay in the DNS. Hence DNS results at t = 3.2 are compared with ILES results at t = 3.0 when the mixing zone widths are similar. Plane sections for volume fraction distributions are shown in figure 2. It is apparent that much more fine scale structure is present in the ILES. It should be noted that different random numbers were used in each simulation. Hence the large scale structures do not appear in the same locations.

Zoom In Zoom Out Reset image size

The main purpose of this section is to compare results for the amount of molecular mixing, as defined by the parameters θ (5) and Θ (6). Figures 3 and 4 compare the ILES and DNS results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows values of ${\bar{f}}_{1}\,{\rm{and}}\,\theta$ plotted against the scaled distance variable $x/W$ where $W=\int \,{\bar{f}}_{1}(1-{\bar{f}}_{1}){\rm{d}}x$ is an integral measure of the mixing zone width. It is evident that DNS and ILES give very similar results. The two ILES results shows that there are small differences in results due to the choice of random numbers for the initial perturbations.

Figure 4 shows that at the end of the simulations, the global mixing parameter reaches very similar values, ${\rm{\Theta }}\sim \mathrm{0.75-0.8},$ for DNS and ILES. However, the early time behavior is different, At the start of the DNS, Θ is close to unity when the initially sharp interface diffuses away and then drops to a lower value as the instability evolves. Finally, as small scale structure develops, Θ rises to a limiting value of about 0.8. This is, of course the physically correct behavior for the chosen values of ν and D. Figure 4 also shows a plot of the quantity Δ, as defined by equation (4). This never exceeds 5% and this indicates that the DNS is adequately resolved. For the ILES there is initially a sharp dip in Θ as the instability evolves without the development of fine scale structure. Then, when sufficient fine scale structure (inertial range) can be resolved, Θ rises to a limiting value of 0.75–0.8. In this case, the early time behavior is determined by the numerics; the ILES results for the degree of molecular mixing are only meaningful when the late-time limiting value is approached. This occurs when W ∼ 0.02. The width of the mixing zone, ∼8 W according to figure 3, is then equal to 160 times the mesh size. For the DNS higher resolution is needed to reach the limiting value of Θ. The dynamic range is not very high in the DNS shown here but should be sufficient to estimate the limiting value of Θ. Αs a limiting value has been reached, it is considered appropriate to describe the flow as 'turbulent' by the end of the simulation. The dynamic range is not high enough to calculate the high-Reynolds behavior of some of the fine-scale properties.

The molecular mixing parameter, θ, varies from θ = 0.6 on the bubble side to θ = 0.95 on the spike side for both DNS and ILES; this is in good agreement with the DNS results of Livescu et al (2011). The spike/bubble ratio ${h}_{{\rm{s}}}/{h}_{{\rm{b}}}\sim 1.7$ is also the same for both DNS and ILES and in agreement with previous results, for example Burton (2011).

An appropriate estimate of the Reynolds number, based on the bubble distance and velocity, is

$$\begin{eqnarray*}&&\,\begin{array}{l}{Re}\,=\,\displaystyle \frac{{h}_{{\rm{b}}}{\dot{h}}_{{\rm{b}}}}{\nu }\sim \displaystyle \frac{2{h}_{{\rm{b}}}\sqrt{\alpha Ag{h}_{{\rm{b}}}}}{\nu }\\ \,\,\,\text{}=2500\,{\rm{at}}\,{\rm{the}}\,{\rm{end}}\,{\rm{of}}\,{\rm{the}}\,{\rm{simulation}}.\end{array}\end{eqnarray*}$$

The Reynolds number achieved is close to that in the RT experiments described by Mueschke et al (2009) in which Θ reaches a plateau level for gaseous mixing at Sc ∼ 0.7. For water channel experiments, Sc ∼ 1000, a significantly higher Reynolds number, probably of order 10⁴ (not achieved in the experiments) is needed to reach the plateau level. This high Sc case is beyond current DNS capability.

The comparison of the DNS and ILES results leads to the following conclusions. The two approaches give very similar results for the late-time degree of molecular mixing, when a high Reynolds number is achieved. This should dispel concerns about the ability of ILES to calculate molecular mixing. If the high Reynolds number behavior of the global properties of the mixing zone is the main point of interest then the ILES approach is able to estimate this with significantly less computer resources than DNS. If the influence of finite Reynolds number is required then DNS is, of course, essential. Moreover, DNS is required for many small scale properties of the mixing zone.

4.3. Three-fluid mixing

In this sub-section a more complex RT mixing problem is considered: the break-up of a dense fluid layer. In some ICF capsule designs, for example Amendt et al (2002), multiple shells are used and this has provided the motivation for a three-layer test case. This is significantly more difficult than the previous test case. For self-similar mixing at a plane boundary the main requirement is to reach the self-similar state before the end of the simulation. For the three-fluid case, turbulent flow needs to established before the intermediate layer breaks up, and for DNS this has proved difficult. It would be desirable to show DNS for significantly higher Re than has been possible at present. Nevertheless, some interesting results for the influence of viscosity on the break-up of the dense layer are given.

DNS results are given for one of the cases for which ILES results were given in Youngs (2009). The initial configuration, with x measured downwards in the direction of gravity has

$$\begin{eqnarray*}&&\text{Fluid}\,\text{1}:\,0\lt x\lt \tfrac{1}{2}H\,{\rm{density}}\,{\rho }_{1}\\ &&\text{Fluid}\,\text{2}:\,\tfrac{1}{2}H\lt x\lt H\,{\rm{density}}\,{\rho }_{2}\\ &&\text{Fluid}\,\text{3}:\,H\lt x\lt 2H\,{\rm{density}}\,{\rho }_{3}.\end{eqnarray*}$$

A range of density ratios was considered in Youngs (2009). For the comparisons shown here ${\rho }_{2}/{\rho }_{1}={\rho }_{2}/{\rho }_{3}=3$ is chosen. Gravity is chosen to give Ag = 1 and H = 1 here. Results are given as a function of a scaled time $\tau \,=\sqrt{Ag/H\,}t.$ Initial random long wavelength perturbations are added at the lower unstable boundary with a power spectrum given by

$$\begin{eqnarray*}P(k)=\left\{\begin{array}{ll}C/{k}^{3} & 4{\rm{\Delta }}x\lt \lambda =2\pi /k\lt \,\tfrac{1}{2}H\\ 0 & {\rm{otherwise}}.\end{array}\right.\end{eqnarray*}$$

The value of C was chosen so that the standard deviation of the initial perturbation was $\sigma =0.00025H.$ For ILES this gave approximate self-similar growth, before dense layer break-up, with an enhanced value of α ∼ 0.06, a typical experimental value.

The mesh used for ILES was 1024 × 512². According to Youngs (2009), this gave near-converged results. For DNS, three mesh sizes were used: 512 × 256², 1024 × 512², 2028 × 1024². In each case the values of $\nu =D$ were chosen to give ${\lambda }_{m}=16{\rm{\Delta }}x$ so that the resolution of the viscous and diffusive scales at early time was similar to that obtained in the last sub-section. A convenient estimate of the Reynolds number for the flow is

$$\begin{eqnarray*}&&{Re}\,=\,\displaystyle \frac{{h}_{{\rm{b}}}\dot{{\mathop{h}\limits^{}}_{{\rm{b}}}}}{\nu }\,{\rm{when}}\,{{h}}_{{\rm{b}}}={\rm{\text{width}}}\,{\rm{\text{of}}}\,{\rm{\text{dense}}}\,{\rm{\text{layer}}},\,\tfrac{1}{2}{H}\\ &&=\,\displaystyle \frac{2hbH\sqrt{\tfrac{1}{2}\alpha AgH}}{\nu }.\end{eqnarray*}$$

This gives Re ∼ 500, 1400, 3400 for the three DNS mesh sizes. For the lowest value of Re ∼ 500 the calculation has been repeated using the intermediate mesh size, thus giving ${\lambda }_{m}=32{\rm{\Delta }}x,$ to show the effect of mesh size at given Re.

Calculations are run to τ = 10, when the amount of mixing is close to the final value. Figures 5 and 6 show volume fraction distributions for the ILES and the highest resolution DNS. It is evident that the overall time scale for the mixing process is very similar for the two calculations. However, more fine-scale structure is present in the ILES, in spite of the lower resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows the evolution of the integral mix width for the dense fluid layer, ${W}_{2}={\textstyle \int {\bar{f}}_{2}(1-{\bar{f}}_{2}){\rm{d}}x.}$ For the DNS at the two higher values of Re, 1400 and 3400, the mixing rate is similar to the ILES. However, for the lowest value of Re = 500 there is a significant difference in behavior. Initially, τ < 0.5, there is rapid mixing due to diffusion at the initially sharp interface. Then as the instability evolves, τ ∼ 2, the mixing rate is significantly higher than for the lowest values of viscosity. This is attributed to two factors. For this case mixing is, to a large extent, controlled by the long wavelength initial perturbations which are not significantly affected by viscosity. Fine scale structure is suppressed by viscosity and this leads to less molecular mixing (see figure 10). The increased density contrast then tends to increase the mixing rate. Figure 7 also confirms that the DNS results for Re = 500 are insensitive to mesh size.

Zoom In Zoom Out Reset image size

Distributions of volume fractions averaged over (y, z)-planes are shown in figures 8 and 9. At τ = 4, the distributions for the highest Re DNS and the ILES are very similar. This is also true for the final state at τ = 10. It is interesting to note that for the final state, figure 9, the effect of Reynolds number is surprisingly small.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 shows the evolution of Θ₂, the global mixing parameter for the dense fluid, as defined by (6). The DNS results show that the Reynolds number has a significant effect. As for the two-fluid mixing case (section 4.2) the interface initially diffuses away and gives values of Θ₂ close to unity. Then as mixing evolves Θ₂ drops to a value of ∼0.55 and finally rises as fine scale structure develops. The minimum value changes little with Reynolds number but time at which the minimum occurs increases as viscosity increases. For high Reynolds number self-similar mixing, Θ₂ should reach a plateau level, as implied by figure 4. For the range of Reynolds numbers considered here this is not achieved before mixing reaches the top of the dense layer at τ ∼ 3. Hence DNS at considerable higher Reynolds number is needed to model the high-Reynolds number behavior. For the ILES a plateau level of Θ₂ ∼ 0.68 is reached for τ = 1–3 and this is typical of experimental results for Sc ∼ 1, Mueschke et al (2009), which indicate Θ ∼ 0.70.

Zoom In Zoom Out Reset image size

For the DNS calculations shown in figure 10 values of Δ₂, as defined by (4), at τ  = 10 are as follows. For the calculations at Re = 500, ${{\rm{\Delta }}}_{2}=0{\rm{.04}}\,{\rm{with}}\,{\lambda }_{m}=16{\rm{\Delta }}x$ and ${{\rm{\Delta }}}_{2}=0{\rm{.007}}\,{\rm{with}}\,{\lambda }_{m}=32{\rm{\Delta }}x.$ Hence the calculation at the higher resolution is well resolved but it is evident from figures 7 and 10, that the calculation with ${\lambda }_{m}=16{\rm{\Delta }}x$ gives very similar results. For the higher Reynolds number simulations with ${\lambda }_{m}=16{\rm{\Delta }}x,$ ${{\rm{\Delta }}}_{2}=0\text{.06}\,\text{for}\,{\text{Re}}\,=\,1400$ and ${{\rm{\Delta }}}_{2}=0\text{.065}\,\text{for}\,{\text{Re}}\,=\,3400.$ The resolution of the DNS is not quite so good but should be adequate for the results given. For the DNS shown here the resolution of the initial viscous scale is kept the same and this implies ${\rm{\Delta }}x\sim {\nu }^{\mathrm{2/3}}.$ The variation of Δ₂ with Re suggests that in order to maintain the accuracy of the DNS at late time, the resolution of the Kolmogorov scale should be kept the same and this implies ${\rm{\Delta }}x\sim {\nu }^{\mathrm{3/4}}.$ In order to model the high-Reynolds number behavior for this test case with DNS calculations at much higher Re are desirable, a factor 10 higher perhaps. If the computational time step is determined by the Courant number, ${\rm{\Delta }}x\sim {\nu }^{\mathrm{3/4}}$ scaling then implies an increase in computer time of a factor 1000, which is impractical at present.

The conclusions drawn from the DNS and ILES results for the three-fluid test case are as follows. The overall progress of the mixing process, as indicated for the volume fraction distributions, is very similar for the ILES and the highest resolution DNS. The ILES appears to give an adequate approximation to the high-Re behavior. DNS is able to assess the influence of viscosity at moderate Reynolds number and at Schmidt number Sc ∼ 1.

4.4. A tilted-rig Rayleigh–Taylor experiment

The final problem considered is the tilted-rig test case described by Andrews et al (2014). This is based on an experiment described by Youngs (1989) for which a tank containing hexane (top, ρ₂ = 0.66 g cm⁻³) and NaI solution (bottom, ρ₁ = 1.89 g cm⁻³) was accelerated downwards. The experimental rig was inclined at an angle β = 5°46' to the vertical. This tilted the initial interface relative to the sides of the tank and led to a large scale overturning motion in addition to the growth of a turbulent mixing zone. This flow is, on average, two-dimensional and the main purpose of the test case is to provide a benchmark problem for 2D Reynolds-averaged Navier–Stokes (RANS) models and was used for this purpose by Denissen et al (2014). It is a challenging test case for DNS; the viscosity used both here and by Andrews et al (2014) needs to be increased by a factor 10 in order to resolve the dissipation scales.

In the simulations, figure 11(a), the tank is static and a gravitational force is applied in the vertical direction. The domain size used in the simulations is ${L}_{x}={L}_{y}\,=15\,{\rm{cm}},\,{L}_{z}=24\,{\rm{cm}}.$ Reflective boundary conditions are used in the x and z directions. Periodic boundary conditions are used in the y direction. The value of L_y is increased from the experimental value of 2.5 cm in order to reduce fluctuations for y-averaged quantities.

Zoom In Zoom Out Reset image size

In the experiment the acceleration takes about 20 ms to rise to a near-constant plateau level. The incompressible ILES and DNS reported in Andrews et al (2014) used the measured acceleration history for g_z(t). However, for the compressible TURMOIL simulations this was impractical and a constant value of g_z = 0.035 cm ms⁻² was used. A correction was made to the time scale for comparisons between the simulations. Results are given here for a scaled time of

$$\begin{eqnarray*}&&\tau =\sqrt{\displaystyle \frac{A{g}_{z}}{{L}_{x}}}t=2.117\end{eqnarray*}$$

which corresponds approximately to figure 11(d).

The DNS of Livescu, shown in Andrews et al (2014) began with a slightly diffuse interface. Moreover, as the simulations began with g_z = 0, the interface further diffused away before the instability evolved. This led to a significant delay in the development of the instability. The TURMOIL DNS results with constant g_z and a sharp initial interface reduce this effect and it is argued that this gives a more direct comparison of ILES with DNS.

The mesh sizes used for the TURMOIL simulations shown here are 600 × 600 × 960 (ILES) and 1600 × 1600 × 2560 (DNS). The value used for ν = D was, 1.57 × 10⁻⁴, the same as the lowest value used by Livescu. Then for the DNS ${\lambda }_{m}=15{\rm{\Delta }}x$ and the resolution of the early stages is similar to that in the previous sub-sections. As described in Andrews et al (2014), random initial perturbations have a k⁻² power spectrum with wavelength range 0.2–7.5 cm and standard deviation 0.0075 cm. Calculations with the same initial perturbations but without the tilt give α ∼ 0.05, typical of the experimental values.

Two estimates of the Reynolds number for the DNS are given. Andrews et al (2014) used a value based on the maximum perturbation wavelength, ${\lambda }_{\max }=7.5\,{\rm{cm}}:$

$$\begin{eqnarray*}&&{Re}\,=\,{\lambda }_{{\rm{\max }}}\sqrt{\displaystyle \frac{Ag{\lambda }_{{\rm{\max }}}}{A+1}}\phantom{\rule{}{2.25ex}}/\nu =14000.\end{eqnarray*}$$

The Reynolds number based on the bubble distance and velocity (measured in the center of tank) at the time selected, t = 2.117, is approximately

$$\begin{eqnarray*}&&{Re}\,=\,\displaystyle \frac{{h}_{{\rm{b}}}{\dot{h}}_{{\rm{b}}}}{\nu }\approx 2400\end{eqnarray*}$$

and this is similar to the value achieved for DNS shown in section 4.2.

For this flow ensemble averages are functions of x and z and are approximated by averaging in the y-direction. Figure 12 shows comparisons of volume fractions and turbulence kinetic energy. As in Andrews et al (2014) the latter is define by

$$\begin{eqnarray*}&&k=\displaystyle \frac{\overline{\rho \{{({u}_{x}-{\tilde{u}}_{x})}^{2}+{u}_{y}^{2}+{({u}_{z}-{\tilde{u}}_{z})}^{2}\}}}{\bar{\rho }}\\ &&\,\,{\rm{where}}\,{\tilde{u}}_{x}=\displaystyle \frac{\overline{\rho {u}_{x}}}{\bar{\rho }},{\tilde{u}}_{z}=\displaystyle \frac{\overline{\rho {u}_{z}}}{\bar{\rho }}.\end{eqnarray*}$$

The ILES and DNS give very similar distributions. The width of the mixing zone is slightly greater in the DNS. As for the three-layer problem, this is attributed to reduced fine-scale mixing in the early stages of the DNS which gives greater density contrast and increased buoyancy.

Zoom In Zoom Out Reset image size

Figure 13 shows comparisons of the kinetic energy dissipation rate per unit mass, ε, and the molecular mixing parameter, θ. Again it is found that ILES and DNS give similar results. For the DNS, $\varepsilon =\overline{{S}_{ij}{e}_{ij}}/\bar{\rho },$ and for the ILES, ε is the implicit numerical dissipation calculated as described in section 3. It is of particular interest to see how this compares with the viscous dissipation in the DNS. At the time shown, the quantity Δ, as defined by (3), is 0.07. Hence the DNS is reasonably well resolved.

Zoom In Zoom Out Reset image size

The evolution of the global mixing parameter, Θ, is shown in figure 14. The trend is similar to that shown in figures 4 and 10. For the DNS, Θ reaches values close to unity when the initial interface diffuses away, then dips to ∼0.5 and finally rises to ∼0.7. From τ = 1.75 onwards, when the effect of finite Reynolds number is presumed small, values of Θ for ILES and DNS are near-identical. Figure 14 also indicates the expected behavior at very high Reynolds number; a high plateau level of Θ ∼ 0.7 should quickly be reached and Θ should increase slightly towards the end of the experiment. The ILES gives an approximation to this behavior at an earlier stage that the DNS. The previous DNS of Livescu showed similar behavior. Θ dipped to a value of ∼0.47 and the rose to values approaching the ILES results. However, because of the time delay (explained earlier), the previous DNS did not reach the stage where the ILES and DNS results agreed.

Zoom In Zoom Out Reset image size

The total computer time (processor–hours) for the DNS was a factor 50 greater than that for the ILES. Hence ILES provides the simplest way of estimating the high-Reynolds number behavior. As noted earlier, if the effect of finite Reynolds number is required, DNS is essential.

Three dimensional simulation, DNS, ILES and some examples of LES with explicit dissipation models, has made a major contribution to understanding the process of turbulent mixing by Rayleigh–Taylor instability. The detailed structure of the mixing zone is difficult to obtain from experiments, especially at high Atwood number, and 3D simulation is able to fill this gap. However, one important result is learnt from experiments. Typical experiments give α ∼ 0.06 whereas simulations, using a wide range of techniques, with ideal initial conditions (random short wavelength perturbations) give α ∼ 0.026. A likely explanation of this is the presence of low levels of long wavelength initial perturbations in the experiments. The experimental results show that mixing in real situations will usually be significantly greater than that for the ideal theoretical case.

The main purpose of the paper has been to give a direct comparison of results for DNS and ILES using the same numerical framework. It has been shown here that, provided the mesh resolution is sufficient to resolve fine-scale structure within the mixing zone, ILES and DNS give very similar results. This applies to quantities such as distributions mean volume fractions, molecular mixing parameter or turbulence kinetic energy. The Navier–Stokes results given here are not fully-converged DNS. However, the deviation from fully-converged DNS has been quantified and it argued that the mainly global results given here are accurate. For quantities such as enstrophy, which are determined by the small scales present in the flow, DNS is required to get meaningful results (probably with higher resolution or with higher order numerics being required to obtain accurate results). DNS may also be needed if more complex physics is present (combustion for example). If the main purpose of the simulations is to obtain the high-Reynolds behavior of global properties, then this can be achieved using ILES with much less computer resources than with DNS. The early stages of the mixing process are influenced by viscosity and diffusivity and DNS is essential for understanding the influence of Reynolds number and Schmidt number on the mixing process.

The DNS results given here have used a relatively simple numerical approach. It is argued that satisfactory results have been obtained for the influence of finite Reynolds number on the mainly global results presented. Further research is needed to assess the suitability of the method for the accurate calculation of fine-scale properties.

For many practical applications 3D simulation is not feasible at present and some form of engineering modeling, for example RANS, is needed. Both DNS and ILES have an essential contribution to make to calibration and validation of such models. At present, DNS results may be used for model calibration/validation in simple situations, such as one-dimensional mixing layers. However, as engineering models need to be applied to much more complex problems, validation is needed for a wide range of benchmarks. It is advocated that comparisons of RANS and ILES have an essential role here.

The computational results shown in this paper were generated when the author was employed by AWE, Aldermaston, United Kingdom. The paper contains material © British Crown Owned Copyright 2017/AWE.