A Fokker–Planck approach to a moment closure for mixing in variable-density turbulence

ABSTRACT

We develop a theory for the cascade mixing terms in a moment closure approach to binary active scalar mixing in variable-density turbulence. To address the variable-density complications we apply, as a principle and constraint, the conservation of the probability density function (PDF) through a Fokker–Planck equation with bounded sample space whose attractor is the beta PDF with skewness. Mixing is related to a single-point PDF as a realisability principle to provide mathematically rigorous expressions for the small scale statistics in terms of largescale moments. The problem of the unknown small-scale mixing is replaced with the determination of the drift and diffusion terms of a Fokker–Planck equation in a beta-PDF-convergent stochastic process. We find that realisability of a beta-convergent process requires the mixing time-scale ratio, taken as a constant in passive scalar mixing, to be a function of the mean mass fraction, mean fluid density, the Atwood number, the density-volume correlation and moments of the density field. We develop and compare the new model with direct numerical simulations data of non-stationary homogeneous variable-density turbulence.

KEYWORDS:

• Favre averages
• mean turbulent reaction rates
• variable-density turbulence

1. Introduction

In conventional second moment closure approaches the mixing process is modelled as a drag with a time-scale set from which we derive the same second moment equation in which the mixing terms can be identified as a model for the mixing in a Navier–Stokes-based moment closure approach.

We specifically address the modelling in a moment-based approach to a material mixing cascade in variable-density turbulence (VDT) in which the density field and the hydrodynamics are tightly coupled.

Unlike passive-scalar mixing, mixing between fluids of very different densities exhibits a high degree of skewness even when there are equal amounts of the two pure fluids (in the symmetric initial condition case), see Livescu and Ristorcelli [1,2]. The probability density function (PDF) at the stages of a mixing process by VDT, more general than those given in [2], are given below with data from direct numerical simulations (DNS) from Livescu [3]. For the reader not familiar with this model problem, homogenous Rayleigh–Taylor transition and turbulence, direct numerical simulations are presented with graphics, parameters, and analysis by Livescu and Ristorcelli [1,2] for initially equal amounts of fluids and for lower density ratios, $A=0.05\text{--}0.5$. For DNS of initially unequal amounts of fluids and at higher density ratios, i.e. A = 0.75, the reader is referred to a paper under review [4]. For an excellent description on the flow and its evolution, see Sec. 4 of [1]. Our concern here is with modelling the mixing process.

1.1. Problem statement

The Navier–Stokes-based moment equations for a second moment of the material field is (1) $\frac{D}{Dt}b=-{\epsilon }_b,$(1) where $b=-⟨\rho v⟩$ with v the fluctuating specific volume and ρ the density fluctuation about their respective means. Here the mixing rate is ${\epsilon }_b=⟨vd⟩$ where d is the dilatation [1,5,6]. The standard phenomenological ${\epsilon }_b$ model is wrong and there is little we have to guide us in a model of ${\epsilon }_b$. One can derive an equation for ${\epsilon }_b$ [7], but there are too many unclosed terms with new fluid physics.

To model the small scale mixing in terms of large scale quantities we invoke realisability as a modelling principle: we require the mixing process be related to an underlying PDF. We effectively add a conservation law – the conservation of the PDF – to the description of material mixing. We invoke a Fokker–Planck equation (FPE) that allows for a skewed mixing process (Figure 1). Beta distributions has been used in assumed PDF methods to describe turbulent mixing in combustion (Girimaji [8] and Fox [9]). Here we invoke, not a beta PDF, but a skewed beta process first derived by Bakosi and Ristorcelli [10].

Figure 1. Time evolution of the density PDF extracted from DNS [2,4] for various Atwood numbers and initial conditions. (Colour online.)

From the FPE [10] we derive the equations for various second-order moments. In particular, we derive the FPE analog to the Navier–Stokes (NS) version of ${\epsilon }_b$. The unknown ${\epsilon }_b$ of the NS problem is thus transferred to specifying the drift and diffusion coefficients of the FPE about which we know considerably more. By this procedure of synthesising second moment equations from the FPE and NSE perspectives, we create a formal mathematical procedure that shows how the small scales are to be slaved to the large scales if the cascade process is realisable. Our goal is a closure model for the cascade using a Fokker Planck equation for the mixing in a second moment equation. Our goal is not a PDF method to predict a PDF to compute in reacting flows; at best this work is a first step in that direction dealing with new mixing phenomena seen in variable density turbulence.

1.2. Outline of the paper

The paper proceeds as follows:

1. In Section 2 we discuss the peculiar issues of mixing in variable-density turbulence (VDT) using the mathematical structure of the VDT equations. Using the PDFs from the DNS we indicate the unusual skewness features of the VDT mixing process.

2. In Section 3, we present the model problem – Rayleigh Taylor in a box – and discuss its special features; non-equilibrium issues and skewness of the density and specific volume fields.

3. Our discussion in Section 4 then turns to some background on the passive scalar-hydrodynamic-cascade time-scale ratio and passive scalar models for mixing. The time-scale ratio will play an important role as a mix-rate metric in variable-density turbulence (VDT). We then further refine the problem goal.

4. In Section 5 we present the nomenclature peculiar to VDT and present some analytic results that are used in the derivations.

5. Sections 6 and 7 are the more theoretical sections of the work in which we synthesise information from the perspective of a Fokker–Planck equation with information from Reynolds-averaged moments of Navier–Stokes. After proceeding as far as possible with the mathematics, in Section 7 we add phenomenology, from the DNS, to the modelling problem.

6. We close with a discussion and summary in Sections 8 and 9.

2. Special features of variable-density mixing

The new physics of variable-density mixing can be seen by inspection of, e.g. the instantaneous mass fraction equation of species α, $Y_{\alpha }^{\ast }$, (2) $\left({\rho }^{\ast }{Y}_{\alpha }^{\ast }\right){,}_t+\left({\rho }^{\ast }{Y}_{\alpha }^{\ast }{U}_k^{\ast }\right){,}_k=-J_{k,k}^{\alpha },\alpha =1,\dots ,N,$(2) where ${\rho }^{\ast }$, $U_k^{\ast }$, and $J_k^{\alpha }$ are the instantaneous mixture-density, mixture velocity, and diffusive mass fluxes, respectively. The problem described in terms of the mass fraction variables in Equation (2(2) $\left({\rho }^{\ast }{Y}_{\alpha }^{\ast }\right){,}_t+\left({\rho }^{\ast }{Y}_{\alpha }^{\ast }{U}_k^{\ast }\right){,}_k=-J_{k,k}^{\alpha },\alpha =1,\dots ,N,$(2) ), is highly nonlinear: (1) in the advective terms of the scalar and momentum equations as ${\rho }^{\ast }=\rho \left({\mathbf{Y}}^{\mathbf{\ast }}\right)$, and (2) the diffusive flux, $J_k^{\alpha }$ (unlike if ${\rho }^{\ast }=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$) depends, very importantly, on ${\rho }^{\ast }$ and thus on all components of $Y_{\alpha }^{\ast }$.

2.1. Differential acceleration

For binary mixing with densities ${\rho }_1$ and ${\rho }_2$ the fluid density is given by (3) $\frac{1}{{\rho }^{\ast }}=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{1-Y_1^{\ast }}{{\rho }_2}=\frac{1+r{Y}_1^{\ast }}{{\rho }_2}r=\frac{{\rho }_2}{{\rho }_1}-1,A=\frac{{\rho }_2-{\rho }_1}{{\rho }_2+{\rho }_1}=\frac{r}{2+r},$(3) where A is the Atwood number. The acceleration and advective terms in the mass fraction and momentum equations are (4) $\frac{D}{\mathrm{D}t}\left(\frac{{\rho }_2}{1+r{Y}_1^{\ast }}{Y}_1^{\ast }\right),\frac{D}{\mathrm{D}t}\left(\frac{{\rho }_2}{1+r{Y}_1^{\ast }}{u}_i^{\ast }\right).$(4) If $r>10$ the terms vary considerably for $0<Y_1<1$. This is every different from the passive scalar mixing problem for which $r\to 0$ and $⟨{\rho }^2⟩/{\overline{\rho }}^2\ll 1$. Thus with density varying an order of magnitude in the acceleration terms one can expect some sizable differences from passive scalar mixing.

2.2. Differential diffusion

The diffusive flux in variable-density mixing is represented as (5) $J_k^{\alpha }=-{\rho }^{\ast }\left({\mathbf{Y}}^{\mathbf{\ast }}\right)\mathcal{D}{Y}_{\alpha ,k}^{\ast },$(5) where the diffusion coefficient $\mathcal{D}$ is assumed constant for all species. The diffusive flux of species 1 in binary mixing is (6) $J_k^1=-\frac{{\rho }_2}{1+r{Y}_1^{\ast }}\mathcal{D}{Y}_{1,k}^{\ast }.$(6) If $r>10$ the diffusive flux varies considerably for $0<Y_1<1$. In the binary case one can write the diffusion flux in terms of the mixture-density, (7) $J_k^1=-\frac{{\rho }_1{\rho }_2}{r}\mathcal{D}\frac{{\rho }^{\ast }{,}_k}{{\rho }^{\ast }}$(7) and the nonlinear nature of the diffusion is clear. This leads to phenomena not seen in the $\rho =\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$ case as heavy fluid mixing into light is qualitatively different than light mixing into heavy, and the PDF becomes skewed [11] even if it has no initial skewness. This mechanism does not occur in passive scalar mixing and can be seen in the equation for the density skewness given Section 3.4.

In Figure 1 the density PDFs from DNS [3] highlight these new and peculiar features of variable-density mixing. Figure 1(a) shows a low-Atwood-number mixing (A = 0.05) starting with equal amounts of each pure fluid. The initial conditions and the mixing process are symmetric. Figure 1(b) shows the Atwood number A = 0.75, reflecting a pure fluid density ratio of 7, with symmetric initial conditions (equal amounts of pure fluids). Immediately we see the generation of skewness in the PDF at high Atwood number. At the first time shown after the initial condition (green line) there is no pure light fluid left while there is plenty of pure heavy fluid, which reflects the earlier observations in [2]: due to inertia the heavy fluid mixes more slowly than the light fluid.

However, this is not simply an Atwood number effect. Comparing the green lines in Figure 1(c,d) we see that the mean density (or the amount of heavy fluid) plays a role. Observe the different amounts of heavy fluid (blue line, right-hand axis). There is a sizable amount of pure heavy fluid left (green and blue lines) and no light fluid left as is consistent with the high-Atwood inertial argument. However note that in Figure 1(d) the amount of pure light fluid (green line right-hand side), there is so little heavy fluid that it mixes rapidly. In all cases at late time mixing approaches a passive scalar Gaussian-like process as seen by the single more and more symmetric peak of the distribution.

3. The model mixing problem: Rayleigh–Taylor in a box

