2.1 Stratification and initial conditions

We have performed calculations of mixing layer growth with a number of forms for the stratification and initial conditions. While the experiments were performed with nearly incompressible aqueous solutions with Atwood numbers between $4\times 10^{-3}$ and $7\times 10^{-2}$ (Davies Wykes & Dalziel Reference Davies Wykes and Dalziel2014), the present calculations are at a moderate Atwood number of $A_{0}=0.091$ , where $A_{0}=(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})/(\unicode[STIX]{x1D70C}_{2}+\unicode[STIX]{x1D70C}_{1})$ (i.e. a density ratio of 1.2:1), and assume that the materials have a compressible perfect gas equation of state with $\unicode[STIX]{x1D6FE}=5/3$ . The maximal Mach number of the flows in our calculations is ${\sim}0.1$ , so the calculations are close to incompressible; however, it should be noted that the calculations are being performed with the aim of exploring the physics of Rayleigh–Taylor mixing between stratified mixing layers, rather than with the intention of achieving a detailed match to the experimental results. Given the explicit hydrodynamics scheme used by TURMOIL, it is far more numerically efficient to perform this study using flows at a Mach numbers ${\sim}0.1$ where the incompressible flow assumption is still expected to be accurate.

The initial stratified structure is assumed to obey the equation of hydrostatic equilibrium

where $p$ is the pressure, $h$ is the vertical coordinate (this will also be referred to $x$ in the context of the three-dimensional numerical calculations), $\unicode[STIX]{x1D70C}$ is the density and $g$ is the gravitational acceleration. This means that for an imposed density stratification, the pressure field is also determined. The equation of state of the material is important in determining the stability of the stratification. If the material is assumed to have a polytropic equation of state with index $\unicode[STIX]{x1D6FE}$ (as for example $5/3$ in the present calculations), then the density stratification is convectively stable if it satisfies the Schwarzschild criterion (Schwarzschild Reference Schwarzschild1906, Reference Schwarzschild1958; Kaniel & Kovetz Reference Kaniel and Kovetz1967)

where $\unicode[STIX]{x1D6FE}$ is the local value of the effective adiabatic coefficient, $c_{v}$ the specific heat and $T$ the temperature, i.e.

where $c_{s}$ is the adiabatic sound speed.

The interpretation of this is that if a parcel of gas is raised adiabatically through the pressure stratification, for stability its density must decrease less than the density decrease in the surrounding static material (assuming that thermal transport between parcels of material is negligible) – otherwise, it will become buoyant and continue to rise. Note that, when specified in terms of the density gradient, this criterion is equivalent to the Ledoux criterion for convective instability (Ledoux Reference Ledoux1947), which takes into account variations in the molecular mass of the material: as typically applied, these criteria are based on the temperature stratification rather than the density stratification, and the canonical form of the Schwarzschild criterion results if the dependence of the equation of state on the mean molecular mass is not taken into account.

This stability criterion given by (2.3) is similar to that for Rayleigh–Taylor stability (Rayleigh Reference Rayleigh1882; Taylor Reference Taylor1950)

i.e.

but is more restrictive on the zone of stability (i.e. it is possible for Rayleigh–Taylor stable flows to be convectively unstable if the fluid is sufficiently compressible). The corresponding correction to the buoyancy source term for compressible turbulence models has been applied on occasion (Andronov et al. Reference Andronov, Bakhrakh, Mokhov, Nikiforov and Pevnitskii1979; Cloutman Reference Cloutman1987), but is not in widespread use.

In this paper, we will assume that the initial stratification to be isothermal. For an arbitrary density profile, the molecular mass must then vary with position. The easiest way of achieving this within TURMOIL is to specify that all the materials are initially mixtures of materials with different molecular masses. This may be achieved by using two materials, but (as in the experiment) an additional advected passive scalar field is then required to track the interfacial mixing. The initial interface is established by a discontinuous change in mean molecular mass. On either side of the interface, the materials will be convectively stable for this isothermal stratification so long as there is not an extreme gradient in molecular mass $\unicode[STIX]{x1D707}_{m}$ ,