The problem we treat is a mixture of two different-density fluids, initially quiescent in a box, driven by the sudden application of a gravitational field. We call this homogeneous variable-density turbulence (HVDT) or homogeneous Rayleigh–Taylor mixing. DNS for this process have been given in [1,2]. The ‘equation of state’ is given by a binary material field as given above in (3(3) $\frac{1}{{\rho }^{\ast }}=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{1-Y_1^{\ast }}{{\rho }_2}=\frac{1+r{Y}_1^{\ast }}{{\rho }_2}r=\frac{{\rho }_2}{{\rho }_1}-1,A=\frac{{\rho }_2-{\rho }_1}{{\rho }_2+{\rho }_1}=\frac{r}{2+r},$(3) ). (8) ${\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1=\frac{2A}{1-A}.$(8) The Boussinesq, or passive-scalar, mixing problem corresponds to infinitesimal A and r. With increasing A (or r), under the application of a mean pressure gradient, the importance of differential acceleration of the different fluids, seen mathematically above, is demonstrated in the DNS simulations below. The DNS data set [3] of our model problem is at A = 0.05, 0.5, 0.75 corresponding to r = 0.1, 2.0, 6.0 and fluids whose densities are different by a ratio of 1.1, 3, and 7. The subscripts 2, and 1 reflect the heavy and the light fluids, respectively.

The Favre mean mass fraction, $\tilde{Y}$, or equivalently, the mean density, $\overline{\rho }$, is a second independent metric defining the problem. For this Rayleigh–Taylor in a box problem the Favre and Reynolds means, $\tilde{Y}$ and $\overline{\rho }$, are constant in time.

At t = 0 the PDF of the material field is approximately a double-delta of different heights corresponding to unequal amounts of the two pure fluids. The height of the spikes will be changed in different simulations to study the effects of the mean mass fraction and, importantly, the initial skewness of the PDF. At each density ratio there are three simulations in which the Favre mean mass fraction, $\tilde{Y}$, and thus the initial PDF skewnesses, are varied. The DNS database for the symmetric density PDF case is given in [1,2], while Ref. [3] discusses non-symmetric cases.

3.1. Special features of Rayleigh–Taylor in a box

The model problem is a non-equilibrium mixing process that begins as laminar, transitions to a turbulent flow, then decays as the potential energy diving the flow is reduced by the mixing process. The special features of the flow are (1) non-equilibrium and transitional, and (2) highly skewed PDFs at the initial state. The skewness is not monotonically decaying (as in the passive scalar case) and changes sign in some configurations.

3.2. The mixing process is non-equilibrium: a new time-scale

The initial turbulence energy is zero, increases rapidly as transition takes place, peaks in scaled time at $\tau \sim \sqrt{A}t=2$, then decays.

In Figure 2(a,b) we show the kinetic energy and the production-dissipation ratio, $P/\epsilon -1$, for the database of the 9 DNS runs. When $P/\epsilon =1$, near $\tau =\sqrt{A}t=2$, the flow transitions to decaying turbulence. We note that at early time the Boussinesq (A = 0.05) flow has the largest production-to-dissipation ratio due to the fact that the flow is not, at the outset, nonlinear enough to create a cascade to supply the dissipation. The higher-A flows go nonlinear and establish a cascade much earlier. During this early time the PDF of the material field is bimodal and becomes more Gaussian-like after $\tau \sim \sqrt{A}t=2$.

Figure 2. Time evolution of (a) the turbulent kinetic energy and (b) the production–dissipation ratio from DNS [3,4] for various Atwood numbers and initial conditions.

In equilibrium turbulence modelling one uses $\epsilon /k$ as the time-scale rate. This presumes a single important time-scale related to the established equilibrium cascade, where production is approximately equal to dissipation. In non-equilibrium turbulence the dissipation does not track the production or the cascade as inertial range begins development. It is useful to think in terms of an additional time-scale that represents the non-equilibrium nature of the flow and the development of a cascade. We identify the two time scales as (9) ${\tau }_{ne}\sim \frac{k}{\dot{k}}{\tau }_e\sim \frac{k}{\epsilon }.$(9) Thus besides the equilibrium time-scale, ${\tau }_e$, ${\tau }_{ne}$ represents a second time-scale of the flow during which the small scale dissipation does not track the cascade which track the production. During this period – when the turbulence energy is rapidly increasing and its dissipation is negligible – the turbulence serves primarily to stir but not mix the material field. This is an important difference from fully developed turbulence.

In keeping with the positivity of a time-scale a characteristic time-scale rate is chosen as (10) $\frac{1}{{\tau }_c}=\sqrt{{\left(\frac{\epsilon }{k}\right)}^2+{\left(\frac{\dot{k}}{k}\right)}^2}=\frac{\epsilon }{k}\sqrt{1+{\left(\frac{P}{\epsilon }-1\right)}^2},$(10) where we have used the equation $\dot{k}=P-\epsilon$ where P is the production, (the product of the mass flux and the pressure gradient), see [1,2]. The latter expression is relevant to homogeneous flows. For inhomogeneous flows there are additional transport terms in $\dot{k}$ that might be included. We shall call the term in the square-root the non-equilibrium factor and it will be used to scale the cascade during the non-equilibrium stages in which stirring dominates mixing. When the flow returns to equilibrium, $P\approx \epsilon$, the characteristic time-rate returns to its usual definition.

3.3. The density and specific volume PDF skewness

Mixing in variable-density fluids with large density-differences have some new and interesting features as shown in the density PDF in Figure 1. In [2] all the simulations begin with symmetric ‘double-delta’ PDFs for the density field. As time evolves the PDF becomes highly asymmetric as the heavy fluid – due to its inertia and resistance to deformation – mixes more slowly, leaving pure heavy fluid long after the pure light fluid has vanished. See also Figure 4 of [2].

3.4. The production of skewness

In homogeneous binary mixing the density skewness, $S_{\rho }=⟨{\rho }^3⟩/⟨{\rho }^2{⟩}^{3/2}$, is governed by [2] (11) $\frac{\mathrm{d}{S}_{\rho }}{\mathrm{d}t}=6\frac{⟨\rho {,}_k\rho {,}_k⟩}{R{e}_{0}S{c}_0}\left(\frac{S_{\rho }}{2⟨{\rho }^2⟩}-\frac{⟨\rho \rho {,}_k\rho {,}_k⟩}{⟨\rho {,}_k\rho {,}_k⟩⟨{\rho }^2⟩}\right).$(11) The equation is not of a simple relaxation type. There is a skewness source; even if the skewness is initially zero, skewness will be generated by $⟨\rho \rho {,}_k\rho {,}_k⟩$. The simplest examples of skewness generation in density-dependent mixing are the simulations of homogeneous variable-density turbulence by Livescu and Ristorcelli [1,2]. As shown in Figure 4 of [2], the initially symmetric distribution rapidly becomes skewed, reflecting the different and much slower mixing rate of the heavy fluid. This phenomena is not seen in passive scalar mixing where the mixing process if it starts with a symmetric PDF remains symmetric and the process is usefully modelled as a Gaussian mixing process. In passive scalar mixing - if the initial PDF is skewed the subsequent skewness merely decays monotonically. The behaviour is very different in VDT as the Figures 1 and 3 highlight.

Figure 3. Time evolution of (a) the third moment, (b) the skewness of the mixture density, (c) the third specific volume moment, (d) the skewness of the specific volume from DNS [2,4] for various Atwood numbers and initial conditions.

In Figure 3 we show the third moments and skewness measures for the density, ρ, and specific volume, v, fields. Several things are worthy of note: (1) the third moments of ρ versus v do not scale with each other, (2) have opposite signs, also importantly, (3) unlike the second moments, the third moments are not monotonic and can increase or decrease from the initial value, and (4) they change sign during the flow evolution.

4. Background and problem statement

Formal developments in moment closure approaches to VDT mixing are a fairly recent development [5]. It is useful to review how mixing is modelled in the passive scalar problem. The asymptotic behaviour of mixing in VDT is, of course, the passive scalar problem.

4.1. Passive scalar mixing – slaved to the hydrodynamical cascade

The conventional modelling assumption for small scale mixing of a passive scalar for turbulence with a developed inertial range is that the small scales are slaved to the large scales by a rate determining cascade and a model of the form ${\epsilon }_b={\epsilon }_b(\epsilon ,k,b)$ is assumed. In moment closure approaches to passive scalar mixing the mixing cascade ${\epsilon }_b$ of b reflecting the material field is (in a homogeneous field) modelled as (12) $\frac{D}{Dt}b=-{\epsilon }_b=-C_{b2}\frac{\epsilon }{k}b,$(12) where $C_{b2}$ is taken as a constant, $C_{b2}\approx 1.8$, as observed in empirical data. This mixing closure is a drag with time-scale set by the hydrodynamical cascade. The above closure can be shown to be derivable from the Fokker–Planck equation for a Langevin process – which models the mixing in Equation (1(1) $\frac{D}{Dt}b=-{\epsilon }_b,$(1) ) as a Gaussian process.

Such a closure implies that the time-scale ratio of the material cascade tracks the hydrodynamical cascade, (13) $C_{b2}=\frac{k}{\epsilon }\frac{\dot{b}}{b}\approx 1.8$(13) see [12] and references therein for alternate views and exact solutions that invalidate this model. In [12] it was demonstrated that the time-scale ratio $C_{b2}$ is better understood as the normalised mixing rate or the mixing rate in time normalised by $k/\epsilon$. The constant time-scale ratio appears a useful engineering approximation for stationary (or asymptotic) passive scalar mixing in a fully developed equilibrium turbulence. However, even in such a simple case, the problem is complex. Reynolds-number-dependent features have been observed [13] that are not consistent with the usual notions of an inviscid cascade, conceptually enforced in LES mixing models, [14]. We show the time-scale ratio in Figure 4 for the 9 DNS cases and see that the time-scale ratio varies substantially during the course of the flow evolution. Clearly, $C_{b2}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$ is not a good approximation at any time. Interestingly, the behaviour of the time-scale ratio is remarkably similar for all the Atwood numbers and different skewnesses and mean densities in the initial conditions. The 9 DNS are for three different Atwood numbers, each  with three different initial heights of the initial double-delta mass fraction distribution.

Figure 4. Time evolution of the time-scale ratio from DNS [2,4] for various Atwood numbers and initial conditions.

4.2. Problem definition and goal

Our goal is to develop an expression for the material mixing cascade, ${\epsilon }_b$, that reflects a variable-density mixing process over the course of which the PDF can be bimodal, uni-modal, skewed, and importantly, has bounded sample space.

There are three primary portions of the flow that a model for the cascade ${\epsilon }_b$ in VDT mixing must emulate: (1) the early time behaviour in which the PDF is bimodal and $P/\epsilon$ is large, implying that the mixing is primarily stirring, (2) the peak of the mixing rate as the cascade builds an inertial subrange and $P/\epsilon -1$ goes through zero and the decay process starts, and (3) the long-time behaviour in which the flow is primarily a passive scalar mixing. These three stages are seen in the time-evolution of the normalised mixing rates in Figure 4.

To close the small-scale dissipation term one makes arguments that there is a developed cascade and that the small scale mixing is slaved to, and can be written in terms of, quantities related to the large scales of the flow. This modelling closure is done heuristically.

Here we add a conservation law – the conservation of a PDF – as a constraint to the mixing description. We call this realisability which is also a formal requirement if the mixing process is modeled by a stochastic diffusion process. The application of this additional principle – that the statistical description of mixing must be related to a PDF – provides for a mathematically rigorous development that enables developing an expression for the enslavement of the small-scale to the large-scale statistics. There is, in fact, a valid unclosed PDF equation for the mixing process that (unfortunately) involves two-point correlations. Based on the heuristic that the small scales are slaved to the large scales, we use a single-point FPE governing the single-point PDF. To do this we presume that the mixing process is usefully modelled as a second-order diffusion process (in sample space) and derive a set of moment equations from the FPE. As will be seen, this simple requirement produces a cascade model for ${\epsilon }_b$ dependent on large scale variablesas (14) ${\epsilon }_b=f\left(b,\frac{\epsilon }{k},\frac{P}{\epsilon }-1,\tilde{Y},r,\overline{\rho },⟨{\rho }^2⟩,⟨{\rho }^3⟩\right).$(14) The form of Equation (14(14) ${\epsilon }_b=f\left(b,\frac{\epsilon }{k},\frac{P}{\epsilon }-1,\tilde{Y},r,\overline{\rho },⟨{\rho }^2⟩,⟨{\rho }^3⟩\right).$(14) ) comes from the mathematical requirement that the cascade forms consistently with the evolution of a Fokker–Planck diffusion process for the material PDF and satisfies variable-density constraints related to Equation (8(8) ${\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1=\frac{2A}{1-A}.$(8) ).

5. Variable density nomenclature and preliminaries

To develop the mathematics we use, as a model, the isothermal mixing between a heavy fluid of mass fraction $Y_2^{\ast }=1-Y_1^{\ast }$, and a light fluid of mass fraction $Y_1^{\ast }$, (15) $\frac{1}{{\rho }^{\ast }}=v^{\ast }=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{Y_2^{\ast }}{{\rho }_2},{\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}_1^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1>0,$(15) where ${\rho }_1,{\rho }_2$ are constant the densities of the pure light and heavy fluids, A is the Atwood number, and $v^{\ast }$ is the specific volume. The isothermal mixing problem, given by Equation (15(15) $\frac{1}{{\rho }^{\ast }}=v^{\ast }=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{Y_2^{\ast }}{{\rho }_2},{\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}_1^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1>0,$(15) ), is a specific class of variable-density mixing problems that Chassaing et al. [15] in their Section 6.3.1 or 4.1.3 call isothermal and isobaric mixing and implies the neglect of temperature and pressure fluctuations on the fluid density, ${\rho }^{ast}$. The statement in Equation (15(15) $\frac{1}{{\rho }^{\ast }}=v^{\ast }=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{Y_2^{\ast }}{{\rho }_2},{\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}_1^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1>0,$(15) ) is mathematically analogous to a moment-based treatment of turbulent combustion in which a single-progress-variable formulation is possible. In the case of premixed reactants the density is similarly written using the ideal gas law with r playing the role of a heat release parameter [16].

5.1. Reynolds and Favre decomposition nomenclature

The following nomenclature for density, mass fraction, and specific volume is used to treat the statistical problem: (16) $\begin{array}{rl}{\rho }^{\ast }& =\overline{\rho }+\rho y^{\ast }=\overline{Y}+y=\tilde{Y}+y^{\mathrm{\prime }\mathrm{\prime }};v^{\ast }=V+v\end{array}$(16) (17) $\begin{array}{rl}⟨y^{\ast }⟩& =\overline{Y}=\tilde{Y}+⟨y^{\mathrm{\prime }\mathrm{\prime }}⟩,\end{array}$(17) (18) $\begin{array}{rl}⟨y⟩& =0,⟨{\rho }^{\ast }{y}^{\mathrm{\prime }\mathrm{\prime }}⟩=0,\overline{\rho }⟨y^{\mathrm{\prime }\mathrm{\prime }}⟩=-⟨\rho y⟩,\end{array}$(18) (19) $\begin{array}{rl}\rho & ={\rho }^{\ast }-\overline{\rho },y=y^{\mathrm{\prime }\mathrm{\prime }}-⟨y^{\mathrm{\prime }\mathrm{\prime }}⟩.\end{array}$(19) Here the customary definitions are used: for first-order moments the overbar is a Reynolds ensemble average and the tilde a Favre average. The angle brackets used for higher order moments denote ensemble averages.

The following VDT relations and identities from [17] will be useful: (20) $V=\frac{1}{{\rho }_2}(1+r\overline{Y}),\overline{\rho }=\frac{{\rho }_2}{1+r\tilde{Y}}.$(20) The following definitions and relations between second-order moments will be used in the developments below (21) $\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}=\frac{⟨{\rho }^{\ast }{y}^{\mathrm{\prime }\mathrm{\prime }2}⟩}{\overline{\rho }},⟨v^2⟩={\left(\frac{r}{{\rho }_2}\right)}^2⟨y^2⟩={\left(\frac{r}{{\rho }_2}\right)}^2\left[⟨y^{\mathrm{\prime }\mathrm{\prime }2}⟩-⟨y^{\mathrm{\prime }\mathrm{\prime }}{⟩}^2\right].$(21) From the above relations for the mean and the variance it is clear that the mass fraction and the volume fraction are proportional. We follow the Besnard–Harlow–Rauenzahn (BHR) approach to VDT in [5]. In the BHR approach the primary second-order moment for the mixing process is $b=-⟨\rho v⟩$. BHR first identified b as an important moment due to the central role it plays in the production of the mass flux in pressure-gradient-driven flows.

We will need relations between b and $\overline{Y}$, $\tilde{Y}$, $⟨y^{\mathrm{\prime }\mathrm{\prime }}⟩$, and $\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}$. From [17] we have (22) $\overline{\rho }V=1-⟨\rho v⟩=1+b=\frac{1+r\overline{Y}}{1+r\tilde{Y}}.$(22) The above results can be re-arranged to obtain an identity that is frequently applied (23) $\overline{Y}=\tilde{Y}+\frac{1+r\tilde{Y}}{r}b.$(23) The term is not singular as $b\sim r^2$. The quantity b is a measure of the difference between Favre and Reynolds means (24) $b=r\frac{\overline{Y}-\tilde{Y}}{1+r\tilde{Y}}.$(24) As density and specific volume (or mass fraction) are inversely related to each other, the correlation is always negative. The result (25) $b=\frac{r^2}{(1+r\tilde{C}{)}^2}\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}={\left(\frac{r\overline{\rho }}{{\rho }_2}\right)}^2\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}},$(25) provides a connection to the Favre variance used in other areas of mixing, see also [17].

6. The Fokker–Plank equation and its moments

As a model for mixing in VDT we invoke a stochastic diffusion process that enforces bounded mass fractions and converges to a beta pdf. We require a process with bounded mass fractions that can develop non-zero skewness in the PDF, and satisfies a conservation law in the instantaneous (i.e. not only in a finite number of statistical moments) sense [18]. A FPE for the PDF of such mixing process for the mass or volume fraction $\hat{y}$ is (26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) We call the FPE (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) a beta process because its statistically stationary solution, its attractor or invariant, is a beta distribution. The function $\hat{y}(1-\hat{y})$ in the diffusion term keeps the process bounded in $\hat{y}$ space, [10,18]. Boundedness of the process is crucially important for non-dilute mixing. If $\hat{y}(1-\hat{y})$ where replaced with an expression that is not a function of $\hat{y}$, one would have an Ornstein-Uhlenbeck process with a Gaussian as its attractor, and a model that is used for symmetric, passive scalar mixing. Note that if the initial PDF f is arbitrary and not a beta PDF, Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) evolves the PDF such that $\hat{y}$ remains bounded and the solution, f, only asymptotically becomes a beta PDF with a non-zero skewness: a beta-convergent stochastic process. The stationary solution (if β, κ, and S are constants) to the FPE is the beta PDF, [10], with arbitrary non-zero skewness.

We call $\beta ,\kappa$, and S the drift, diffusion, and relaxation coefficients. They are, at this point, unspecified functions of the statistics of the PDF which are single-point, large scale quantities. Their specification will lead to a model for the small scale ${\epsilon }_b$ in terms of large scale quantities as is consistent with the usual notion that the small scales are slaved to the large scales. Requiring the mixing process to follow Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) produces a formal mathematical expression for how small scale mixing and its time-evolution are represented in terms of the large scales by a realisable mixing process. The physical mixing process is a small scale process involving two point statistics [19,20].

It is worth contrasting this work to an earlier Fokker–Planck approach to passive scalar mixing by Fox [21] based on models for complex mixing microstructure. Fox models the small two-point structure of the exact PDF equation from Navier-Stokes. Fox uses earlier experimental studies on mixing that use a narrative relating mixing to a lamella-scale microstructure and after some original and complex intuitive reasoning, scaling and normalisation invokes a bivariate Fokker–Planck equation for proxies for the concentration and scaled lamella length. Our approach is more pedestrian: we ignore the two-point structure of the small-scale mixing and assume that the small scales are slaved to the large scales such that a single point FPE is relevant.

At the level of principle Fox's [21] Fokker–Planck equation has a similar structure and solution that of our Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) with, in addition to the bimodal treatment, important differences. His mixing process treats zero-skewness, symmetric, passive scalar mixing; in which case our S becomes his $\overline{Y}$. The more general problem of mixing between unequal amounts of pure fluids (in a box) or in mixing between two different streams of fluids cannot be treated by the symmetric beta model even for passive scalar mixing and in practical problems involving inhomogeneity such as the planetary boundary layer or the bubbles-, and spikes-skewness in the Rayleigh–Taylor mixing layer. In its simplest form this is the case of mixing in a box of homogenous stationary turbulence in which the initial PDF is two unequal-height delta spikes.

In such cases the skewness (Figure 3) of the PDF – in variable density turbulence (VDT) – plays a hugely important role which, for our current discussion, results in different drift coefficients. In the beta process Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) $S-\overline{Y}\ne 0$ is responsible for generating the skewness of the PDF. Skewness is important in general passive scalar mixing processes that can be generated by mixing between two fluids with unequal mass fractions.

In summary, the differences between Fox's problem and approach and ours are as follows.

• Fox treats symmetric passive scalar mixing with zero skewness. We treat active scalar mixing.

• We treat the case in which skewness may arise due to initial conditions or due to variable-density effects. In our FPE skewness will be connected to $S-\overline{Y}\ne 0$.

• For mixing in VDT we see that symmetric initial conditions of the PDF (no initial skewness) become skewed. The behaviour of skewness in VDT is not monotonic as it is in the case of passive scalar mixing.