We simulate a box of size 2 units in each dimension. In the direction of gravity, we use reflective boundaries at $x=\pm 1$ (so that there is zero fluid velocity through the boundaries, but free slip along them), while in the $y$ and $z$ directions we assume periodic boundary conditions. We take the gravitational acceleration in the test problem to be $g_{x}=g=1$ , and the interface pressure to be $p=300$ , which means that the calculations should be close to incompressible, since $gH/c^{2}\simeq 1/250$ , with $H=2$ the total height of the domain.

We consider four distinct profiles, which are most clearly seen as the initial time profiles plotted in figure 5. These were:

1. (a) a linear density profile on above and below the interface;

2. (b) quadratic profiles above and below the interface, tending from uniform at the external boundary to a steep dependence close to the interface;

3. (c) as (b), except with the slower density variation in the lower part of the domain close to the interface;

4. (d) as (a), except with the density variation close to the interface clipped to $\pm 0.02$ units of the mean density.

These cases can be understood more clearly by referring to figure 5. We also consider additional linear density profile calculations with an additional regular $\unicode[STIX]{x1D706}=0.2$ perturbation across one dimension of the mesh, with different amplitudes, (e) and (f); in case (f), the perturbation amplitude saturated at a certain level (a case without this saturation behaviour was also run, but the results with the saturated profile seemed more interesting, as we will discuss below). Cases (a) and (b) relate to experimental conditions considered by Davies Wykes & Dalziel (Reference Davies Wykes and Dalziel2014), the comparison with which validates our numerical modelling. The other profiles have been chosen to allow us to study the robustness of our conclusions on the energetics and flow structure for a wider class of distributions than were considered in the experimental work, i.e. to address questions on whether the final structure depends on the symmetry, initial Atwood number or level of surface perturbation in the initial conditions.

The density stratification is defined to run from $\unicode[STIX]{x1D70C}=1$ to $\unicode[STIX]{x1D70C}=1.2$ across each layer for all cases apart from the truncated profile. For example, for the linearly stratified case,

and

Together with the perfect gas equation of state

where $n$ is the particle density and $k_{B}$ the Boltzmann constant, and the assumption that the initial stratification is isothermal, this specifies the profile of molecular mass $\unicode[STIX]{x1D707}_{m}$ , and hence the mass fraction $m_{1}$

We passively advect a marker fluid which tags material which originated in the top of the tank, so this can unambiguously identified.

We initiate the instability by perturbing the initial interface with an initial spectrum of modes with power spectrum $P(\unicode[STIX]{x1D705})\propto \unicode[STIX]{x1D705}^{-3}$ which is dominated by long wavelengths, and has been argued by Inogamov (Reference Inogamov1978) to produce self-similar mixing for simple Rayleigh–Taylor flows. We add perturbing modes in a range of wavelengths between $\unicode[STIX]{x1D706}_{min}=4\unicode[STIX]{x0394}x$ , where $\unicode[STIX]{x0394}x$ is the mesh resolution, and $\unicode[STIX]{x1D706}_{max}=0.2$ , i.e.  $1/10$ of the width of the domain, with an overall standard deviation of $10^{-3}\unicode[STIX]{x1D706}_{max}$ .

2.2 Numerical scheme

The calculations in this paper use the three-dimensional Eulerian hydrodynamics code TURMOIL (Youngs Reference Youngs1991). As applied in the present paper, this code implements a MILES formalism for the calculation of the turbulent mixing flows in the high Reynolds number limit, with dissipation due to the numerical and artificial viscosity of the scheme, but no explicit subgrid model or treatment of physical viscosity. The staggered scheme utilized by TURMOIL is also found not to be subject to the high levels of dissipation at low Mach numbers, which can be an issue for centred Godunov-type schemes (Thornber et al. Reference Thornber, Drikakis, Williams and Youngs2008). The results of all these studies suggest that MILES is indeed accurate in its predictions, for integral properties at sufficiently high Reynolds number, and substantially more economical in the compute resource required to achieve these results.

Figure 1.

Numerical shadowgrams for the case with an initially linear stratification at (a) $t=20$ , (b) $t=40$ , (c) $t=60$ , (d) $t=80$ . This plot is of the two-dimensional Laplacian of mean projected density. As the problem is specified in dimensionless units, the colour scale is also dimensionless: if the domain were scaled into dimensional units using a factor $L$ , the values in the colour map would scale as $L^{-2}$ . The colour map values would also scale linearly with material density.

Figure 2.

Experimental shadowgrams for linear stratification at (a) $t=5~\text{s}$ , (b) $t=10~\text{s}$ , (c) $t=20~\text{s}$ , (d) $t=80~\text{s}$ , from Davies Wykes & Dalziel (Reference Davies Wykes and Dalziel2014) with permission. The experimental time scales can be related to the dimensionless scales used here by dividing them by $\unicode[STIX]{x1D70F}=0.302H/(Ag)^{1/2}=0.254~\text{s}$ , where $2H$ is the total height of the tank and $A$ is the Atwood number in the experiment, although we have not attempted to match the level of initial perturbations which would be expected to influence this comparison. Frame (c) of this plot, at an experimental time of $t=20~\text{s}$ , i.e. $t/\unicode[STIX]{x1D70F}\simeq 80$ , corresponds to a dimensionless time close to that shown in figure 1(d).

Figure 3.

Passive scalar 0.1 % and 99.9 % surfaces at (a) $t=20$ , (b)  $t=80$ .

4.1 Probability density function

The probability density function (p.d.f.) of the initial profile is

The sorted profile satisfies

where we use $h_{s}$ for the height corresponding to each density in the sorted profile and $A$ for the transverse area of the domain. Hence

The sorted profile is monotonic, so $h_{s}(\unicode[STIX]{x1D70C})$ can be inverted uniquely to give $\unicode[STIX]{x1D70C}_{s}(h)$ .

If we specialise to initial profiles with a single unstable interface at $h_{i}$ , then within the sorted mixing layer

where $h_{{<}}(\unicode[STIX]{x1D70C})$ and $h_{{>}}(\unicode[STIX]{x1D70C})$ are the height profiles either side of the interface, which are stable and hence invertible, and at the edges of the mixing layer

4.2 Variational model

In this section, we develop a variational formalism for the structure in the relaxed mixing layer inspired by the work of Lynden-Bell (Reference Lynden-Bell1967) and Tremaine, Hénon & Lynden-Bell (Reference Tremaine, Hénon and Lynden-Bell1986) on the equilibrium of stellar dynamical systems. These papers treated the solution of the equilibrium problem as an entropy maximization problem subject to constraints, as is also discussed in a different context by, e.g. Coolen, Kühn & Sollich (Reference Coolen, Kühn and Sollich2005), Chap. 13. Related approaches have been used before to study the vorticity field, (Carnevale, Frisch & Salmon Reference Carnevale, Frisch and Salmon1981; Montgomery Reference Montgomery, Ray Smith and Grandy1985; Robert & Sommeria Reference Robert and Sommeria1991). This has proved most useful for two dimensional turbulence and its analogues, as the non-uniqueness of the representation in terms of point vortices led to fundamental difficulties when applied in three dimensions (Robert & Sommeria Reference Robert and Sommeria1991). In the present case, the density stratification of the flow provides an alternative functional space in which we can apply this type of method.

In the present case, the probability density function is the volume occupied by material at each density, rather than the mass per unit volume in a six-dimensional phase space (with spatial and velocity coordinates). However, in both cases, the phase space variable is conserved in the absence of mixing processes. For the Vlasov equation applicable to stellar dynamics,

with $\boldsymbol{v}$ the velocity, $\boldsymbol{a}$ the acceleration due to body forces, $\unicode[STIX]{x1D735}$ the spatial vector derivative and