• Fox's analysis is useful, interesting, and connects to other methods for passive scalar mixing [22], and is mathematically important as a special limiting case for more general models to attain.

6.1. The mean density equation from Fokker–Planck

Our benchmark mixing problem is the homogeneous, buoyantly-driven VDT problem; turbulence in box for which the mean density and mean Favre mass fraction are constant. The FPE has three unknown coefficients, S, β, and κ. In this section we apply the constraint of constant mean density to find that S is a function of $\beta ,\kappa$, and the closure problem is reduced to the problem of defining $\beta ,\kappa$ as S is not an independent object.

Using the definition of ${\rho }^{\ast }\left(\hat{y}\right)$ and taking the moment of Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) we find that the mean density equation is governed by (27) $\frac{D}{Dt}\overline{\rho }=\frac{\beta }{2}⟨{\rho }_{y^{\ast }}^{\ast }(S-y^{\ast })⟩+\frac{\kappa }{2}⟨{\rho }_{y^{\ast }{y}^{\ast }}^{\ast }{y}^{\ast }(1-y^{\ast })⟩.$(27) The angle brackets reflect ensemble averages (moments) and ${\rho }_{y^{\ast }}^{\ast }=\mathrm{d}{\rho }^{\ast }/\mathrm{d}{y}^{\ast }$. The above equation is general and reflects the statistical nature of the process and not any physics. We now add the physics of the variable-density relations, Equation (15(15) $\frac{1}{{\rho }^{\ast }}=v^{\ast }=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{Y_2^{\ast }}{{\rho }_2},{\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}_1^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1>0,$(15) ). Note the following relations, derivable from Equation (15(15) $\frac{1}{{\rho }^{\ast }}=v^{\ast }=\frac{Y_1^{\ast }}{{\rho }_1}+\frac{Y_2^{\ast }}{{\rho }_2},{\rho }^{\ast }=\frac{{\rho }_2}{1+r{Y}_1^{\ast }},r=\frac{{\rho }_2}{{\rho }_1}-1>0,$(15) ) (28) ${\rho }_{y^{\ast }}^{\ast }=-\frac{r}{{\rho }_2}{\rho }^{\ast 2}{\rho }_{y^{\ast }{y}^{\ast }}^{\ast }=2\frac{r^2}{{\rho }_2^2}{\rho }^{\ast 3}.$(28) For the current model problem $\overline{\rho }$ is constant and the mean density equation becomes (29) $-\frac{\beta }{2}⟨{\rho }_{y^{\ast }}^{\ast }(S-y^{\ast })⟩=\frac{\kappa }{2}⟨{\rho }_{y^{\ast }{y}^{\ast }}^{\ast }{y}^{\ast }(1-y^{\ast })⟩$(29) and an exact expression for S in terms of β and κ and various moments is determined. Performing the indicated operations specifies S as (30) $(1+rS)\frac{\widehat{{\rho }^2}}{\hat{\rho }}=1-\frac{2\kappa }{\beta }\left[1-(2+r)\frac{\widehat{{\rho }^2}}{\hat{\rho }}+(1+r)\frac{\widehat{{\rho }^3}}{\hat{\rho }}\right]=1-\frac{2\kappa }{\beta }G.$(30) The higher order moments above can be related to their dimensional equivalents as (31) $\begin{array}{rl}\hat{\rho }& =\frac{\overline{\rho }}{{\rho }_2}=\frac{1}{1+r\tilde{Y}},\end{array}$(31) (32) $\begin{array}{rl}\widehat{{\rho }^2}& =\frac{⟨{\rho }^{\ast 2}⟩}{{\rho }_2^2}=\frac{{\overline{\rho }}^2}{{\rho }_2^2}\left[1+\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\right],\end{array}$(32) (33) $\begin{array}{rl}\widehat{{\rho }^3}& =\frac{⟨{\rho }^{\ast 3}⟩}{{\rho }_2^3}=\frac{{\overline{\rho }}^3}{{\rho }_2^3}\left[1+3\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}+\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}\right].\end{array}$(33) Our Monte Carlo (MC) simulations verify that the above choice of S in Equation (30(30) $(1+rS)\frac{\widehat{{\rho }^2}}{\hat{\rho }}=1-\frac{2\kappa }{\beta }\left[1-(2+r)\frac{\widehat{{\rho }^2}}{\hat{\rho }}+(1+r)\frac{\widehat{{\rho }^3}}{\hat{\rho }}\right]=1-\frac{2\kappa }{\beta }G.$(30) ) produces a constant mean density. Using the above definitions G is written as (34) $G=\frac{1}{(1+r\tilde{Y}{)}^2}\left\{-r^2\tilde{Y}(1-\tilde{Y})+\left[1+2r-(2+r)r\tilde{Y}\right]\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}+(1+r)\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}\right\}.$(34) The FPE coefficient S is now given in terms of the drift and diffusion coefficients and moments of the density field (35) $1+rS=\frac{\hat{\rho }}{\widehat{{\rho }^2}}\left[1-\frac{2\kappa }{\beta }G\right].$(35) With G given by Equation (34(34) $G=\frac{1}{(1+r\tilde{Y}{)}^2}\left\{-r^2\tilde{Y}(1-\tilde{Y})+\left[1+2r-(2+r)r\tilde{Y}\right]\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}+(1+r)\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}\right\}.$(34) ) the FPE coefficient S reflects the fluid physics of VDT through Equation (28(28) ${\rho }_{y^{\ast }}^{\ast }=-\frac{r}{{\rho }_2}{\rho }^{\ast 2}{\rho }_{y^{\ast }{y}^{\ast }}^{\ast }=2\frac{r^2}{{\rho }_2^2}{\rho }^{\ast 3}.$(28) ). The problem of the small scale ${\epsilon }_b$ has been replaced with the specification of $\beta ,\kappa$ which will be functions of single point large scale statistics.

A few limits are of interest. The Boussinesq limit $(r\to 0)$ is (36) $G_B=-\left[r^2\tilde{Y}(1-\tilde{Y})-\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}-\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}\right]$(36) and thus $G\sim r^2\sim 0$. The asymptotic fully mixed state for which the density moments vanish is (37) $G_{\mathrm{\infty }}=-\frac{r^2\tilde{Y}(1-\tilde{Y})}{(1+r\tilde{Y}{)}^2}$(37) for arbitrary $r\ge 0$. Using the results from Appendix 2 in the no-mix state one has (38) $\begin{array}{rl}\frac{⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^2}& =\frac{r^2}{1+r}\tilde{Y}(1-\tilde{Y})\end{array}$(38) (39) $\begin{array}{rl}\frac{⟨{\rho }^3{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^3}& =-\frac{r}{1+r}\frac{⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^2}[1-(2+r\left)\tilde{Y}\right].\end{array}$(39) (40) $\begin{array}{rl}{G}_{\mathrm{n}\mathrm{m}}& =-\frac{r^3}{1+r}\tilde{Y}(1-\tilde{Y})=-r\frac{⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^2}.\end{array}$(40)

6.2. The Fokker–Planck $⟨{\mathbf{y}}^{\mathbf{2}}⟩$ equation

Having the FPE it is now possible to generate evolution equations for the cascade of any of the moments. The equation for the variance of the Reynolds mass fraction fluctuation, derived from the FPE, is (41) $\begin{array}{rl}\frac{D}{Dt}⟨y^2⟩& =-\beta ⟨y^2⟩+\kappa \overline{Y}(1-\overline{Y}){\mathrm{\Theta }}_{yt}\end{array}$(41) (42) $\begin{array}{rl}{\mathrm{\Theta }}_{yt}& =1-\frac{⟨y^2⟩}{\overline{Y}(1-\overline{Y})}\overline{Y}=\tilde{Y}+\frac{1+r\tilde{Y}}{r}b.\end{array}$(42) The $⟨y^2⟩$ equation is the simplest second moment equation from the FPE that reflects the variable-density mixing process. To insure the long term fully mixed state and that $⟨y^2⟩\to 0$ from above requires that the diffusion coefficient scales as (at the very least) $\kappa ={\kappa }_y^{\mathrm{\prime }}⟨y^2⟩$. The variance equation is then written as (43) $\frac{D}{Dt}⟨y^2⟩=-\beta \left[1-\frac{{\kappa }_y^{\mathrm{\prime }}}{\beta }\overline{Y}(1-\overline{Y}){\mathrm{\Theta }}_{yt}\right]⟨y^2⟩=-{\epsilon }_y.$(43) This is our first interesting result. The first equality shows how ${\epsilon }_y$ might be modelled if mixing proceeded as a realisable diffusion process (with a beta attractor for the mass fractions). We see that the mix rate depends on a measure of the mix state, ${\mathrm{\Theta }}_{yt}$. Realisability for mixing requires a monotone decay of $⟨y^2⟩$, which can be ensured if $\beta \ge {\kappa }_y^{\mathrm{\prime }}/4$.

6.3. The Fokker–Planck $\mathbf{b}$ equation

Our primary goal is to develop a model for the mixing cascade ${\epsilon }_b$ in Equation (1(1) $\frac{D}{Dt}b=-{\epsilon }_b,$(1) ), suitable for early time mixing, non-equilibrium conditions, as well as asymptotically. The FPE moment equation for $b=-⟨\rho v⟩=r/{\rho }_2⟨\rho y⟩$ is found by taking the moments of Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) with ρ and v: (44) $\frac{D}{Dt}⟨\rho v⟩=\frac{\beta }{2}⟨(v{\rho }_{y^{\ast }}^{\ast }+\rho v_{y^{\ast }}^{\ast })(S-y^{\ast })⟩+\frac{\kappa }{2}⟨(v{\rho }_{y^{\ast }{y}^{\ast }}^{\ast }+2v_{y^{\ast }}^{\ast }{\rho }_{y^{\ast }}^{\ast })y^{\ast }(1-y^{\ast })⟩.$(44) If we apply Equations (28(28) ${\rho }_{y^{\ast }}^{\ast }=-\frac{r}{{\rho }_2}{\rho }^{\ast 2}{\rho }_{y^{\ast }{y}^{\ast }}^{\ast }=2\frac{r^2}{{\rho }_2^2}{\rho }^{\ast 3}.$(28) ) and (30(30) $(1+rS)\frac{\widehat{{\rho }^2}}{\hat{\rho }}=1-\frac{2\kappa }{\beta }\left[1-(2+r)\frac{\widehat{{\rho }^2}}{\hat{\rho }}+(1+r)\frac{\widehat{{\rho }^3}}{\hat{\rho }}\right]=1-\frac{2\kappa }{\beta }G.$(30) ) and take the moments we find (45) $\frac{D}{Dt}b=-\frac{\beta }{2}\left[b+\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2+⟨{\rho }^2⟩}+\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2+⟨{\rho }^2⟩}\frac{2\kappa }{\beta }G\right].$(45) In anticipation of subsequent derivations the equation is rewritten as (46) $\frac{D}{Dt}b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)b+\frac{2\kappa }{\beta }G\right]\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\eta b=\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}.$(46) The analysis of the $⟨y^2⟩$ equation indicated that κ vanishes asymptotically. Without loss of generality we take $\kappa ={\kappa }^{\mathrm{\prime }}b$ and (47) $-\frac{D}{Dt}b={\epsilon }_b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }G\right]b.$(47) Equation (47(47) $-\frac{D}{Dt}b={\epsilon }_b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }G\right]b.$(47) ) is a FPE-based model for ${\epsilon }_b$ and is our central result: it is a large scale parameterisation of small scale mixing based on realisability as linked to the FPE diffusion process.¹