this implies that the phase space density, $f_{p}$ , is conserved on a trajectory satisfying

In a similar manner, the mass conservation equation for incompressible flow

corresponds to conservation of density on a trajectory $\dot{\boldsymbol{x}}=\boldsymbol{v}$ from an initial position $\boldsymbol{x}=\boldsymbol{x}_{0}$ at $t=0$ . (For compressible flow in the absence of dissipation or shocks, a similar relation holds for the specific entropy of the fluid.)

The density p.d.f. at any time may in principle be calculated using

This probability density function corresponds to the volume occupied by material at each density,

In the absence of mixing and as a result of incompressibility, $\unicode[STIX]{x1D70C}(\boldsymbol{x},t)=\unicode[STIX]{x1D70C}(\boldsymbol{x}_{0},0)$ , and the Jacobian $\det (\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\boldsymbol{x}_{0})=1$ , so that $f(\unicode[STIX]{x1D70C},t)=f(\unicode[STIX]{x1D70C}_{1},0)$ . However, coarse graining in the evaluation of $f$ – or, in reality, fine-scale diffusive processes once the transport of material has led to different densities becoming finely intermingled – will mean that the p.d.f. can in fact change as a result of mixing (Lynden-Bell Reference Lynden-Bell1967).

Applying the formalism of Lynden-Bell (Reference Lynden-Bell1967) to the current problem, we define a variational integral

The first term is a configurational entropy in the form derived by Shannon (Reference Shannon1948), $f\ln f$ , to within a constant term which plays no role in the eventual solution, where $f$ is given in terms of the density distribution by (4.12). The integral $\int f(\unicode[STIX]{x1D70C})\ln f(\unicode[STIX]{x1D70C})\,\text{d}\unicode[STIX]{x1D70C}$ will increase monotonically with time as a result of mixing processes. In the absence of other constraints, the variational maximum of this function will lead to a uniform density distribution. In statistical physics, a function of this form was proposed by Boltzmann when formulating his H-theorem, dependent on the phase space density $f(\boldsymbol{x},\boldsymbol{v},t)$ (Boltzmann Reference Boltzmann1964). For the nearly incompressible flow treated in the present paper, we have marginalized the probability density function to be a function of the local fluid density, $f(\unicode[STIX]{x1D70C},t)$ .

The Lagrange multipliers $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ are used to impose constant mass

and potential energy

respectively. It is appropriate to require the potential energy to be constant, as we are interested in the structure of an equilibrium profile, and so by the principle of virtual work all relevant variations must release no further potential energy. There seems to be a natural relation between the physical interchange of fluid parcels during turbulent mixing and the infinitesimal interchanges which are modelled by the functional variation.

The configurational entropy term corresponds to a net attraction towards statistically well-represented states within the ensemble of possible configurations. The logarithmic form chosen used here is standard in statistical physics, and does not require a characteristic scale to be defined for the density, if that is assumed to be strictly positive. Any direct relationship to the physical entropy is, however, obscured by the incompressible limit we are assuming, where the sound speed is assumed to be very large compared to any scale of interest.

Note that the integrals in (4.13) may be transformed between integrals over position and integrals over density because the stability of the equilibrium profile to interchange instability implies that the function $h(\unicode[STIX]{x1D70C})$ is monotonic.

Varying $h(\unicode[STIX]{x1D70C})\rightarrow h+\unicode[STIX]{x1D6FF}h$ gives

where the functional extremum of $I$ satisfies $\unicode[STIX]{x1D6FF}I=0$ . Integrating the first component by parts subject to boundary conditions gives

so that either $\unicode[STIX]{x1D6FF}h=0$ and the initial solution is unperturbed, or

This equation can be integrated once to give

and again to give the equilibrium density profile in implicit form as

a simple, universal prediction for the density profile across the mixing zone. Note that without the potential energy constraint, $\unicode[STIX]{x1D707}\rightarrow 0$ , the variational solution is

i.e. the density becomes uniform, as expected.