Parenthetically, for a homogeneous problem, for which $\tilde{Y}$ is constant, the first moment equation for the Reynolds-averaged mass fraction (48) $\frac{D}{Dt}\overline{Y}=-\frac{\beta }{2}[\overline{Y}-S]$(48) can be used to derive an equation for b by applying Equation (23(23) $\overline{Y}=\tilde{Y}+\frac{1+r\tilde{Y}}{r}b.$(23) ). The result is (49) $\frac{D}{Dt}b=-\frac{\beta }{2}[b+r\hat{\rho }(\tilde{Y}-S\left)\right].$(49) Equation (47(47) $-\frac{D}{Dt}b={\epsilon }_b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }G\right]b.$(47) ) predicts, a priori, a curious phenomena: within a Markovian framework the asymptotic mixing rate is potentially a function of the initial conditions. The asymptotic b equation is (50) $\frac{D}{Dt}b=-\beta \left[1-r^2\frac{\tilde{Y}(1-\tilde{Y})}{(1+r\tilde{Y}{)}^2}\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }\right]b$(50) suggesting the possibility that the mean initial mass fraction and density ratio determine the asymptotic mixing rate. This is not seen in the DNS data, suggesting that the second term is small compared to the first. The comparison is not conclusive as the DNS data represents a single realisation of the flow for each Atwood number and each $\tilde{Y}$; the theory above is for an ensemble of realisations. Using the above results the next step is to determine $\beta ,\kappa$, and S in order to create a realisable skewed non-Gaussian process for the mixing cascade, ${\epsilon }_b$ which we then we test against DNS data using Monte Carlo simulations.

7. Determining the drift $\beta$ and diffusion $\kappa$ coefficients

The Fokker–Plank b equation (47(47) $-\frac{D}{Dt}b={\epsilon }_b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }G\right]b.$(47) ) shows how the mixing depends on the drift and diffusion coefficients. We have proceeded as far as the mathematics will allow for the drift β and diffusion κ coefficients. Our procedure now becomes phenomenological. We have two non-dimensional second moment metrics that describe the mixing process: the time-scale ratio or mix rate metric and a mix state metric. We begin by looking for collapses of the data with variable-density metrics.

7.1. The time-scale ratio is a mix rate metric

In the passive scalar turbulence moment closure problem the conventional model for the mixing is a drag model with time-scale set by the ‘eddy turnover’ rate (51) ${\epsilon }_b=C_{b2}\frac{\epsilon }{k}b,$(51) where $C_{b2}\approx 1.8$ and, as discussed above, implies a Gaussian model for the mixing process. At large times we expect passive-scalar Gaussian-like mixing and β to scale as (52) $\beta \sim \frac{\epsilon }{k}.$(52) It should be made clear that the passive scalar limit with $C_{b2}\approx 1.8$ is an approximation that is known to not be consistent with observed behaviours, [12,13] but qualitatively useful for engineering models.

There are a number of reasons to work with some form of the time-scale ratio [12] (53) $-\frac{k}{\epsilon }\frac{\dot{b}}{b}$(53) as a ‘mix rate’ metric reflecting as it does the coupling between the mixing and hydrodynamic cascades. In non-dimensional time, $\tau =\left(\epsilon /k\right)t$ the time-scale ratio is the mixing rate [12].

The FPE b equation (47(47) $-\frac{D}{Dt}b={\epsilon }_b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }G\right]b.$(47) ) suggests that the time-scale ratio of Figure 4 might better collapse when scaled as (54) $tsc{r}_{VD}=-\left(1+\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\right)\frac{k}{\epsilon }\frac{\dot{b}}{b}.$(54) Figure 5 shows this to be the case and we have an indication supporting the FPE strategy. The scatter in the time-scale data is reduced near the peak and also early time and during flow transition. This scaling collapses the time-scale ratio to approximately unity at the initial time.

Figure 5. Time evolution of the time-scale ratio, scaled as in Equation (54(54) $tsc{r}_{VD}=-\left(1+\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\right)\frac{k}{\epsilon }\frac{\dot{b}}{b}.$(54) ) from DNS [2,4] for various Atwood numbers and initial conditions.

7.2. A second-order mix state metric for variable-density turbulence

There are several second-order mix state metrics [17]. As the primary second-order moment in the BHR strategy in VDT is b [17] and our primary goal is a model for ${\epsilon }_b$, it is consistent to choose a second moment mix state metric as (55) $\mathrm{\Theta }=1-\frac{b}{b_{\mathrm{n}\mathrm{m}}},$(55) where $b_{\mathrm{n}\mathrm{m}}$ is given in Appendix 2. We now query the DNS if there is a variable-density scaling, ${\mathrm{\Theta }}^{\mathrm{\prime }}=f\mathrm{\Theta }$, that better collapses the data. We will construct variable-density functions and look to see if we can find more universal behaviour. There are several VDT functions that naturally appear from the Fokker–Planck analysis: (56) $\begin{array}{rl}{f}_1& =\frac{1}{2}(1+\eta b+\eta ),\end{array}$(56) (57) $\begin{array}{rl}{f}_2& =1+b,\end{array}$(57) (58) $\begin{array}{rl}{f}_3& =1+\eta b,\end{array}$(58) (59) $\begin{array}{rl}{f}_4& =\eta ,\end{array}$(59) (60) $\begin{array}{rl}{f}_5& =\frac{1}{2}\left(\frac{1+\eta b+\eta }{1+\eta b}\right).\end{array}$(60) We expect these functions to collapse the mixing process during the stage when the variable-density effects are important. We have arranged them in forms that become, asymptotically, unity and thus reflect the late-stage passive-scalar-mixing limit. We find that the best collapse for Θ is with $f_5$ thus we adopt a new variable-density mix metric ${\mathrm{\Theta }}^{\mathrm{\prime }}$ as (61) ${\mathrm{\Theta }}^{\mathrm{\prime }}=\frac{1}{2}\left(1+\frac{\eta }{1+\eta b}\right)\mathrm{\Theta }.$(61) The quantities Θ and $\mathrm{\Theta }(1-\mathrm{\Theta })$ and their VDT equivalents ${\mathrm{\Theta }}^{\mathrm{\prime }}$ and ${\mathrm{\Theta }}^{\mathrm{\prime }}(1-{\mathrm{\Theta }}^{\mathrm{\prime }})$ are shown in Figure 6. The collapse will be useful.

Figure 6. Time evolution of the second-order mix state metrics, given by Equations (55(55) $\mathrm{\Theta }=1-\frac{b}{b_{\mathrm{n}\mathrm{m}}},$(55) ) and (61(61) ${\mathrm{\Theta }}^{\mathrm{\prime }}=\frac{1}{2}\left(1+\frac{\eta }{1+\eta b}\right)\mathrm{\Theta }.$(61) ) from DNS [2,4] for various Atwood numbers and initial conditions.

7.3. Closing the drift and diffusion coefficients

We take β, and for simplicity κ, to scale as (62) $\beta ={\beta }_1\frac{\epsilon }{k},\kappa ={\kappa }^{\mathrm{\prime }}b={\kappa }_1\frac{\epsilon }{k}{b}^2.$(62) It is possible to show that ${\kappa }_1$ is second-order (in powers of second moments) relative to ${\beta }_1$.

The b Equation (47(47) $-\frac{D}{Dt}b={\epsilon }_b=-\frac{\beta }{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }^{\mathrm{\prime }}}{\beta }G\right]b.$(47) ), reflecting a realisable non-Gaussian mixing process, becomes (63) $-\frac{D}{Dt}b={\epsilon }_b=\frac{{\beta }_1}{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }_1}{{\beta }_1}Gb\right]\frac{\epsilon }{k}b.$(63) For simplicity and because our interest is in capturing the second moments, we take ${\kappa }_1/{\beta }_1$ as constant. To capture the behaviour of other higher statistics we expect that a deeper study of κ is required. For now our interest is a proof of concept adequately representing second-order statistics. A more complex dependence of $C_{b2}$ is now established: (64) $C_{b2}=\frac{{\beta }_1}{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }_1}{{\beta }_1}Gb\right].$(64) Equation (64(64) $C_{b2}=\frac{{\beta }_1}{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }_1}{{\beta }_1}Gb\right].$(64) ) can be written in terms of the variable-density time-scale ratio as (65) $tsc{r}_{VD}=-(1+\eta b)\frac{k}{\epsilon }\frac{\dot{b}}{b}=\frac{1}{2}{\beta }_1(1+\eta +\eta b)+{\kappa }_{1}Gb\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\eta b=\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}$(65) and thus (66) $\frac{2tsc{r}_{VD}}{1+\eta +\eta b}={\beta }_1+\frac{2{\kappa }_{1}b}{1+\eta +\eta b}G.$(66) To leading order, we can think of ${\beta }_1$ as the time-scale ratio (the non-dimensional mix rate).²

The DNS data shows that the flow can be divided into 3 regimes: early time stirring and transitional mixing, intermediate-time mixing where the kinetic energy is maximum, and a late time decay and asymptotic Gaussian mixing.

1. Early time, $tsc{r}_{VD}\approx 1$, $\sqrt{A}t<1$: In this regime the time-scale ratio is consistent with the idea that as the flow begins the fluids are unmixed and primarily undergoing stirring, Θ is small and the time-scale of the material field tracks the time-scale of the hydrodynamic field, $tsc{r}_{VD}\approx 1$. The PDF bimodal and can be highly skewed.

2. Maximum kinetic energy mixing, $\sqrt{A}t\approx 2$: Maximum kinetic energy is coincident with maximum mixing at time $\sqrt{A}t\approx 2$ at which $\mathrm{\Theta }=1/2$. In this regime the PDF develops a single peak, the skewness (depending on the initial conditions) is substantially reduced.

3. Asymptotic late stage mixing, $\sqrt{A}t\sim 5$: $\mathrm{\Theta }\approx 1$ and for which the variable-density time-scale ratio $tsc{r}_{VD}\approx 2$. In this regime the PDF is well on its way to becoming more Gaussian-like.

A simple phenomenological approximation consistent with these notions and relating the mix rate metric to the mix state metric is (67) ${\beta }_1=\frac{2}{1+\eta +\eta b}\left[{\beta }_{10}^{\mathrm{\prime }}+{\beta }_2^{\mathrm{\prime }}{\mathrm{\Theta }}^{\mathrm{\prime }}+{\beta }_3^{\mathrm{\prime }}{\mathrm{\Theta }}^{\mathrm{\prime }}(1-{\mathrm{\Theta }}^{\mathrm{\prime }})\right].$(67) In this expression ${\beta }_{10}^{\mathrm{\prime }}$ captures the time-scale ratio behaviour at early time when there is little mixing and ${\mathrm{\Theta }}^{\mathrm{\prime }}$ is small. The ${\beta }_3^{\mathrm{\prime }}$ term reflects a non-equilibrium peak where ${\mathrm{\Theta }}^{\mathrm{\prime }}\approx 1/2$. The ${\beta }_2^{\mathrm{\prime }}$ term captures the asymptotic passive scalar limit. The various β's are unknown and we model them according to the following notions consistent with the DNS data.

• We expect ${\beta }_2^{\mathrm{\prime }}$ and ${\beta }_3^{\mathrm{\prime }}$ to be small in the highly non-equilibrium early time stirring phase when little mixing occurs and turbulence production in the large scales dominates dissipation. We expect the mixing rate to be suppressed according to the magnitude of the non-equilibrium metric. During this non-equilibrium time-period the relevant time-scale rate is $\dot{k}/k$ as discussed above.

• The most non-equilibrium portion of the mixing is near where $\mathrm{\Theta }(1-\mathrm{\Theta })$ is maximum. It is where k is largest and is demarcated as where $P=\epsilon$ after which the turbulence decays.

• At late time the flow is quasi-equilibrium and $(P/\epsilon -1{)}^2$ is small. We expect ${\beta }_1\to 1.8$ and thus ${\beta }_{10}^{\mathrm{\prime }}+{\beta }_2^{\mathrm{\prime }}\to 2$ as ${\mathrm{\Theta }}^{\mathrm{\prime }}\to 1$.

• We have seen that scaling the time by $\sqrt{A}$ collapses the different-Atwood-number flows. We need to scale time during the non-equilibrium portion of the flow so we choose a proxy for the Atwood number that reflects the evolution of the mixing – a local-time-scale Atwood-number proxy is $A\left(t\right)\sim \sqrt{b}$. Thus we incorporate the effect of the current Atwood number with the factor $\sqrt{A\left(t\right)}\sim b^{1/4}$ as suggested from the no-mix result, $A\sim r\sim \sqrt{b}$. While we have provided physical reasoning as to why this is plausible, our experiments have shown that this causes the MC simulation results to compare move favourably with DNS compared to other scalings.

The above points suggest (68) ${\beta }_1=\frac{2}{1+\eta +\eta b}\left[{\beta }_{10}+{\beta }_2\frac{{\mathrm{\Theta }}^{\mathrm{\prime }}}{\sqrt{1+(P/\epsilon -1{)}^2{b}^{1/4}}}+{\beta }_3\frac{4{\mathrm{\Theta }}^{\mathrm{\prime }}(1-{\mathrm{\Theta }}^{\mathrm{\prime }})}{\sqrt{1+(P/\epsilon -1{)}^2{b}^{1/4}}}.\right]$(68) We have, in principle, posed a relation between the time-scale ratio – a mix rate metric and a VDT mix state metric, ${\mathrm{\Theta }}^{\mathrm{\prime }}$ based on intuitive arguments and DNS data with a mathematical structure imposed by a realisability requirement.

7.4. Monte Carlo solutions

Every FPE, such as Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ), is statistically equivalent to a system of stochastic differential equations [23]. We develop our ideas using Monte Carlo (MC) simulations in which we solve, for an ensemble of Lagrangian particles, the stochastic equation (69) $D{y}^{\ast }\left(t\right)=\frac{1}{2}\beta (S-y^{\ast })\mathrm{d}t+\sqrt{\kappa y^{\ast }(1-y^{\ast })}\mathrm{d}W\left(t\right),$(69) where $\mathrm{d}W$ is a Wiener process with independent increments, zero mean and $\mathrm{d}t$ variance, governing each stochastic particle [10] and β, κ, and S are as given above. The complete model is given in Appendix 1.

Figure 7. Time evolution of two normalised mix rates (left), the normalised mix state, b, (down-right) and its rate (up-right) for A = 0.05. Dotted lines – DNS data [2,4], solid lines – Monte Carlo model predictions.

Figure 8. Time evolution of two normalised mix rates (left), the normalised mix state, b, (down-right) and its rate (up-right) for A = 0.5. Dotted lines – DNS data [2,4], solid lines – Monte Carlo model predictions.

Figure 9. Time evolution of two normalised mix rates (left), the normalised mix state, b, (down-right) and its rate (up-right) for A = 0.75. Dotted lines – DNS data [2,4], solid lines – Monte Carlo model predictions.

The model predictions using MC are compared to DNS data in Figures 7^{Figure 8}–9 for three different Atwood numbers each for three different initial conditions. The MC drift and diffusion are specified as above with the hydrodynamical quantities, ϵ, k, and P taken from DNS. We have plotted the two time-scale ratios as well as the normalised b and it's time-scale. It should be emphasised that we have predicted the second moments of the mixing from the quiescent (stationary stochastic blobs of fluid) initial condition, through transition, to maximum kinetic energy, and then decay. The agreement is good especially when one considers we are comparing single DNS realisations with a theory for ensemble realisations.

8. General discussion

Having determined the drift and diffusion coefficients of a variable-density mixing process, we have also determined the mixing rate of any second-order moment of the material field. In Figures 10^{Figure 11}–12 we plot the time-evolutions of the time-scale ratios of and their normalised fields for $⟨y^2⟩$ and $⟨{\rho }^2⟩$ as predicted by the model using MC compared to DNS data. The same drift and diffusion coefficients also work for these other second-order fields.

Figure 10. Time evolution of the normalised mix rates of the specific volume (up-left), the mixture density (down-left), and the normalised mix states as measured by the specific volume variance (up-right) and the mixture density (down-right) for A = 0.05. Dotted lines – DNS data [2,4], solid lines – Monte Carlo model predictions.

Figure 11. Time evolution of the normalised mix rates of the specific volume (up-left), the mixture density (down-left), and the normalised mix states as measured by the specific volume variance (up-right) and the mixture density (down-right) for A = 0.5. Dotted lines – DNS data [2,4], solid lines – Monte Carlo model predictions.

Figure 12. Time evolution of the normalised mix rates of the specific volume (up-left), the mixture density (down-left), and the normalised mix states as measured by the specific volume variance (up-right) and the mixture density (down-right) for A = 0.75. Dotted lines – DNS data [2,4], solid lines – Monte Carlo model predictions.

Additionally one has an expression for the mixing cascade of any quantity one can derive an equation from FPE (not a trivial task for many moments). Equation (64(64) $C_{b2}=\frac{{\beta }_1}{2(1+\eta b)}\left[(1+\eta +\eta b)+\frac{2{\kappa }_1}{{\beta }_1}Gb\right].$(64) ) is a moment closure model for the mixing of b and equivalently (in homogeneous flow) $\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}$ due to Equation (25(25) $b=\frac{r^2}{(1+r\tilde{C}{)}^2}\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}={\left(\frac{r\overline{\rho }}{{\rho }_2}\right)}^2\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}},$(25) ). From the FPE for the variance $⟨y^2⟩$, Equation (43(43) $\frac{D}{Dt}⟨y^2⟩=-\beta \left[1-\frac{{\kappa }_y^{\mathrm{\prime }}}{\beta }\overline{Y}(1-\overline{Y}){\mathrm{\Theta }}_{yt}\right]⟨y^2⟩=-{\epsilon }_y.$(43) ), we have a moment expression for the mixing rate of $⟨y^2⟩$ as (70) ${\epsilon }_y={\beta }_1\left[1+\frac{{\kappa }_1}{{\beta }_1}\frac{{\mathrm{\Theta }}_{yt}}{{\mathrm{\Theta }}_{yt}-1}{b}^2\right]\frac{\epsilon }{k}⟨y^2⟩$(70) with the drift and diffusion coefficients as determined. Second-order moment models for the mixing of $⟨{\rho }^2⟩$ and for $\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}$ in inhomogeneous flows requires the FPE form equation.

8.1. Ensemble models versus individual realisations

The 9 DNS date we have used represent  individual realisations of the flow for 3 Atwood numbers and 3 different mean densities at each Atwood number at a fairly low ${256}^3$ grid resolution. In comparing MC to DNS it is crucial to keep in mind that the Fokker–Plank development for the cascade model is for an ensemble of realisations. One cannot expect the MC ensemble to reflect the individual DNS realisations and as such we do not expect perfect agreements. We have available data from three additional DNS of this problem with ${1024}^3$ grid resolution and in a subsequent work (now in progress) following Schwarzkopf et al. [24] we are also studying the model's performance in inhomogeneous flows.

8.2. Comment on application of Monte Carlo simulations as a modelling tool

Statistical approaches to turbulence face the moment closure problem that requires expressions for unclosed moments in terms of known moments. This modelling process is generally done intuitively and is usually only tested by attempting to compute real physics with the model and seeing if it does better for some statistical moments for which one has data. Here our goal of using MC has been to aid model development. We have conducted MC simulations of Equation (26(26) $\frac{D}{Dt}f=-\frac{1}{2}\frac{\mathrm{\partial }}{\mathrm{\partial }\hat{y}}\left[\beta \right(S-\hat{y}\left)f\right]+\frac{1}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\hat{y}}^2}\left[\kappa \hat{y}\right(1-\hat{y}\left)f\right].$(26) ) which, by definition, evolves a valid PDF from which we numerically estimatetatistical moments of the density, specific volume, or mass fraction fields and thus can test closure assumptions directly, without a moment closure or DNS calculation of the coupled hydrodynamics + mix problem. We find this particularly useful in VDT where the relation between moments of the specific volume and density fields involves non-Gaussian moments of all orders [17] and thus guessing at closures has many challenges.

9. Concluding comments

Mixing in VDT is a substantially more complicated process than passive scalar mixing, since VDT involves additional effects: (1) active scalar tightly coupled to the momentum equation, (2) non-Gaussian and skewed PDFs, and for our benchmark problem, (3) non-equilibrium turbulence.

We apply the idea that a model for the small scale mixing slaved to the large scales can be based on a suitable underlying single-point PDF. We essentially add a conservation law – the conservation of the PDF – as a constraint to the description of material mixing. This realisability requirement produces a model for the small scale ${\epsilon }_b$ that depends on large scale variables. At the level of principle this requirement is then a mathematically formal expression of the usual notion that in a turbulence where the production is in the large scales the small scale mixing is slaved to the large scales. At the level of practice we have replaced the single small-scale unknown ${\epsilon }_b$, about which we know very little, with two unknowns, the FPE drift and diffusion coefficients, about which we know considerably more.

In summary we highlight our findings as, we have