This leaves four as-yet free parameters $\unicode[STIX]{x1D706}$ , $\unicode[STIX]{x1D707}$ , $\unicode[STIX]{x1D70C}^{\star }$ and $\unicode[STIX]{x1D705}$ , which determine the position of the inflection point of the profile in height and density, the width and maximum slope of the profile. Two of these parameters may be determined by requiring the density profile to match the initial profile at the edges of the mixing region $h_{1}$  and  $h_{2}$

where $\unicode[STIX]{x1D70C}_{0}(h)$ is the initial density profile. While these nonlinear equations may have multiple solutions, many of these may be eliminated by the requirement that at the edge of the mixing layer, the gradient of the mixed profile must be larger than that of the unmixed profile, as entrainment tends to move the density of the material towards the perfectly mixed density.

A third parameter is determined by requiring that the mass within the mixing region is conserved,

Here the integral over the model for the final profile may be performed analytically

These give three equations constraining the form of the mixing layer based on the structure of the initial solution.

For this model, the remaining parameter is the only way in which the internal structure within the final mixing layer can depend on the structure of the initial stratification. If the dissipation of the convectively released kinetic energy were perfect, then $\unicode[STIX]{x1D705}\rightarrow 0$ , and the value of $\unicode[STIX]{x1D70C}_{0}(h_{1})=\unicode[STIX]{x1D70C}_{0}(h_{2})=\unicode[STIX]{x1D70C}^{\star }$ is determined by the mass conservation constraint

This is the perfect mixing profile of Davies Wykes & Dalziel (Reference Davies Wykes and Dalziel2014). Hence the remaining parameter is effectively a measure of how efficiently the mixing process dissipates energy, rather than using it to reorder stably stratified parts of the density profile.

Figure 8.

Model comparison: dotted is initial state, solid is final state, dashed is model fit. (a) Double linear, (b) symmetric quadratic, (c) asymmetric quadratic, (d) truncated double linear, (e) double linear with long-wavelength perturbation, (f) double linear with larger long-wavelength perturbation. The model parameters have been chosen to give a good fit to the mixing layer profile in (a–d). In (e) and (f), the fit has been kept as that shown in (a), to make it easier to compare these panels, however it should be clear that by reducing energy overshoot, the model can be fit to the density profiles from the TURMOIL calculations at a similar accuracy to the other cases.

4.3 Comparison with results

In figure 8, we compare the densities found at the final time step for each of the initial density profiles we consider, shown as a solid curve, with the tanh-profile form predicted by the variational model, shown as a dashed curve. The agreement between these profiles is generally very good, apart from the furthest extrema of the overshoot in the mixing layer. The fit has not been adjusted for the final two profiles, with the enhanced initial perturbations: while it seems clear that an improved tanh-profile fit could be obtained, including the fit tuned to the case shown in figure 8(a) serves as a useful datum for the reduced level of overshoot found in these cases.

The profiles for which calculations have been performed have been chosen to address the obvious questions which arose when the good agreement of the variational formula with the simple two-linear profile was observed: is the form robust to changes in the source of buoyancy and boundary conditions of final mixing layer as a result of density stratification, figure 8(b,d), surface perturbations, figure 8(e,f) or the symmetry of the initial conditions, figure 8(c). In each case, the prediction of the simple formula for the structure of the final mixing layer appears to be robust, subject only to changes in the total extent of the mixing layer which is an undetermined constant in the model (so long as the calculations are run with sufficient numerical resolution to well resolve the final mixing layer). There is a small amount of mixing beyond the extent of layer predicted by the model, which corresponds to a blurring of the sharp change in gradient predicted at the edge of the mixing region. However, this appears to result from the late-time effects of internal waves and anisotropic turbulence in the stratified envelope, once the nearly isotropic turbulence within the mixing layer itself has died down. There is other evidence in the distribution of the scalar variable within the mixing layer, figure 6, and the structure of the edge of the mixing layer observed at mid to late times in experimental and calculational shadowgraphy data, figures 1 and 2, that this sharp edge does reflect a physical reality.