• Used moment equations derived from the Fokker–Planck equation together with moment equations from the Navier–Stokes equation to derive models for small-scale dissipative moments in terms of large-scale moments using a realisability principle.

• Found useful variable-density scalings for mix state and mix rate metrics.

• Developed and applied MC as a theoretical developmental tool for closure validation.

• Created an engineering closure for a second-order moment, b, applicable in moment approaches to non-equilibrium situations.

• Created a closure for a PDF method treatment of variable-density turbulence mixing.

Acknowledgements

We thank Daniel Livescu for providing extensive data sets from his direct numerical simulations and John Schwarzkopf for advise and input on the manuscript and general creative questions and insights.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

J. Bakosi http://orcid.org/0000-0002-0604-5555

Additional information

Funding

This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396, and supported by the Laboratory Directed Research and Development Program. Research presented in this article was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20150498ER. LANL Report LA-UR-18-23839.

Notes

1 For the problem to model mixing we require $\beta (1+\eta +\eta b)>-2{\kappa }^{\mathrm{\prime }}G$. Note that $\beta >0$ and $\kappa >0$.

2 The denominator does not deviate very substantially from 2.

Appendices

Appendix 1. Final model

The new b equation and model for the time-scale ratio is summarised as follows. (A1) $\begin{array}{rl}\frac{D}{Dt}b& =-C_{b2}\left(r,\tilde{Y},b,\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2},\frac{P}{\epsilon }\right)\frac{\epsilon }{k}b,\end{array}$(A1) (A2) $\begin{array}{rl}{C}_{b2}& =\frac{{\beta }_1}{2(1+\eta b)}\left[(1+\eta +\eta b)+2{\kappa }_{0}bG\right],\end{array}$(A2) where (A3) $\begin{array}{rl}\eta b& =\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}=\left[\frac{1-S{k}_{\beta }{b}^{1/2}+(2+\chi )b}{(1+b{)}^2}\right]b,\chi =\frac{(1-\overline{Y}{)}^3+{\overline{Y}}^3}{\overline{Y}(1-\overline{Y})},\end{array}$(A3) (A4) $\begin{array}{rl}\eta & =\frac{⟨{\rho }^2⟩}{b{\overline{\rho }}^2}=\left[\frac{1-S{k}_{\beta }{b}^{1/2}+(2+\chi )b}{(1+b{)}^2}\right],\overline{Y}=\tilde{Y}+\frac{1+r\tilde{Y}}{r}b.\end{array}$(A4) The skewness is taken from the beta PDF and is given by (A5) $S{k}_{\beta }=\frac{⟨y^{\mathrm{\prime }\mathrm{\prime }3}⟩}{(\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}{)}^{3/2}}=2\frac{(1-2\tilde{Y})(\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}{)}^{1/2}}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}+\tilde{Y}(1-\tilde{Y})}=2\frac{1-2\tilde{Y}}{\left(\tilde{Y}\right(1-\tilde{Y}){)}^{1/2}}\frac{(1-\mathrm{\Theta }{)}^{1/2}}{2-\mathrm{\Theta }}.$(A5) The above expression, given in [17], was tested in [25]. Coefficients and definitions are (A6) $\begin{array}{rl}{\beta }_1& =\frac{2}{1+\eta +\eta b}\left[{\beta }_{10}+{\beta }_2\frac{{\mathrm{\Theta }}^{\mathrm{\prime }}}{\sqrt{1+(P/\epsilon -1{)}^2{b}^{1/4}}}+{\beta }_3\frac{4{\mathrm{\Theta }}^{\mathrm{\prime }}(1-{\mathrm{\Theta }}^{\mathrm{\prime }})}{\sqrt{1+(P/\epsilon -1{)}^2{b}^{1/4}}}\right]\end{array}$(A6) (A7) $\begin{array}{rl}{\mathrm{\Theta }}^{\mathrm{\prime }}& =\frac{1}{2}\left(1+\frac{\eta }{1+\eta b}\right)\mathrm{\Theta },\mathrm{\Theta }=1-\frac{b}{b_{\mathrm{n}\mathrm{m}}},b_{\mathrm{n}\mathrm{m}}=\frac{r^2}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y})\end{array}$(A7) (A8) $\begin{array}{rl}{\beta }_{10}& =(\frac{1+b}{1+\eta b}{)}^a,{\beta }_2=[0.8\sim 1.0\left]\right(1+\eta b),{\beta }_3=(1+\eta b)\end{array}$(A8) (A9) $\begin{array}{rl}G& =\frac{-1}{(1+r\tilde{Y}{)}^2}\left[r^2\tilde{Y}(1-\tilde{Y})-(1+r)\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\right].\end{array}$(A9) The exponent a plays a small role at the initial condition at high A; in simulations both $a=1/2$ and a = 1 work well. The constant in ${\beta }_2$ sets the asymptotic passive scalar time-scale ratio. We typically take the constant in ${\beta }_2$ as 1 which produces an asymptotic time-scale ratio $C_{b2}=2$ for proof of concept; the value 0.8 compares more favorably to DNS data. The objects in Equation (A7(A7) $\begin{array}{rl}{\mathrm{\Theta }}^{\mathrm{\prime }}& =\frac{1}{2}\left(1+\frac{\eta }{1+\eta b}\right)\mathrm{\Theta },\mathrm{\Theta }=1-\frac{b}{b_{\mathrm{n}\mathrm{m}}},b_{\mathrm{n}\mathrm{m}}=\frac{r^2}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y})\end{array}$(A7) ) can be tweaked according to the users whim.

We have taken the function ${\kappa }_1={\kappa }_0{\beta }_1$ as our concern is with a second-order moments. If higher order moments are a concern then a deeper inquiry into the functional form of ${\kappa }_1$ is required. The coefficient ${\kappa }_0\approx 0.1$ does not appear to play a strong role for second-order moments.

In working with inhomogeneous flows it is better to express $1+r\tilde{Y}$ in terms of the mean density (A10) $\begin{array}{rl}1+r\tilde{Y}& =\frac{{\rho }_2}{\overline{\rho }},r=\frac{{\rho }_2}{{\rho }_1}-1=\frac{2A}{1-A},A=\frac{r}{2+r}\end{array}$(A10) (A11) $\begin{array}{rl}{b}_{\mathrm{n}\mathrm{m}}& =\frac{(\overline{\rho }-{\rho }_1)({\rho }_2-\overline{\rho })}{{\rho }_1{\rho }_2},\frac{b}{b_{\mathrm{n}\mathrm{m}}}=\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}_{\mathrm{n}\mathrm{m}}}=1-\mathrm{\Theta }.\end{array}$(A11) To close the equation at the second moment level we need to close $⟨{\rho }^3⟩$ in G. First let us see if this is necessary. (A12) $G=\frac{1}{(1+r\tilde{Y}{)}^2}\left\{-r^2\tilde{Y}(1-\tilde{Y})+\left[1+2r-(2+r)r\tilde{Y}\right]\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}+(1+r)\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}\right\}.$(A12) Do we need a closure? This question can be answered by answering the following two questions: 1) is $G_m\approx G$ and is $⟨{\rho }^3⟩$ small? (A13) $\begin{array}{rl}{G}_m& =\frac{1}{(1+r\tilde{Y}{)}^2}\left\{-r^2\tilde{Y}(1-\tilde{Y})+\left[1+2r-(2+r)r\tilde{Y}\right]\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\right\}.\end{array}$(A13) (A14) $\begin{array}{rl}\frac{⟨{\rho }^3⟩}{\overline{\rho }⟨{\rho }^2⟩}& \ll 1.\end{array}$(A14) If not, then we need to test the following two models (A15) $\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}=-\frac{r}{1+r}\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}[1-(2+r\left)\tilde{Y}\right].$(A15) Alternatively if one works with Reynolds variables one finds (A16) $\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}=-\frac{r}{1+r}\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\left[\frac{1-2\tilde{Y}-2b(1+r\tilde{Y})/r}{1+r\tilde{Y}}\right]$(A16) or go to the beta PDF. It appears that the first ‘no-mix’ model with the variance given by the mixed varince works best and the expression for G becomes (A17) $G=\frac{1}{(1+r\tilde{Y}{)}^2}\left[-r^2\tilde{Y}(1-\tilde{Y})+(1+r)\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\right].$(A17)

Appendix 2. No-mix results

The no mix results given in the text are easily obtained using the double-delta distribution (A18) $P_{\delta }\left(\hat{y}\right)={\alpha }_1\delta (1-\hat{y})+(1-{\alpha }_1)\delta \left(\hat{y}\right).$(A18) Taking the moments (A19) $\overline{Y}=\int \hat{y}{P}_{\delta }\left(\hat{y}\right)\mathrm{d}\hat{y}\overline{\rho }\tilde{Y}=\int \rho \left(\hat{y}\right)\hat{y}{P}_{\delta }\left(\hat{y}\right)\mathrm{d}\hat{y}=\int \frac{{\rho }_2}{1+r\hat{y}}\hat{y}{P}_{\delta }\left(\hat{y}\right)\mathrm{d}\hat{y}$(A19) one finds the value of ${\alpha }_1$. Using the usual $\overline{\rho }={\rho }_2/(1+r\tilde{Y})$ [17] one obtains the following relations between the Favre and Reynolds means (A20) ${\alpha }_1=\overline{Y}\stackrel{\mathrm{n}\mathrm{m}}{=}\tilde{Y}\frac{1+r}{1+r\tilde{Y}}1-{\alpha }_1=1-\overline{Y}\stackrel{\mathrm{n}\mathrm{m}}{=}\frac{1-\tilde{Y}}{1+r\tilde{Y}}.$(A20) Thus (in the no-molecular-mix limit) one has the following relations between the Favre and Reynolds means (A21) $\overline{Y}(1-\overline{Y})\stackrel{\mathrm{n}\mathrm{m}}{=}\frac{1+r}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y}),$(A21) which allows one to express no-mix identities in either Favre or Reynold variables. In general mixing situations $\overline{Y}=\tilde{Y}+b(1+r\tilde{Y})/r$. The mean specific volume is similarly found from the sample space integral to be (A22) $V=\int \frac{1+r\hat{y}}{{\rho }_2}{P}_{\delta }\left(\hat{y}\right)\mathrm{d}\hat{y}=\frac{1+r\overline{Y}}{{\rho }_2}.$(A22) The specific volume density variance is computed from (A23) $b_{\mathrm{n}\mathrm{m}}=\overline{\rho }{V}_{\mathrm{n}\mathrm{m}}-1=\frac{1+r\overline{Y}}{1+r\tilde{Y}}-1=\frac{r^2}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y})$(A23) after substituting in for $\overline{Y}$ in terms of $\tilde{Y}$ from the above no mix relations in Equation (A20(A20) ${\alpha }_1=\overline{Y}\stackrel{\mathrm{n}\mathrm{m}}{=}\tilde{Y}\frac{1+r}{1+r\tilde{Y}}1-{\alpha }_1=1-\overline{Y}\stackrel{\mathrm{n}\mathrm{m}}{=}\frac{1-\tilde{Y}}{1+r\tilde{Y}}.$(A20) ). The enumerator is maximum at $\tilde{Y}=1/2$ and $b_{\mathrm{n}\mathrm{m}}\le 1$ as the density ratio gets large. This is the same result one gets if one substitutes $\overline{\rho }$ in the no mix relation, ${\rho }_1{\rho }_2{b}_{\mathrm{n}\mathrm{m}}=(\overline{\rho }-{\rho }_1)({\rho }_2-\overline{\rho })$.

The no-mix values of other relevant second moments are then (A24) $\begin{array}{rl}{b}_{\mathrm{n}\mathrm{m}}& =\frac{(\overline{\rho }-{\rho }_1)({\rho }_2-\overline{\rho })}{{\rho }_1{\rho }_2}=\frac{r^2}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y})=\frac{r^2}{(1+r)}\overline{Y}(1-\overline{Y})\end{array}$(A24) (A25) $\begin{array}{rl}⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}& =(\overline{\rho }-{\rho }_1)({\rho }_2-\overline{\rho })={\overline{\rho }}^2\frac{r^2}{1+r}\tilde{Y}(1-\tilde{Y})={\rho }_2^2\frac{r^2}{(1+r{)}^2}\overline{Y}(1-\overline{Y})\end{array}$(A25) (A26) $\begin{array}{rl}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}_{\mathrm{n}\mathrm{m}}& =\tilde{Y}(1-\tilde{Y})=\frac{(\overline{\rho }-{\rho }_1)({\rho }_2-\overline{\rho })}{r^2{\overline{\rho }}^2}(1+r)\end{array}$(A26) (A27) $\begin{array}{rl}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }3}}}_{\mathrm{n}\mathrm{m}}& =\tilde{Y}(1-\tilde{Y})[1-2\tilde{Y}]\end{array}$(A27) (A28) $\begin{array}{rl}S{k}_{\mathrm{n}\mathrm{m}}& =\frac{1-2\overline{Y}}{\left(\overline{Y}\right(1-\overline{Y}){)}^{1/2}}\end{array}$(A28) (A29) $\begin{array}{rl}⟨v^2{⟩}_{\mathrm{n}\mathrm{m}}& ={\left(\frac{r}{{\rho }_2}\right)}^2\overline{Y}(1-\overline{Y})={\left(\frac{r}{{\rho }_2}\right)}^2\frac{1+r}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y})\end{array}$(A29) (A30) $\begin{array}{rl}⟨y^2{⟩}_{\mathrm{n}\mathrm{m}}& =\overline{Y}(1-\overline{Y})=\frac{1+r}{(1+r\tilde{Y}{)}^2}\tilde{Y}(1-\tilde{Y})\end{array}$(A30) (A31) $\begin{array}{rl}⟨{\rho }^{\ast }{y}^2{⟩}_{\mathrm{n}\mathrm{m}}& ={\rho }_2\overline{Y}(1-\overline{Y})\frac{1+r\overline{Y}}{1+r}\end{array}$(A31) (A32) $\begin{array}{rl}⟨y^3{⟩}_{\mathrm{n}\mathrm{m}}& =\overline{Y}(1-\overline{Y})[1-2\overline{Y}]\end{array}$(A32) (A33) $\begin{array}{rl}⟨{\rho }^3{⟩}_{\mathrm{n}\mathrm{m}}& =-(\overline{\rho }-{\rho }_1)({\rho }_2-\overline{\rho })[{\rho }_1+{\rho }_2-2\overline{\rho }].\end{array}$(A33) The results satisfy the Cauchy-Schwarz inequality, e.g. $⟨{\rho }^2⟩⟨v^2⟩\le ⟨\rho v{⟩}^2$.

The above reflect the maximum values of the second moments and are useful for scaling and ensuring boundedness requirements. In the no-mix limit the density intensity normalised by r is writtenas (A34) $\begin{array}{rl}\frac{⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}}{r^2{\overline{\rho }}^2}& =\frac{\tilde{Y}(1-\tilde{Y})}{1+r}=\overline{Y}(1-\overline{Y})\frac{(1+r\tilde{Y}{)}^2}{(1+r{)}^2}\end{array}$(A34) (A35) $\begin{array}{rl}\frac{⟨{\rho }^3{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^3}& =-\frac{r}{1+r}\frac{⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^2}[1-(2+r\left)\tilde{Y}\right].\end{array}$(A35) Alternatively, if one works with Reynolds variables one finds (A36) $\frac{⟨{\rho }^3{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^3}=-\frac{r}{1+r}\frac{⟨{\rho }^2{⟩}_{\mathrm{n}\mathrm{m}}}{{\overline{\rho }}^2}\left[\frac{1-2\overline{Y}}{1+r\tilde{Y}}\right],$(A36) which suggests the closure model under mixing conditions for the cubic moment as (A37) $\frac{⟨{\rho }^3⟩}{{\overline{\rho }}^3}=-\frac{r}{1+r}\frac{⟨{\rho }^2⟩}{{\overline{\rho }}^2}\left[\frac{1-2\tilde{Y}-2b(1+r\tilde{Y})/r}{1+r\tilde{Y}}\right].$(A37) The Atwood number and r are related as (A38) $A=\frac{r}{2+r}r=\frac{2A}{1-A},$(A38) Our DNS data set is at A = 0.05, 0.5, 0.75, corresponding to $r=0.1,2.0,6.0$, respectively. Following the same methods the fourth-order moment appearing in the constitutive relation for the density variance can be found to be, in the no-mix limit, as (A39) $\begin{array}{rl}⟨{\rho }^2{v}^2{⟩}_{\mathrm{n}\mathrm{m}}& =\int (\hat{\rho }-\overline{\rho }{)}^2(\hat{v}-V{)}^2\mathrm{d}P=\frac{r^4}{(1+r{)}^2}\overline{Y}(1-\overline{Y})\left[\right(1-\overline{Y}{)}^3+{\overline{Y}}^3]\end{array}$(A39) (A40) $\begin{array}{rl}& =\frac{r^4}{(1+r{)}^2}\overline{Y}(1-\overline{Y})[1-3\overline{Y}(1-\overline{Y}\left)\right].\end{array}$(A40)

Appendix 3. Select $\beta$-PDF facts

Facts regarding the beta PDF. The problem solution requires dealing with both Favre and Reynolds averages. The Favre-, and Reynolds-averaged PDFs are related by (A41) $\tilde{P}(\hat{y},x)=\frac{{\rho }^{\ast }\left(\hat{y}\right)}{\overline{\rho }\left(x\right)}P(\hat{y},x).$(A41) The Favre PDF is a two-parameter PDFs such that (A42) $\tilde{P}(\hat{y},x)=\tilde{P}(\hat{y},x;\tilde{Y},\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}})\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}=\frac{⟨{\rho }^{\ast }{y}^{\mathrm{\prime }\mathrm{\prime }2}⟩}{\overline{\rho }}.$(A42) One PDF in that class is the Favre β-PDF: (A43) $\begin{array}{rl}\tilde{P}(\hat{y};\alpha ,\beta )& ={\hat{y}}^{\alpha -1}(1-\hat{y}{)}^{\beta -1}{B}^{-1}{\int }_0^1{\hat{y}}^{\alpha -1}(1-\hat{y}{)}^{\beta -1}\mathrm{d}\hat{y}=B\end{array}$(A43) (A44) $\begin{array}{rl}\alpha & =\tilde{Y}\frac{\tilde{Y}(1-\tilde{Y})}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}\left[1-\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{\tilde{Y}(1-\tilde{Y})}\right]=\tilde{Y}\frac{\theta }{1-\theta }\end{array}$(A44) (A45) $\begin{array}{rl}\beta & =(1-\tilde{Y})\frac{\tilde{Y}(1-\tilde{Y})}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}\left[1-\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{\tilde{Y}(1-\tilde{Y})}\right]=(1-\tilde{Y})\frac{\theta }{1-\theta }\end{array}$(A45) (A46) $\begin{array}{rl}\theta & =1-\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{\tilde{Y}(1-\tilde{Y})}\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}=⟨{\rho }^{\ast }{y}^{\mathrm{\prime }\mathrm{\prime }2}⟩/\overline{\rho }\end{array}$(A46) (A47) $\begin{array}{rl}\tilde{P}(\hat{y},x;\alpha ,\beta )& =\tilde{P}(\hat{y},x;\tilde{Y},\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}).\end{array}$(A47) The quantity θ is sometimes called (for better or worse) the mixedness. It is more useful to think of θ as a measure of the relative scalar variance: the variance is maximum when $\theta =0$ which corresponds to the no-mix limit. For mixed fluids, i.e. low variance, $\theta \approx 1$.

The results above are found by applying the beta PDF identities (A48) $\tilde{Y}=\frac{\alpha }{\alpha +\beta }\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}=\frac{\alpha \beta }{(\alpha +\beta {)}^2(\alpha +\beta +1)}.$(A48) We note that (A49) $\begin{array}{rl}b& =\frac{r^2}{(1+r\tilde{Y}{)}^2}\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}={\left(\frac{r\overline{\rho }}{{\rho }_2}\right)}^2\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}.\end{array}$(A49) (A50) $\begin{array}{rl}\frac{b}{b_{\mathrm{n}\mathrm{m}}}& =\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}_{\mathrm{n}\mathrm{m}}}\end{array}$(A50) The centred Favre-skewness is givenby (A51) $\frac{⟨y^{\mathrm{\prime }\mathrm{\prime }3}⟩}{(\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}{)}^{3/2}}=\frac{2(\beta -\alpha )(\alpha +\beta +1{)}^{1/2}}{(\alpha +\beta +2)(\alpha \beta {)}^{1/2}}=2\frac{(1-2\tilde{Y})(\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}{)}^{1/2}}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}+\tilde{Y}(1-\tilde{Y})}=2\frac{1-2\tilde{Y}}{\left(\tilde{Y}\right(1-\tilde{Y}){)}^{1/2}}\frac{(1-\theta {)}^{1/2}}{2-\theta }.$(A51) In the last equation above relations (A44(A44) $\begin{array}{rl}\alpha & =\tilde{Y}\frac{\tilde{Y}(1-\tilde{Y})}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}\left[1-\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{\tilde{Y}(1-\tilde{Y})}\right]=\tilde{Y}\frac{\theta }{1-\theta }\end{array}$(A44) ), (A45(A45) $\begin{array}{rl}\beta & =(1-\tilde{Y})\frac{\tilde{Y}(1-\tilde{Y})}{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}\left[1-\frac{\widetilde{y^{\mathrm{\prime }\mathrm{\prime }2}}}{\tilde{Y}(1-\tilde{Y})}\right]=(1-\tilde{Y})\frac{\theta }{1-\theta }\end{array}$(A45) ) have been used.