Self-similar Reynolds-averaged mechanical–scalar turbulence models for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz instability-induced mixing in the small Atwood number limit

Schilling, Oleg

doi:10.1063/5.0055193

Analytical self-similar solutions to two-, three-, and four-equation Reynolds-averaged mechanical–scalar turbulence models describing incompressible turbulent Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz instability-induced mixing in planar geometry are derived in the small Atwood number (Boussinesq) limit. The models are based on the turbulent kinetic energy K and its dissipation rate ε, together with the scalar (heavy-fluid mass fraction) variance S and its dissipation rate $χ$ modeled either differentially or algebraically. The models allow for a simultaneous description of mechanical and scalar mixing, i.e., mixing layer growth and molecular mixing, respectively. Mixing layer growth parameters and other physical observables relevant to each instability are obtained explicitly as functions of the model coefficients. The turbulent fields are also expressed in terms of the model coefficients, with their temporal power-law scalings obtained by requiring that the self-similar equations are explicitly time-independent. The model calibration methodology is described and discussed. Expressions for a subset of the various physical observables are used to calibrate each of the two-, three-, and four-equation models, such that the self-similar solutions are consistent with experimental and numerical simulation data corresponding to these values of the observables and to specific canonical Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz turbulent flows. A calibrated four-equation model is then used to reconstruct the mean and turbulent fields, and late-time turbulent equation budgets for each instability-induced flow across the mixing layer. The reference solutions derived here can provide systematic calibrations and better understanding of mechanical–scalar turbulence models and their predictions for instability-induced turbulent mixing in the very large Reynolds number limit.

I. INTRODUCTION

Advances in computing capability over the past decades have enabled increasingly detailed direct numerical simulation (DNS), large-eddy simulation (LES), and implicit large-eddy simulation (ILES) studies of turbulent mixing layers arising from classical hydrodynamic Rayleigh–Taylor (RT), Richtmyer–Meshkov (RM), and Kelvin–Helmholtz (KH) instabilities. These advances have engendered better fundamental understanding of the linear, nonlinear, transitional, and turbulent flow regimes arising from these interfacial instabilities.^1,2 However, these flows are complex and depend on a large number of dimensionless flow parameters (e.g., Atwood, Mach, Reynolds, Schmidt, and Prandtl numbers), fluid properties (e.g., molecular transport coefficients of gases and liquids), and initial conditions (e.g., modal amplitudes and wavelengths of interfacial perturbations). The flow complexity of such turbulent mixing layers is further exacerbated in the high-energy-density regime in which plasma physics and radiation hydrodynamics must also be considered.³ Furthermore, many of these turbulent flows are found in astrophysical, geophysical, and engineering applications in which the range of temporal and spatial scales and other physics preclude experiments or simulations that reach the conditions required to achieve the transition to a fully turbulent state. As a result, there remains a need to develop, verify, and validate computationally less expensive, reduced-order models that adequately represent the effects of turbulence on the mean hydrodynamic flows arising in these applications. Reynolds-averaged turbulence models are thus used to render the simulation of large Reynolds number turbulent mixing computationally tractable in such applications.

A large variety of Reynolds-averaged turbulence models based on different choices of mechanical turbulent fields (e.g., turbulent kinetic energy dissipation rate or lengthscale) and of varying degrees of complexity ranging from two-equation to Reynolds stress models have been developed to model RT, RM, and KH mixing (see Ref. 2). These models contain a large number of dimensionless coefficients that are broadly specified using two different approaches: (1) most models are calibrated by iteratively modifying a subset of the coefficients until the models adequately predict some chosen set of data, and (2) a much smaller number of models are calibrated using analytic solutions to the Reynolds-averaged equations in the self-similar approximation to constrain the model coefficients in order to predict key observables, such as the mixing layer growth parameters and exponents. The models so calibrated are then applied to other, more diverse flows to quantitatively evaluate their overall predictive capability. In the second approach, the emphasis is on determining the values of the model coefficients that give specified growth parameters and other key quantities that have been measured experimentally or obtained numerically. Complete expressions for the turbulent fields, terms in the turbulent transport equations, and other derived self-similar quantities are typically not given. Moreover, previous self-similarity analyses of turbulence models for the applications of interest here primarily considered mechanical turbulence, thereby precluding them from predicting scalar mixing properties, such as molecular mixing parameters, scalar production-to-dissipation/destruction ratios, and mechanical-to-scalar timescale ratios. A recent review of simulation and modeling of RT instability and mixing experiments, including the application of Reynolds-averaged models, is available.⁴

Self-similar solutions to the two-equation K–ε model were derived^5–7 for RT, RM, and KH mixing. Most effort has considered calibrating the coefficients in the buoyancy/shock production terms in the K and ε equations for RT and RM mixing, which are the most important, least understood, and least universal model coefficients.⁸ Self-similarity analysis was used to examine the growth of a small Atwood number RT mixing layer developing in a water channel experiment,⁶ solving a mean temperature equation using coefficients similar to those used in the standard K–ε model:^9,10 the value of the buoyancy production coefficient $C_{ε 0}$ was determined from similarity and the experimental value of the RT growth parameter α. It was shown that a K–ε model with calibrated buoyancy production terms could reproduce experimentally measured reshocked RM mixing layer growth.¹¹ This model was later extended to a Reynolds stress model to capture local anisotropy effects.¹² Transport equations for the mass flux velocity and density variance were included to better capture buoyancy effects and turbulent mixing. The conditions needed for and implications of self-similarity of RT, RM, and KH mixing modeling in the context of single- and two-fluid K–ε models were extensively discussed.⁸ Model calibration strategies were also presented in Ref. 8, in which the self-similar equations were integrated over the mixing layer to reduce the equations to easily solvable ordinary differential equations (“zero-dimensional models”) in time to obtain constraints between model coefficients and physical observables. These studies only considered mechanical turbulence.

A K–L formulation, where L is a turbulent lengthscale—the Besnard–Harlow–Rauenzahn (BHR) model—was developed to model variable-density and buoyancy-driven turbulence,¹³ where gradient-diffusion and similarity closures were used to model the turbulent transport equations. Transport equations for the mass flux velocity a_i and negative density fluctuation–specific volume fluctuation correlation b were included in the model. A single-fluid K–L model was also developed,¹⁴ in which the buoyancy production term in the K equation was modeled via a relation with buoyancy–drag models. Self-similar solutions to the model equations for RT and RM mixing in the Boussinesq approximation were used to calibrate the model coefficients and demonstrate the ability of the model to reproduce experimentally measured RT and RM mixing layer widths and scaling parameters. A variant of this model that included shear production terms was applied to low energy density¹⁵ and high energy density¹⁶ turbulent shear layers. A self-similar analysis of the K–L–a version of the BHR model was later performed for small Atwood number RT flow, and the model was applied to predict RT, RM, and KH mixing quantities (with many model coefficient values adjusted for each instability case).¹⁷

Self-similar solutions to an alternative formulation of the K–L–a model were presented¹⁸ and later extended to Kelvin–Helmholtz flows using a second turbulent lengthscale.¹⁹ The BHR and BHR-2²⁰ (i.e., K–L–a–b) models were subsequently extended to include transport equations for the Reynolds stress tensor²¹ and for a second lengthscale.²² Conditions for the attainment of self-similarity of turbulent fields and statistics were introduced and investigated using artificial fluid large-eddy simulation of small Atwood RT mixing.²³ The mean heavy fluid mass fraction and self-normalized K, a, and b profiles from the simulation and from the K–L–a and BHR-2 models were compared. The predictions of several related lengthscale-based models applied to RT mixing were compared.^24,25 Most recently, the Dimonte–Tipton model¹⁴ was further improved for RT, reshocked RM, and KH mixing.^26–28 In particular, these studies also proposed systematic self-similar methodologies to obtain a unified set of model coefficients for different mixing cases and made several modifications to improve the accuracy and robustness of the models. The K–L–a models^17,18 using an algebraic two-fluid expression for b describe immiscible mixing and are not appropriate for fluids that molecularly mix.

More recently, an advanced compressible, multicomponent K–ε model was developed and shown to well-predict several reshocked RM instability experiments.^29,30 Three- and four-equation scalar generalizations of the K–ε model were subsequently proposed and calibrated a priori for RT mixing.³¹ Note that a mechanical turbulence-based K– $ε$ or K–L model can only predict large-scale quantities, such as the bubble and spike front widths and their growth rates: additional equations describing transport and diffusion of a scalar are needed to predict higher-order mixing statistics, such as the degree of molecular mixing or mechanical-to-scalar timescale ratio. The four-equation model framework discussed here allows for a simultaneous description of large-scale (K and S) and small-scale (ε and $χ$ ⁠) fields, while a two-equation K–L or three-equation K–L–a model only describes large-scale fields.^14,18

The self-similar model analysis, calibration, and model applications presented here extend and augment previous studies as follows. First, the K and ε equations describing RT, RM, and KH mixing are augmented by S (scalar variance) and χ (scalar variance dissipation rate) equations that model passive scalar mixing. Thus, three- and four-equation models that treat mechanical and scalar turbulence and can model scalar statistics are considered. Expressions are derived for constant self-similar parameters not considered previously, such as molecular mixing parameters, production-to-dissipation/destruction ratios, and the mechanical-to-scalar timescale ratios. Inequalities that must be satisfied by the model coefficients that follow from positivity constraints are noted to avoid unphysical coefficient choices. Explicit space–time-dependent expressions for all turbulent fields are presented in terms of the model coefficients.

Second, new calibrations of the two- and three-equation K–ε models based on analytic solutions to the model coefficients as functions of the growth parameters and exponents for RT, RM, and KH mixing treated simultaneously are presented. Several calibrations are described and discussed in detail, including the methodologies for calibrating four-equation models. The predicted values of the self-similar quantities are summarized for each model and discussed. Third, the analytic expressions for the key mean and turbulent fields are evolved in space and time numerically as applied to particular small Atwood number cases for each instability to illustrate the rate of spreading in space and growth or decay in time of the fields, which has not been shown and discussed in previous studies. Furthermore, the analytic solutions are used to reconstruct the production, dissipation/destruction, and turbulent diffusion terms for each turbulent transport equation and shown at the latest evolution time to illustrate the relative importance of these contributions to the equation budgets, which has not been shown and discussed in previous studies.

Finally, a clear, concise, and complete discussion of the assumptions used to formulate the self-similar approximation, as well as of the steps required to derive the solutions, is provided. The detailed derivation of the self-similar solutions is formulated compactly and in a unified manner for all three instabilities. The steps in the derivation are clearly described to allow reproducing the results and to apply the same procedure to other fluid dynamics problems. The procedure for determining the power-laws of each mean and turbulent field is presented in detail, rather than assumed using dimensional analysis. Therefore, the emphasis of the present analysis is the physics of the self-similar solutions rather than numerical applications of calibrated Reynolds-averaged models as presented in previous work.

This paper is organized as follows. The formulation of the two-, three-, and four-equation turbulence models analyzed here is presented in Sec. II. The analytical self-similar solutions (including fields, growth parameters, and other self-similar quantities) of the models for RT, RM, and KH mixing are given in Sec. III, IV, and V, respectively. To illustrate the utility of the solutions and discuss their behavior, the two-, three-, and four-equation model coefficients are calibrated in Sec. VI and a recommended four-equation model is applied to each instability case in Sec. VII. A summary of the results and discussion are given in Sec. VIII. The self-similar solutions of the model equations for RT, RM, and KH mixing, for decaying homogeneous isotropic turbulence, and for shockless rapid compression are derived in detail in Appendices A–C, respectively.

II. MECHANICAL–SCALAR TURBULENCE MODEL FORMULATIONS AND SIMILARITY ANALYSIS

The three- and four-equation Reynolds-averaged Navier–Stokes (RANS) models evaluated a priori³¹ using direct numerical simulation (DNS) data^32–34 modeling a very small Atwood number, Schmidt number Sc $= 7$ hot/cold water channel Rayleigh–Taylor mixing experiment³⁵ are briefly summarized here. These models are appropriate for miscible binary fluid mixing and are closely related to the models that have been used for reacting turbulent flows. The scalar fields are passive in these models as the mechanical turbulence equations do not depend on the scalar variance or its dissipation rate S or χ.

A. Transport equations

For a general incompressible field with Reynolds decomposition $ϕ_{α} = {\bar{ϕ}}_{α} + ϕ_{α}^{'}$ and compressible field with Favre decomposition $ϕ_{α} = {\tilde{ϕ}}_{α} + ϕ_{α}^{″}$ ⁠, the Reynolds- and Favre-averaged fields are ${\bar{ϕ}}_{α}$ and ${\tilde{ϕ}}_{α} = {\bar{ρ ϕ}}_{α} / \bar{ρ}$ ⁠, respectively, and $ϕ_{α}^{'}$ and $ϕ_{α}^{″}$ are the Reynolds and Favre fluctuating fields, respectively. The mean scalar, mean momentum, mean internal energy, turbulent kinetic energy, turbulent kinetic energy dissipation rate, scalar variance, and scalar variance dissipation rate transport equations are (four-equation model and assuming the summation convention for repeated vector indices)

\bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) \tilde{s} = \frac{\partial}{\partial x_{j}} [(\frac{\bar{μ}}{Sc} + \frac{μ_{t}}{σ_{s}}) \frac{\partial \tilde{s}}{\partial x_{j}}],

(1)

\bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) {\tilde{v}}_{i} = \bar{ρ} g_{i} - \frac{\partial \bar{p}}{\partial x_{i}} + \frac{\partial {\bar{σ}}_{i j}}{\partial x_{j}} - \frac{\partial τ_{i j}}{\partial x_{j}},

(2)

\bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) \tilde{U} = - a_{j} \frac{\partial \bar{p}}{\partial x_{j}} - \bar{p} \frac{\partial {\tilde{v}}_{j}}{\partial x_{j}} + ({\tilde{σ}}_{i j} + \bar{σ_{i j}^{″}}) \frac{\partial {\tilde{v}}_{i}}{\partial x_{j}} - {\tilde{σ}}_{i j} \frac{\partial a_{i}}{\partial x_{j}} + \bar{ρ} ε + \frac{\partial}{\partial x_{j}} (\frac{μ_{t}}{σ_{U}} \frac{\partial \tilde{h}}{\partial x_{j}}),

(3)

\begin{matrix} \bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) K = a_{j} (\frac{\partial \bar{p}}{\partial x_{j}} - \frac{\partial {\bar{σ}}_{i j}}{\partial x_{i}}) - τ_{i j} \frac{\partial {\tilde{v}}_{i}}{\partial x_{j}} \\ - \bar{ρ} ε + \frac{\partial}{\partial x_{j}} [(\bar{μ} + \frac{μ_{t}}{σ_{K}}) \frac{\partial K}{\partial x_{j}}], \end{matrix}

(4)

\begin{matrix} \bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) ε = C_{ε 0} \frac{ε}{K} a_{j} (\frac{\partial \bar{p}}{\partial x_{j}} - \frac{\partial {\bar{σ}}_{i j}}{\partial x_{i}}) \\ - \frac{2}{3} C_{ε 3} \bar{ρ} ε \frac{\partial {\tilde{v}}_{j}}{\partial x_{j}} - C_{ε 1} \frac{ε}{K} τ_{i j}^{dev} \frac{\partial {\tilde{v}}_{i}}{\partial x_{j}} \\ - C_{ε 2} \frac{\bar{ρ} ε^{2}}{K} + \frac{\partial}{\partial x_{j}} [(\bar{μ} + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}], \end{matrix}

(5)

\begin{matrix} \bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) S = \frac{2 μ_{t}}{σ_{s}} {(\frac{\partial \tilde{s}}{\partial x_{j}})}^{2} - 2 \bar{ρ} χ + \frac{\partial}{\partial x_{j}} [\bar{ρ} (\bar{D} + \frac{ν_{t}}{σ_{S}}) \frac{\partial S}{\partial x_{j}}], \end{matrix}

(6)

\begin{matrix} \bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) χ = C_{χ 0} \frac{μ_{t}}{\sqrt{Sc}} \frac{ε}{K} {(\frac{\partial \tilde{s}}{\partial x_{j}})}^{2} - \bar{ρ} (C_{χ 2} \frac{χ}{S} - C_{χ 3} \frac{ε}{K}) χ \\ + \frac{\partial}{\partial x_{j}} [\bar{ρ} (\bar{D} + \frac{ν_{t}}{σ_{χ}}) \frac{\partial χ}{\partial x_{j}}], \end{matrix}

(7)

respectively, where ${\tilde{v}}_{i}$ is the mean velocity vector, $g_{i} = - g_{0} δ_{i z}$ is the acceleration vector, $\bar{ρ}$ is the mean density, $\bar{p}$ is the mean pressure, $τ_{i j} = \bar{ρ v_{i}^{″} v_{j}^{″}} = τ_{i j}^{iso} + τ_{i j}^{dev} = (2 / 3) \bar{ρ} K δ_{i j} + τ_{i j}^{dev}$ is the Reynolds stress tensor (decomposed into its isotropic plus deviatoric part), ${\bar{σ}}_{i j} = 2 \bar{μ} ({\tilde{S}}_{i j} - δ_{i j} {\tilde{S}}_{k}^{k} / 3)$ is the mean viscous stress tensor, ${\tilde{S}}_{i j} = (1 / 2) (\partial {\tilde{v}}_{i} / \partial x_{j} + \partial {\tilde{v}}_{j} / \partial x_{i})$ is the mean strain-rate tensor, $\bar{ν} = \bar{μ} / \bar{ρ}$ is the mean mixture kinematic viscosity, $\bar{D}$ is the mean mixture scalar diffusivity, Sc $= \bar{ν} / \bar{D}$ is the mixture Schmidt number, $a_{i} = - \bar{v_{i}^{″}} = \bar{ρ^{'} v_{i}^{'}} / \bar{ρ}$ is the mass flux velocity, $\tilde{U}$ is the mean internal energy, $\tilde{h} = \tilde{U} + \bar{p} / \bar{ρ}$ is the mean enthalpy, $K = \tilde{v_{i}^{″} v_{i}^{″}} / 2$ is the turbulent kinetic energy, $ε = \bar{σ_{i j} \partial v_{i}^{″} / \partial x_{j}} / \bar{ρ}$ is the turbulent kinetic energy dissipation rate, $S = \tilde{s^{″ 2}}$ is the scalar variance, $χ = \bar{D} \bar{ρ {(\partial s^{″} / \partial x_{j})}^{2}} / \bar{ρ}$ is the scalar variance dissipation rate, $C_{ε 0}$ – $C_{ε 3}$ and $C_{χ 0}$ – $C_{χ 3}$ are positive dimensionless coefficients, and σ_s, σ_U, σ_K, $σ_{ε}$ ⁠, σ_S, and $σ_{χ}$ are positive dimensionless turbulent Schmidt numbers. Here, it is assumed that there is no distinction between the coefficients ${σ_{K}^{*}, σ_{ε}^{*}}$ and ${σ_{K}, σ_{ε}}$ discussed in Ref. 31. Also, Eq. (7) analyzed here differs from that in Ref. 31, which includes $\sqrt{R e_{t}}$ factors multiplying the $C_{χ 2}$ and $C_{χ 3}$ coefficients, where $R e_{t} = ν_{t} / (C_{μ} \bar{ν})$ is the turbulent Reynolds number and

ν_{t} = \frac{μ_{t}}{\bar{ρ}} = C_{μ} \frac{K^{2}}{ε}

(8)

is the turbulent viscosity (μ_t is the dynamic turbulent viscosity and $C_{μ}$ is a positive dimensionless coefficient): as mentioned in Appendix A, such $\sqrt{R e_{t}}$ factors would break self-similarity.

Alternatively, instead of solving Eq. (7), an algebraic closure

χ = C_{χ} \frac{ε}{K} S,

(9)

where $C_{χ}$ is a positive dimensionless coefficient, can be used for the scalar variance dissipation rate in Eq. (6) (three-equation model). A concentration variance equation similar to Eq. (6) was first introduced in the modeling of a round, turbulent free jet and closed using Eq. (9).³⁶ Such three-equation turbulence models based on the K, ε, and S equations have been adopted in most studies of turbulent reacting flows (see Refs. 37–42 for further discussion and details). The mean scalar field can be the mean heavy fluid mass fraction $\tilde{s} = {\tilde{m}}_{1}$ (with ${\tilde{m}}_{2} = 1 - {\tilde{m}}_{1}$ the mean light fluid mass fraction), the mean concentration $\tilde{s} = \tilde{c}$ ⁠, or a mean progress variable $\tilde{s} = \tilde{ζ}$ ⁠. These fields correspond to the heavy fluid mass fraction variance $S = \tilde{m_{1}^{″ 2}}$ ⁠, the density variance $S = \bar{ρ^{' 2}}$ ⁠, the concentration variance $S = \tilde{c^{″ 2}}$ ⁠, or a progress variable variance $S = \tilde{ζ^{″ 2}}$ ⁠. Rather than using Eq. (9), a scalar variance dissipation rate equation similar to Eq. (7) has been used to improve the fidelity of the modeling of turbulent premixed and partially premixed combustion and flames.^43–45

The algebraic closure of the mass flux velocity is (for ideal gases with mixture ratio of specific heats $\bar{γ}$ ⁠)

a_{i} = - \frac{ν_{t}}{σ_{ρ} \bar{ρ}} (\frac{\partial \bar{ρ}}{\partial x_{i}} - \frac{\bar{ρ}}{\bar{γ} \bar{p}} \frac{\partial \bar{p}}{\partial x_{i}}),

(10)

the buoyancy-modified Reynolds stress tensor is³⁴

\begin{matrix} τ_{i j} = \frac{2}{3} \bar{ρ} K δ_{i j} - 2 μ_{t} ({\tilde{S}}_{i j} - \frac{δ_{i j}}{3} \frac{\partial {\tilde{v}}_{k}}{\partial x_{k}}) \\ - C_{A}^{ijkl} ν_{t} \frac{K}{\bar{ρ} ε} (\frac{\partial \bar{ρ}}{\partial x_{k}} \frac{\partial \bar{p}}{\partial x_{l}} + \frac{\partial \bar{ρ}}{\partial x_{l}} \frac{\partial \bar{p}}{\partial x_{k}} - \frac{2}{3} δ_{k l} \frac{\partial \bar{ρ}}{\partial x_{m}} \frac{\partial \bar{p}}{\partial x_{m}}), \end{matrix}

(11)

where $C_{A}^{ijkl}$ are dimensionless tensor coefficients. The term in the last parentheses in Eq. (11) is often called the baroclinic tensor and vanishes for a constant density flow.

The turbulent lengthscale is

L = C_{ε} \frac{K^{3 / 2}}{ε},

(12)

where $C_{ε}$ is a dimensionless coefficient (often referred to as the dimensionless turbulent kinetic energy dissipation rate), the mechanical turbulent timescale is (ω is the turbulent frequency)

τ_{m} = \frac{1}{ω} = \frac{K}{ε},

(13)

and the scalar turbulent timescale is

τ_{s} = \frac{S}{2 χ},

(14)

which is $τ_{s} = τ_{m} / (2 C_{χ})$ for a three-equation model using Eq. (9). For the three-equation model, the mechanical-to-scalar timescale ratio is

R = \frac{τ_{m}}{τ_{s}} = 2 C_{χ}

(15)

(see Ref. 31).

The mean internal energy equation follows that in Ref. 11, and thermal conduction is neglected. Using the identity $- \bar{p} \partial \bar{v_{j}^{″}} / \partial x_{j} = \bar{v_{j}^{″}} \partial \bar{p} / \partial x_{j} - (\partial / \partial x_{j}) (\bar{p} \bar{v_{j}^{″}})$ and the closure $\bar{p} \bar{v_{j}^{″}} = - (μ_{t} / σ_{U}) (\partial / \partial x_{j}) (\bar{p} / \bar{ρ})$ arising from Eq. (10) with γ = 1 and the assumption $σ_{ρ} = σ_{U}$ in this term (which will be verified later using the similarity analysis), the exact equation

\begin{matrix} \bar{ρ} (\frac{\partial}{\partial t} + {\tilde{v}}_{j} \frac{\partial}{\partial x_{j}}) \tilde{U} = - \bar{p} \frac{\partial \bar{v_{j}^{″}}}{\partial x_{j}} - \bar{p} \frac{\partial {\tilde{v}}_{j}}{\partial x_{j}} + ({\tilde{σ}}_{i j} + \bar{σ_{i j}^{″}}) \frac{\partial {\tilde{v}}_{i}}{\partial x_{j}} \\ + {\tilde{σ}}_{i j} \frac{\partial \bar{v_{i}^{″}}}{\partial x_{j}} + \bar{ρ} ε + \frac{\partial}{\partial x_{j}} (\frac{μ_{t}}{σ_{U}} \frac{\partial \tilde{U}}{\partial x_{j}}) \end{matrix}

(16)

reduces to Eq. (3).

The fifteen dimensionless model coefficients are σ_s, $σ_{U}$ ⁠, σ_K, $σ_{ε}$ ⁠, σ_S, $σ_{χ}, C_{A}^{ijkl}, σ_{ρ}, C_{μ}, C_{ε 0}, C_{ε 1}, C_{ε 2}, C_{χ 0}, C_{χ 2}$ ⁠, and $C_{χ 3}$ ⁠. The pressure–dilatation correlation is neglected in the K and ε equations. Although the shear production term $- τ_{i j} \partial {\tilde{v}}_{i} / \partial x_{j}$ in Eqs. (4) and (5) was shown to be very small in small Atwood number RT flow,³⁴ they are included above for completeness and are required for KH mixing. The gradient of the deviatoric part of the Reynolds stress component $τ_{z z}^{dev}$ was shown to be nonnegligible in very small Atwood number RT mixing.³⁴

Equations (1)–(7) allow the model to be applied to flows with different Schmidt numbers. It is expected that solving a transport equation for χ will give better predictions than using an algebraic model in which the mechanical-to-scalar timescale ratio (i.e., $C_{χ}$ ⁠) is a constant independent of the model coefficients. It should be noted that mixing in RT and RM flows cannot be predicted solely using a mechanical turbulence model (i.e., only based on K, ε or L): transport equations describing the evolution of a scalar mixing progress variable (and perhaps its dissipation rate) must also be included. This is also true for higher-order (i.e., second-order or second-moment) models.

B. Self-similarity assumptions

Consider the one-dimensional form of the model equations (1)–(6) and either the algebraic model (9) or the transport equation (7) for which the mixing layer is within $z \in [h_{s} (t), h_{b} (t)]$ and $h_{s} (t) < 0$ ⁠. The mixing layer width is $h (t) = h_{b} (t) - h_{s} (t) \approx 2 h_{b} (t)$ for very small Atwood number At $= (ρ_{1} - ρ_{2}) / (ρ_{1} + ρ_{2}) ≪ 1$ (i.e., the bubble and spike front growth is nearly symmetric).

To analytically derive the self-similar solutions in the Boussinesq approximation, assume

small density difference $ρ_{1} - ρ_{2} ↓ 0$ (or At $↓ 0$ ⁠);
no turbulence outside the mixing layer [ $K = ε = S = χ = 0$ for $| z | > h (t) / ξ$ ⁠, where the scale parameter ξ is discussed in Sec. 1];
no mean vertical velocity field in a Galilean reference frame moving with the mixing layer (⁠ $\bar{w} = \tilde{w} + \bar{w^{″}} = 0$ ⁠, which implies that the mean vertical velocity is given by $\tilde{w} = a_{z}$ ⁠);
incompressible, so that $\partial {\tilde{v}}_{i} / \partial x_{i} = 0$ ⁠;
${\bar{σ}}_{i j} = {\tilde{σ}}_{i j} = 0$ ⁠; and
molecular viscosity and diffusivity contributions to the turbulent transport are negligible compared to the turbulent contributions, i.e., $ν_{t} ≫ \bar{ν}$ and $D_{t} ≫ \bar{D}$ ⁠, where $D_{t} = ν_{t} / {Sc}_{t}$ is the turbulent diffusivity (⁠ ${Sc}_{t}$ is the constant turbulent Schmidt number, taken to be σ_S here).

Therefore, in the small Atwood number limit of a statistically one-dimensional RT, RM, or KH mixing layer, the transport equations can be simplified by neglecting the gradients in the homogeneous (x and y) directions, so that the resulting model only considers variations across the mixing layer (the vertical z-direction in the convention adopted here). The mean pressure gradient in Eqs. (4) and (5) is taken to be in hydrostatic equilibrium

\frac{\partial \bar{p}}{\partial z} = - \bar{ρ} g_{0}

(17)

(⁠ $g_{0} =$ constant > 0 for canonical RT flow). This approximation decouples the mean momentum equation (2) if advection is neglected and was used previously.^5–7 Numerical simulations showed that the mean pressure satisfies hydrostatic equilibrium closely at sufficiently late times (for example, see Ref. 34). The mean pressure gradient term in Eq. (10) is neglected as in the incompressible limit the mean mixture sound speed is $c_{s} = \sqrt{\bar{γ} \bar{p} / \bar{ρ}} ↑ \infty$ ⁠. The resulting simplified equations are (A1)–(A7) for all three instability cases.

C. Self-similarity analysis

Analytic self-similar solutions were obtained for the K–ε model,^5–7 the K–L model,¹⁴ the K–L–a model,^17,18 and the K–2L–a model.¹⁹ The two- and three-equation model studies only considered Rayleigh–Taylor and Richtmyer–Meshkov mixing, and the four-equation model considered mixing induced by all three instabilities considered here. The steps required to obtain these solutions are as follows:

assume mean fields that are linear in the similarity variable $η = η (z, t)$ and inverse parabolic in η for the turbulent fields;
for each turbulent field, assume a separable form consisting of the assumed spatiotemporal profile multiplied by a power-law in time and a prefactor, $ϕ_{α} (z, t) = A_{ϕ_{α}} f (η) t^{p_{ϕ_{α}}}$ ⁠;
substitute these fields into the simplified one-dimensional mean and turbulence model equations;
require each resulting equation to be explicitly time-independent, which determines the exponents ${p_{ϕ_{α}}}$ ⁠;
in these explicitly time-independent equations require the constant terms and terms proportional to powers of η to separately equal zero; and
solve these equations for the turbulent Schmidt numbers ${σ_{ϕ_{α}}}$ ⁠, mixing layer parameters and exponents appearing in the mixing layer widths, and prefactors ${A_{ϕ_{α}}}$ as functions of the model coefficients.

This procedure can be performed in (z, t)-space^8,14,18,19 or in $η$ -space^5–7 and yield equivalent results. As will be shown in Appendix A, strict self-similarity imposes a strong constraint on the turbulent Schmidt numbers: they must be equal, Eq. (A48). For additional discussion of self-similarity assumptions and solution methods, see Refs. 46 and 47.

1. General mixing layer width

For each instability, the self-similar mixing layer width grows as a generalized power-law

h (t) = \frac{ξ}{2} A_{h} t^{θ_{h}},

(18)

where A_h is a constant depending on characteristic parameters for each instability and model coefficients, and the dimensional considerations given the flow parameters in A_h determine the value of the exponent θ_h. The width (18) corresponds to the statistically one-dimensional growth of a three-dimensional mixing layer arising from multimode initial perturbations in a fully developed turbulent state. Note that the corresponding mixing layer velocity and acceleration are

v_{h} (t) = \frac{d h}{d t} = \frac{ξ}{2} θ_{h} A_{h} t^{θ_{h} - 1} = \frac{θ_{h} h (t)}{t} .

(19)

a_{h} (t) = \frac{d^{2} h}{d t^{2}} = \frac{ξ}{2} θ_{h} (θ_{h} - 1) A_{h} t^{θ_{h} - 2} = \frac{θ_{h} (θ_{h} - 1) h (t)}{t^{2}},

(20)

respectively.

The dimensionless similarity variable is

η (z, t) = \frac{ξ z}{h (t)} = \frac{2 z}{A_{h} t^{θ_{h}}} \in [- 1, 1]

(21)

with ξ = 1 and 2 corresponding to the interpretation of h as the layer half-width (i.e., h_b) and width, respectively. From the definition of $η (z, t)$ ⁠, it follows that the temporal and spatial partial derivatives transform as

\frac{\partial}{\partial t} = \frac{\partial η}{\partial t} \frac{d}{d η} = - \frac{θ_{h} η}{t} \frac{d}{d η}, \frac{\partial}{\partial z} = \frac{\partial η}{\partial z} \frac{d}{d η} = \frac{2}{A_{h} t^{θ_{h}}} \frac{d}{d η} .

(22)

2. Self-similar turbulent field profiles

The self-similar turbulent fields have the general form

ϕ_{α} (z, t) = A_{ϕ_{α}} f (\frac{ξ z}{h (t)}) t^{p_{ϕ_{α}}}

(23)

with inverted parabolic profile

\begin{matrix} f (\frac{ξ z}{h (t)}) = {\begin{matrix} 1 - ξ^{2} z^{2} / h {(t)}^{2} & if | z | < h (t) / ξ \\ 0 & if | z | \geq h (t) / ξ \end{matrix} \\ = {\begin{matrix} 1 - η^{2} & if | η | < 1, \\ 0 & if | η | \geq 1, \end{matrix} \end{matrix}

(24)

which is symmetric about the centerline $z = η = 0$ ⁠. The requirement that $f (\cdot)$ is zero outside of the mixing layer is needed for realizability. The temporal and spatial derivatives

\frac{\partial ϕ_{α}}{\partial t} = \frac{p_{ϕ_{α}} ϕ_{α} (z, t)}{t} + 2 A_{ϕ_{α}} θ_{h} \frac{ξ^{2} z^{2}}{h {(t)}^{2}} t^{p_{ϕ_{α}} - 1}

(25)

\frac{\partial ϕ_{α}}{\partial z} = - 2 A_{ϕ_{α}} \frac{ξ^{2} z}{h {(t)}^{2}} t^{p_{ϕ_{α}}}

(26)

follow directly from Eqs. (23) and (24). The gradients of the turbulent fields are linear and are positive for z < 0 and negative for z > 0.

3. Self-similar mean and turbulent fields

The mean and turbulent fields are expressed in the separable form

\begin{matrix} \tilde{s} (z, t) = \tilde{s} (η) f_{s} (t), \bar{v} (z, t) = \bar{v} (η) f_{v} (t), \\ \tilde{U} (z, t) = \tilde{U} (η) f_{U} (t), \\ K (z, t) = K (η) f_{K} (t), ε (z, t) = ε (η) f_{ε} (t), \\ S (z, t) = S (η) f_{S} (t), χ (z, t) = χ (η) f_{χ} (t), \end{matrix}

(27)

where the time-dependent functions ${f_{ϕ_{α}} (t)} = {t^{p_{ϕ_{α}}}}$ are determined by requiring that the transformed self-similar equations be explicitly independent of time (see Appendix A). For the profile Eq. (24), the temporal and spatial derivatives (25) and (26) become

\begin{matrix} \frac{\partial ϕ_{α}}{\partial t} = A_{ϕ_{α}} [f (η) \frac{d f_{ϕ_{α}}}{d t} + 2 θ_{h} \frac{f_{ϕ_{α}} (t)}{t} η^{2}] \\ = A_{ϕ_{α}} [\frac{d f_{ϕ_{α}}}{d t} - θ_{h} \frac{f_{ϕ_{α}} (t)}{t} η \frac{d}{d η}] f (η), \end{matrix}

(28)

\frac{\partial ϕ_{α}}{\partial z} = - 2 ξ A_{ϕ_{α}} f_{ϕ_{α}} (t) \frac{η}{h (t)} = \frac{2 A_{ϕ_{α}} f_{ϕ_{α}} (t)}{A_{h} t^{θ_{h}}} \frac{d f}{d η} .

(29)

The coupled nonlinear equations have compactly supported inverted parabolic solutions

(\begin{matrix} \tilde{U} \\ K \\ ε \\ S \\ χ \end{matrix}) (η) = (\begin{matrix} A_{U} \\ A_{K} \\ A_{ε} \\ A_{S} \\ A_{χ} \end{matrix}) (1 - η^{2})

(30)

with Dirichlet boundary conditions $ϕ_{α} (η = \pm 1) = 0$ ⁠. The constant prefactors ${A_{ϕ_{α}}}$ are functions of the model coefficients and other parameters (e.g., g₀, At etc.). The details of the determination of the relationships between model coefficients and the field prefactors in (30) are given in Appendix A.

The maximum or minimum value of a turbulent field is given by its centerline value

ϕ_{α} (0, t) = A_{ϕ_{α}} t^{p_{ϕ_{α}}} .

(31)

Using Eq. (24) and the profile integrals,

\begin{matrix} \frac{ξ}{2 h (t)} \int_{- h (t) / ξ}^{h (t) / ξ} f (\frac{ξ z}{h (t)}) d z = \frac{2}{3}, \\ \frac{ξ}{2 h (t)} \int_{- h (t) / ξ}^{h (t) / ξ} \sqrt{f (\frac{ξ z}{h (t)})} d z = \frac{π}{4}, \end{matrix}

(32)

it follows that the time-dependent turbulent fields integrated over the mixing layer are

ϕ_{α} (t) = \frac{2}{3} A_{ϕ_{α}} t^{p_{ϕ_{α}}}, L (t) = \frac{π}{4} A_{L} t^{p_{L}} .

(33)

The simplest, lowest order linear mean field profiles common to all three instabilities are the mean scalar for miscible mixing

\tilde{s} (η) = \frac{1 + η}{2}

(34)

satisfying $\tilde{s} (η = - 1) = 0$ and $\tilde{s} (η = 1) = 1$ for $s = m_{1}$ ⁠, c, and ξ. The mean density for miscible mixing is

\bar{ρ} (η) = ρ_{0} + \frac{Δ ρ}{2} η = ρ_{0} (1 + At η)

(35)

with average density $ρ_{0} = (ρ_{1} + ρ_{2}) / 2$ and $Δ ρ = ρ_{1} - ρ_{2} > 0$ satisfying $\bar{ρ} (η = - 1) = ρ_{2}$ and $\bar{ρ} (η = 1) = ρ_{1}$ ⁠. In the self-similar equations, the mean field gradients are

\frac{d \tilde{s}}{d η} = \frac{1}{2}, \frac{d \bar{ρ}}{d η} = \frac{Δ ρ}{2} .

(36)

Using Eq. (18), dimensional analysis gives the general temporal scalings

\begin{matrix} K (t) \propto {(\frac{d h}{d t})}^{2} \propto {(\frac{ξ}{2} θ_{h} A_{h})}^{2} t^{2 θ_{h} - 2}, \\ ε (t) \propto \frac{1}{t} {(\frac{d h}{d t})}^{2} \propto {(\frac{ξ}{2} θ_{h} A_{h})}^{2} t^{2 θ_{h} - 3}, \\ S (t) \propto 1, χ (t) \propto \frac{1}{t} (for a mass fraction variance), \\ a_{z} (t) \propto \frac{d h}{d t} \propto \frac{ξ}{2} θ_{h} A_{h} t^{θ_{h} - 1}, \\ ν_{t} (t) \propto h \frac{d h}{d t} \propto θ_{h} {(\frac{ξ}{2} A_{h})}^{2} t^{2 θ_{h} - 1}, \\ L (t) \propto h, τ_{m} (t), τ_{s} (t) \propto t . \end{matrix}

(37)

For each instability considered here, the mean and turbulent fields can be used to analytically reconstruct the self-similar time-evolution, production, dissipation/destruction, turbulent diffusion, and advection terms in the one-dimensional transport equations

\frac{\partial ϕ_{α} (z, t)}{\partial t} = P_{ϕ_{α}} (z, t) - D_{ϕ_{α}} (z, t) + T_{ϕ_{α}} (z, t) - A_{ϕ_{α}} (z, t) .

(38)

In the gradient-diffusion approximation, the turbulent diffusion term is the negative gradient of the associated turbulent field flux,

\begin{matrix} T_{ϕ_{α}} (z, t) = - \frac{\partial F_{ϕ_{α}} (z, t)}{\partial z}, \\ F_{ϕ_{α}} (z, t) = \bar{ϕ_{α}^{'} w^{'}} = - \frac{ν_{t}}{σ_{ϕ_{α}}} \frac{\partial ϕ_{α}}{\partial z} = - \frac{C_{μ} A_{K}^{2} A_{ϕ_{α}}}{σ_{ϕ_{α}} A_{ε}} f \frac{d f}{d z} t^{2 θ_{K} - θ_{ε} + p_{ϕ_{α}}} \\ = \frac{2 C_{μ} A_{K}^{2} A_{ϕ_{α}}}{σ_{ϕ_{α}} A_{ε}} \frac{ξ^{2} z}{h {(t)}^{2}} [1 - \frac{ξ^{2} z^{2}}{h {(t)}^{2}}] t^{2 θ_{K} - θ_{ε} + p_{ϕ_{α}}}, \end{matrix}

(39)

which has an approximate negative sinusoidal profile with $F_{ϕ_{α}} (z, t) < 0$ for z < 0 and $F_{ϕ_{α}} (z, t) > 0$ for z > 0. The advection term is

\begin{matrix} A_{ϕ_{α}} (z, t) = \tilde{w} \frac{\partial ϕ_{α}}{\partial z} = a_{z} \frac{\partial ϕ_{α}}{\partial z} = - \frac{ν_{t}}{σ_{ρ} \bar{ρ}} \frac{\partial \bar{ρ}}{\partial z} \frac{\partial ϕ_{α}}{\partial z} \\ = \frac{F_{ϕ_{α}} (z, t)}{\bar{ρ} (z, t)} \frac{σ_{ϕ_{α}}}{σ_{ρ}} \frac{ξ Δ ρ}{2 h (t)} \\ \approx \frac{2 At C_{μ} A_{K}^{2} A_{ϕ_{α}}}{σ_{ρ} A_{ε}} \frac{ξ^{3} z}{h {(t)}^{3}} [1 - \frac{ξ^{2} z^{2}}{h {(t)}^{2}}] t^{2 θ_{K} - θ_{ε} + p_{ϕ_{α}}}, \end{matrix}

(40)

which has the same profile shape as $F_{ϕ_{α}} (z, t)$ ⁠.

III. SELF-SIMILAR SOLUTIONS FOR RAYLEIGH–TAYLOR MIXING

The solutions of the self-similar equations are summarized here; see Appendix A for details on their derivation.

A. Equations for Rayleigh–Taylor mixing

As in previous similarity analyses of incompressible Rayleigh–Taylor mixing, the mean velocity $\tilde{w}$ and shear production $- τ_{z z} \partial \tilde{w} / \partial z$ are assumed to be zero, and the mean density in the turbulent diffusion terms and in the denominator of the buoyancy production terms is approximated by $\bar{ρ} \approx ρ_{0}$ ⁠. Thus, Eqs. (1)–(7) reduce to

\frac{\partial \tilde{s}}{\partial t} = \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{s}} \frac{\partial \tilde{s}}{\partial z}),

(41)

\frac{\partial \tilde{U}}{\partial t} = - \frac{g_{0} ν_{t}}{σ_{ρ} ρ_{0}} \frac{\partial \bar{ρ}}{\partial z} + ε + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{U}} \frac{\partial \tilde{U}}{\partial z}),

(42)

\frac{\partial K}{\partial t} = \frac{g_{0} ν_{t}}{σ_{ρ} ρ_{0}} \frac{\partial \bar{ρ}}{\partial z} - ε + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{K}} \frac{\partial K}{\partial z}),

(43)

\frac{\partial ε}{\partial t} = C_{ε 0} \frac{ε}{K} \frac{g_{0} ν_{t}}{σ_{ρ} ρ_{0}} \frac{\partial \bar{ρ}}{\partial z} - C_{ε 2} \frac{ε^{2}}{K} + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{ε}} \frac{\partial ε}{\partial z}),

(44)

\frac{\partial S}{\partial t} = \frac{2 ν_{t}}{σ_{s}} {(\frac{\partial \tilde{s}}{\partial z})}^{2} - 2 χ + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{S}} \frac{\partial S}{\partial z}),

(45)

\begin{matrix} \frac{\partial χ}{\partial t} = C_{χ 0} \frac{ν_{t}}{\sqrt{Sc}} \frac{ε}{K} {(\frac{\partial \tilde{s}}{\partial z})}^{2} - (C_{χ 2} \frac{χ}{S} - C_{χ 3} \frac{ε}{K}) χ + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{χ}} \frac{\partial χ}{\partial z}) . \end{matrix}

(46)

Equations (45) and (46) have not been previously considered in self-similarity analysis, but Eq. (45) has been previously considered in the similarity analysis of the K–L–a–V model with V the mass fraction variance.⁴⁸

B. Self-similar turbulent fields

Explicit expressions for the RT self-similar turbulent fields corresponding to the two-, three-, and four-equation models are presented here as a function of the model coefficients.

1. Two- and three-equation models

For the two- and three-equation models, using the solutions derived in Appendix A, Eqs. (A56)–(A59), the turbulent kinetic energy, turbulent kinetic energy dissipation rate, scalar variance, and scalar variance dissipation rate are functions of $f (\cdot)$ given in Eq. (24),

\begin{matrix} K (z, t) = \frac{C_{μ} {(g_{0} At)}^{2} σ_{s}}{σ_{ρ}^{2}} \frac{{(C_{ε 2} - C_{ε 0})}^{3}}{(4 C_{ε 0} - 3) {(4 C_{ε 2} - 3)}^{2}} f (\frac{ξ z}{h (t)}) t^{2}, \end{matrix}

(47)

ε (z, t) = \frac{C_{μ} {(g_{0} At)}^{2} σ_{s}}{σ_{ρ}^{2}} {(\frac{C_{ε 2} - C_{ε 0}}{4 C_{ε 2} - 3})}^{2} f (\frac{ξ z}{h (t)}) t,

(48)

S (z, t) = \frac{C_{ε 2} - C_{ε 0}}{4 [C_{ε 2} - C_{ε 0} + C_{χ} (4 C_{ε 0} - 3)]} f (\frac{ξ z}{h (t)}),

(49)

χ (z, t) = \frac{C_{χ} (4 C_{ε 0} - 3)}{4 [C_{ε 2} - C_{ε 0} + C_{χ} (4 C_{ε 0} - 3)]} f (\frac{ξ z}{h (t)}) \frac{1}{t} .

(50)

Note that if $4 C_{ε 0} - 3 > 0$ ⁠, then $C_{ε 2} > C_{ε 0}$ in order for $K (z, t) > 0$ ⁠. Morgan et al.⁴⁸ obtained $A_{S} = (1 - θ_{m}) / 4$ (θ_m is the molecular mixing parameter) but not its dependence on the model coefficients.

The mass flux velocity, turbulent viscosity, turbulent lengthscale, mechanical turbulent timescale, and scalar turbulent timescale are

a_{z} (z, t) = - \frac{C_{μ} g_{0} {At}^{2} {(C_{ε 2} - C_{ε 0})}^{2} σ_{s}}{(4 C_{ε 0} - 3) (4 C_{ε 2} - 3) σ_{ρ}^{2}} f (\frac{ξ z}{h (t)}) t,

(51)

\begin{matrix} ν_{t} (z, t) = \frac{{(C_{μ} g_{0} At)}^{2} σ_{s}}{σ_{ρ}^{2}} \frac{{(C_{ε 2} - C_{ε 0})}^{4}}{{(4 C_{ε 0} - 3)}^{2} {(4 C_{ε 2} - 3)}^{2}} f (\frac{ξ z}{h (t)}) t^{3}, \end{matrix}

(52)

L (z, t) = \frac{C_{ε} \sqrt{C_{μ} σ_{s}} {(C_{ε 2} - C_{ε 0})}^{5 / 2} g_{0} At}{σ_{ρ} {(4 C_{ε 0} - 3)}^{3 / 2} (4 C_{ε 2} - 3)} \sqrt{f (\frac{ξ z}{h (t)})} t^{2},

(53)

τ_{m} (t) = \frac{C_{ε 2} - C_{ε 0}}{4 C_{ε 0} - 3} t, τ_{s} (t) = \frac{1}{2 C_{χ}} \frac{C_{ε 2} - C_{ε 0}}{4 C_{ε 0} - 3} t .

(54)

Note that $a_{z} (z, t) < 0$ and the turbulent timescales are independent of z. The scalar fields depend on the coefficients ${C_{ε q}}$ in the ε equation, while the mechanical fields are independent of S and χ (i.e., the scalar fields are passive).

2. Four-equation model

For the four-equation model, not assuming an algebraic model (9), the coefficients A_S and $A_{χ}$ in

S (z, t) = A_{S} f (\frac{ξ z}{h (t)}),

(55)

χ (z, t) = A_{χ} f (\frac{ξ z}{h (t)}) \frac{1}{t}

(56)

are given by Eqs. (A73) and (A74), respectively. The scalar turbulent timescale is

τ_{s} (t) = \frac{A_{S}}{2 A_{χ}} t .

(57)

These expressions are the same for RM and KH mixing, but with the mixing layer width h(t) and coefficients A_S and $A_{χ}$ corresponding to those instabilities.

C. Self-similar parameters

Explicit expressions for the RT self-similar mixing layer growth parameter, ratio of generated turbulent kinetic energy to released potential energy, and molecular mixing parameter corresponding to the three- and four-equation models are presented here as a function of the model coefficients.

1. Mixing layer growth parameter

For small Atwood number Rayleigh–Taylor instability, the self-similar mixing layer width grows as a power-law (18) with $θ_{h} = 2$ for a constant acceleration,

h (t) = \frac{ξ}{2} A_{h} t^{2} = \frac{ξ}{2} α At g_{0} t^{2}

(58)

with $A_{h} = α$ At g₀. The mixing layer velocity and acceleration are

v_{h} (t) = ξ α At g_{0} t, a_{h} (t) = ξ α At g_{0} .

(59)

Using Eq. (A40), the dimensionless mixing layer growth parameter is^5–8

\begin{matrix} α (C_{μ}, C_{ε 0}, C_{ε 2}, σ_{ρ}) = \frac{1}{g_{0} At} \sqrt{\frac{4 C_{μ} A_{K}^{2}}{σ_{s} A_{ε}}} \\ = \frac{2 C_{μ}}{σ_{ρ}} \frac{{(C_{ε 2} - C_{ε 0})}^{2}}{(4 C_{ε 0} - 3) (4 C_{ε 2} - 3)}, \end{matrix}

(60)

which exhibits the explicit dependence on the model coefficients; $α_{b} = α / 2$ ⁠. Note that α is degenerate, in that a given value can be obtained using different combinations of the coefficients. The coefficients must satisfy the inequalities $C_{ε 2}, C_{ε 0} > 3 / 4$ to ensure that $α > 0$ ⁠; the choice $C_{ε 0} = C_{ε 2}$ is excluded as it results in no growth.

It can be shown that the buoyancy production term in the ε equation cannot be set to zero as follows. If $C_{ε 0} = 0$ ⁠, Eq. (60) simplifies to $α = [2 C_{μ} / (3 σ_{ρ})] C_{ε 2}^{2} / (3 - 4 C_{ε 2})$ and positivity requires that $C_{ε 2} < 3 / 4$ ⁠, which violates the requirement that $C_{ε 2} > 1$ for positivity of the RM turbulent timescales [see Eq. (82)] and for decaying isotropic turbulence (see Appendix B). Although similarity requires $σ_{ρ} = σ_{s}$ ⁠, this will not be assumed in the presentation of the results for generality and to allow for additional model calibration options.

As expected from dimensional analysis, $K \propto {(d h / d t)}^{2} \propto {(g_{0} At t)}^{2}, ε \propto {(d h / d t)}^{2} / t \propto {(g_{0} At)}^{2} t$ ⁠, S is dimensionless (for a mass fraction variance), $χ \propto 1 / t, a_{z} \propto d h / d t \propto g_{0}$ At $^{2} t, ν_{t} \propto hdh / d t \propto {(g_{0} At)}^{2} t^{3}, L \propto h \propto g_{0}$ At $t^{2}$ ⁠, and $τ_{m}, τ_{s} \propto t$ ⁠. In particular, ${K, ε, a_{z}, ν_{t}, L} ↓ 0$ as g₀At $↓ 0$ ⁠, as required in the limit of a stable interface. Using Eqs. (53), (58), and (33), the turbulent lengthscale-to-mixing layer width ratio is

\frac{L (t)}{h (t)} = \frac{C_{ε} π}{4 ξ} \sqrt{\frac{(C_{ε 2} - C_{ε 0}) σ_{s}}{(4 C_{ε 0} - 3) C_{μ}}} .

(61)

2. Ratio of generated turbulent kinetic energy to released potential energy

The dimensionless ratio of the generated turbulent kinetic energy to released potential energy, $ρ_{0} K / Δ Φ$ ⁠, has been obtained from RT numerical simulations^49,50 and experiments,^51,52 and has been used to calibrate turbulent lengthscale-based Reynolds-averaged models for RT mixing.^14,18,19

The released potential energy is

\begin{matrix} Δ Φ (t) = Φ (0) - Φ (t) \\ = \frac{ξ g_{0}}{2 h (t)} [\int_{- h / ξ}^{0} ρ_{2} z d z + \int_{0}^{h / ξ} ρ_{1} z d z - \int_{- h / ξ}^{h / ξ} \bar{ρ} z d z] \\ = \frac{ρ_{1} - ρ_{2}}{12 ξ} g_{0} h (t) = \frac{ρ_{0}}{12} α {(g_{0} At t)}^{2}, \end{matrix}

(62)

which is proportional to α. Integrating Eq. (47) over the mixing layer using Eq. (32) gives

K (t) = \frac{2}{3} \frac{C_{μ} {(g_{0} At)}^{2} σ_{s}}{σ_{ρ}^{2}} \frac{{(C_{ε 2} - C_{ε 0})}^{3}}{(4 C_{ε 0} - 3) {(4 C_{ε 2} - 3)}^{2}} t^{2},

(63)

and it follows that the ratio of generated turbulent kinetic energy to released potential energy is

\frac{ρ_{0} K}{Δ Φ} (C_{ε 0}, C_{ε 2}, σ_{s}, σ_{ρ}) = \frac{4 σ_{s}}{σ_{ρ}} \frac{C_{ε 2} - C_{ε 0}}{4 C_{ε 2} - 3} .

(64)

Note that $ρ_{0} K / Δ Φ > 0$ is singular for $C_{ε 2} = 3 / 4$ ⁠, and the inequality $(C_{ε 2} - C_{ε 0}) / (4 C_{ε 2} - 3) > 0$ should be satisfied; as $C_{ε 2} > 3 / 4$ ⁠, it follows that $C_{ε 2} > C_{ε 0}$ ⁠. Both Eqs. (60) and (64) depend on the difference $C_{ε 2} - C_{ε 0}$ ⁠, the magnitude of which is related to the competition between turbulent kinetic energy dissipation rate buoyancy production (i.e., $C_{ε 0}$ ⁠) and destruction (i.e., $C_{ε 2}$ ⁠).

3. Molecular mixing parameter

The molecular mixing parameter

θ_{m} = 1 - \frac{S (z, t)}{\tilde{s} (z, t) [1 - \tilde{s} (z, t)]} = 1 - 4 A_{S}

(65)

is a constant (and positivity requires $A_{S} < 1 / 4$ ⁠), as $\tilde{s} (1 - \tilde{s}) = f (ξ z / h (t)) / 4$ ⁠, and is also equal to the spatially integrated expression (in the self-similar state),

θ_{m} = 1 - \frac{\int_{- h / ξ}^{h / ξ} S (z, t) d z}{\int_{- h / ξ}^{h / ξ} \tilde{s} (z, t) [1 - \tilde{s} (z, t)] d z} .

These expressions also apply to Richtmyer–Meshkov and Kelvin–Helmholtz mixing.

Using Eq. (A58), it follows that the molecular mixing parameter is

θ_{m} (C_{ε 0}, C_{ε 2}, C_{χ}) = \frac{(4 C_{ε 0} - 3) C_{χ}}{C_{ε 2} - C_{ε 0} + (4 C_{ε 0} - 3) C_{χ}}

(66)

for the three-equation model, which must be positive. For specific, fixed values of $C_{ε 0}$ and $C_{ε 2}$ ⁠, θ_m is determined by the value of $C_{χ}$ ⁠.

For the four-equation model,

θ_{m} (C_{ε 0}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = 1 - 4 A_{S}

(67)

with A_S given by Eq. (A73). Expression (66) is considerably simpler than (67) and depends only on three coefficients. This expression is the same for RM and KH mixing, but with the coefficient A_S corresponding to these instabilities.

D. Production-to-dissipation/destruction and mechanical-to-scalar timescale ratios

Explicit expressions for the RT self-similar production-to-dissipation/destruction and mechanical-to-scalar timescale ratios corresponding to the two-, three-, and four-equation models are presented here as a function of the model coefficients.

1. Production-to-dissipation/destruction ratios

The turbulent kinetic energy production-to-dissipation ratio is

\frac{P_{K}}{D_{K}} (C_{ε 0}, C_{ε 2}) = - \frac{g_{0} a_{z}}{ε} = \frac{4 C_{ε 2} - 3}{4 C_{ε 0} - 3},

(68)

and the turbulent kinetic energy dissipation rate production-to-destruction ratio is

\frac{P_{ε}}{D_{ε}} (C_{ε 0}, C_{ε 2}) = \frac{C_{ε 0}}{C_{ε 2}} \frac{P_{K}}{D_{K}} = \frac{C_{ε 0}}{C_{ε 2}} \frac{4 C_{ε 2} - 3}{4 C_{ε 0} - 3} .

(69)

As $C_{ε 2} > C_{ε 0}, P_{ε} / D_{ε} < P_{K} / D_{K} > 1$ ⁠. For the three-equation model, the scalar variance production-to-dissipation ratio is

\frac{P_{S}}{D_{S}} (C_{ε 0}, C_{ε 2}, C_{χ}) = \frac{ν_{t} {(\partial \tilde{s} / \partial z)}^{2}}{σ_{s} χ} = 1 + \frac{C_{ε 2} - C_{ε 0}}{(4 C_{ε 0} - 3) C_{χ}} .

(70)

These quantities must be positive and are singular for $C_{ε 0} = 3 / 4$ ⁠. Increasing the value of $C_{χ}$ decreases the value of $P_{S} / D_{S} > 1$ ⁠.

For the four-equation model,

\frac{P_{S}}{D_{S}} (C_{ε 0}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{1}{4 A_{χ}},

(71)

and the scalar variance dissipation rate production-to-destruction ratio is

\frac{P_{χ}}{D_{χ}} (C_{ε 0}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{1}{C_{χ 2}} [C_{χ 0} \frac{ν_{t} {(\partial \tilde{s} / \partial z)}^{2}}{\sqrt{Sc} χ} + C_{χ 3}] \frac{ε S}{K χ} = \frac{1}{C_{χ 2}} (\frac{C_{χ 0} σ_{s}}{4 \sqrt{Sc} A_{χ}} + C_{χ 3}) \frac{4 C_{ε 0} - 3}{C_{ε 2} - C_{ε 0}} \frac{A_{S}}{A_{χ}}

(72)

with A_S and $A_{χ}$ given by Eqs. (A73) and (A74). These quantities must be positive and are singular if the denominators vanish.

2. Mechanical-to-scalar timescale ratio

For the three-equation model, the mechanical-to-scalar timescale ratio is (15) and for the four-equation model,

R (C_{ε 0}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = 2 \frac{A_{K}}{A_{ε}} \frac{A_{χ}}{A_{S}} = \frac{2 (C_{ε 2} - C_{ε 0})}{4 C_{ε 0} - 3} \frac{A_{χ}}{A_{S}},

(73)

which is singular for $C_{ε 0} = 3 / 4$ ⁠, with A_S and $A_{χ}$ given by Eqs. (A73) and (A74).

IV. SELF-SIMILAR SOLUTIONS FOR RICHTMYER–MESHKOV MIXING

The solutions of the self-similar equations are summarized here; see Appendix A for details on their derivation.

A. Equations for Richtmyer–Meshkov mixing

The self-similar equations describing the growth of incompressible Richtmyer–Meshkov mixing are formally the same as those for Rayleigh–Taylor mixing (41)–(46), but with an impulsive, rather than constant, acceleration:

g_{0} (t) = Δ v_{s} δ (t),

(74)

where $Δ v_{s} > 0$ is the velocity jump due to the impulse (or “shock” passage) across the initial interface, and $δ (\cdot)$ is the Dirac delta function. Thus, the “shock” production terms exist only at t = 0 (the time of “shock” passage through the interface), and the equations subsequently describe a decaying turbulence state. Note that this decaying turbulence is inhomogeneous due to the turbulent diffusion terms (which are absent in homogeneous, isotropic decaying turbulence as discussed in Appendix B 1). The production terms set the initial conditions, i.e., $Δ v_{s}$ ⁠.

B. Self-similar turbulent fields

Explicit expressions for the RM self-similar turbulent fields corresponding to the two-, three-, and four-equation models are presented here as a function of the model coefficients.

1. Two- and three-equation models

For the two- and three-equation models, using the solutions derived in Appendix A, Eqs. (A60)–(A63), the turbulent kinetic energy, turbulent kinetic energy dissipation rate, scalar variance, and scalar variance dissipation rate are functions of $f (\cdot)$ given in Eq. (24),

\begin{matrix} K (z, t) = \frac{h_{0}^{2} θ σ_{s} (2 - 3 θ)}{8 C_{μ}} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{2 θ} f (\frac{ξ z}{h (t)}) t^{2 θ - 2}, \end{matrix}

(75)

ε (z, t) = \frac{h_{0}^{2} θ σ_{s} {(2 - 3 θ)}^{2}}{8 C_{μ}} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{2 θ} f (\frac{ξ z}{h (t)}) t^{2 θ - 3},

(76)

S (z, t) = \frac{θ}{4 [θ + 2 C_{χ} (2 - 3 θ)]} f (\frac{ξ z}{h (t)}),

(77)

χ (z, t) = \frac{C_{χ} (2 - 3 θ) θ}{4 [θ + 2 C_{χ} (2 - 3 θ)]} f (\frac{ξ z}{h (t)}) \frac{1}{t} .

(78)

The mass flux velocity, turbulent viscosity, turbulent lengthscale, mechanical turbulent timescale, and scalar turbulent timescale are

a_{z} (z, t) = - \frac{{At}^{+} h_{0} θ σ_{s}}{4 σ_{ρ}} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{θ} f (\frac{ξ z}{h (t)}) t^{θ - 1},

(79)

ν_{t} (z, t) = \frac{h_{0}^{2} θ σ_{s}}{8} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{2 θ} f (\frac{ξ z}{h (t)}) t^{2 θ - 1},

(80)

L (z, t) = C_{ε} \sqrt{\frac{θ σ_{s}}{8 (2 - 3 θ) C_{μ}}} h_{0} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{θ} \sqrt{f (\frac{ξ z}{h (t)})} t^{θ},

(81)

τ_{m} (t) = (C_{ε 2} - 1) t, τ_{s} (t) = \frac{C_{ε 2} - 1}{2 C_{χ}} t .

(82)

Positivity of the turbulent timescales requires $C_{ε 2} > 1$ ⁠.

2. Four-equation model

For the four-equation model, not assuming an algebraic model (9), the coefficients A_S and $A_{χ}$ in Eqs. (55) and (56) are given by Eqs. (A75) and (A76), respectively. The scalar turbulent timescale is given by Eq. (57) with these coefficients.

C. Self-similar parameters

Explicit expressions for the RM self-similar mixing layer growth exponent, molecular mixing parameter corresponding to the three- and four-equation models, and mechanical and scalar field decay exponents are presented here as a function of the model coefficients.

1. Mixing layer growth exponent

For small Atwood number Richtmyer–Meshkov instability, the self-similar mixing layer width grows as a power-law (18) with $θ_{h} = θ$ ⁠,

h (t) = \frac{ξ}{2} A_{h} t^{θ} = \frac{ξ}{2} h_{0} {(\frac{Δ v_{s} {At}^{+} t}{θ h_{0}})}^{θ}

(83)

with $A_{h} = h_{0} / t_{0}^{θ} = h_{0}^{1 - θ} (Δ v_{s}$ At $^{+} / θ)^{θ}, t_{0} = θ h_{0} / (Δ v_{s}$ At $^{+})$ ⁠, At⁺ is the post-shock Atwood number, and $h (t_{0}) = (ξ / 2) h_{0}$ ⁠. The mixing layer velocity and acceleration are

v_{h} (t) = \frac{ξ}{2} {(\frac{Δ v_{s} {At}^{+} t}{θ h_{0}})}^{θ - 1} Δ v_{s} {At}^{+}, a_{h} (t) = \frac{ξ}{2} (θ - 1) {(\frac{Δ v_{s} {At}^{+} t}{θ h_{0}})}^{θ - 2} \frac{{(Δ v_{s} {At}^{+})}^{2}}{θ h_{0}} .

(84)

The dimensionless mixing layer growth exponent that follows from inverting Eq. (A65) is^5,7,8

θ (C_{ε 2}) = \frac{2 C_{ε 2} - 3}{3 (C_{ε 2} - 1)},

(85)

and only depends on $C_{ε 2}$ ⁠; growth requires $C_{ε 2} > 3 / 2$ ⁠. The growth exponent is singular for $C_{ε 2} = 1$ ⁠.

Note that h(t) depends on a timescale t₀ and initial lengthscale $h_{0}$ ⁠, as distinct from RT mixing; h₀ is interpreted as an initial width of a mixing layer front, and not an amplitude (such as in single-mode linear and nonlinear theory). The dependence on At⁺ ensures that there is no mixing layer growth if the “post-shock” densities are equal (i.e., At $^{+} = 0$ ⁠). The dependence of t₀ on θ is motivated by requiring that the initial velocity due to the impulse is $d h / d t |_{t_{0}} = (ξ / 2) Δ v_{s}$ At⁺.

A relationship between the growth rate of RM mixing and the decay of turbulent kinetic energy was proposed⁵³ using the power-laws predicted by the K–ε model for decaying isotropic turbulence (see Appendix B 4) (where, however, θ was incorrectly interpreted as θ_b and was derived for homogeneous decaying turbulence—different from the expression for RM mixing). The analytic result $θ = 1 / 3$ for three-dimensional RM instability was also derived,⁵⁴ and is consistent with recent numerical simulation⁵⁵ and analytical⁵⁶ results. Note that $a_{h} (t) < 0$ with $θ < 1$ ⁠.

As expected from dimensional analysis, $K \propto {(d h / d t)}^{2} \propto {(h_{0} / t)}^{2} {(t / t_{0})}^{2 θ}, ε \propto {(d h / d t)}^{2} / t \propto (h_{0}^{2} / t^{3}) {(t / t_{0})}^{2 θ}$ ⁠, S is dimensionless (for a mass fraction variance), $χ \propto 1 / t, a_{z} \propto d h / d t \propto (h_{0} / t) {(t / t_{0})}^{θ}, ν_{t} \propto hdh / d t \propto (h_{0}^{2} / t) {(t / t_{0})}^{2 θ}, L \propto h \propto h_{0} {(t / t_{0})}^{θ}$ ⁠, and $τ_{m}, τ_{s} \propto t$ ⁠. In particular, ${K, ε, a_{z}, ν_{t}, L} ↓ 0$ as $Δ v_{s}$ At $^{+} ↓ 0$ ⁠, as required in the limit of a stable interface. Using Eqs. (81), (83), and (33), the turbulent lengthscale-to-mixing layer width ratio is

\frac{L (t)}{h (t)} = \frac{C_{ε} π}{4 ξ} \sqrt{\frac{θ σ_{s}}{2 (2 - 3 θ) C_{μ}}} .

(86)

2. Molecular mixing parameter

Using Eqs. (65) and (A62), the molecular mixing parameter is

θ_{m} (C_{ε 2}, C_{χ}) = \frac{2 (2 - 3 θ) C_{χ}}{θ + 2 (2 - 3 θ) C_{χ}} = \frac{6 C_{χ}}{2 C_{ε 2} - 3 + 6 C_{χ}}

(87)

for the three-equation model [ $C_{ε 2} = C_{ε 2} (θ)$ ] and must be positive.

For the four-equation model,

θ_{m} (C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = 1 - 4 A_{S}

(88)

with A_S given by Eq. (A75). Expression (87) is considerably simpler than (88) and depends only on two coefficients. For a specific, fixed value of $C_{ε 2}$ (or θ), θ_m is determined by the value of $C_{χ}$ ⁠.

3. Decay exponents of mechanical and scalar turbulent fields

Following “shock” passage at t = 0, the production terms vanish and the turbulence is in a decaying state described by dissipation/destruction and diffusion. The power-law decay exponents of the mechanical–scalar turbulent fields follow by integrating Eqs. (75)–(78) over the mixing layer using Eq. (32) so that

K_{0} (t) = \frac{h_{0}^{3} θ σ_{s} (2 - 3 θ)}{12 C_{μ}} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{3 θ} t^{3 θ - 2},

(89)

ε_{0} (t) = \frac{h_{0}^{3} θ σ_{s} {(2 - 3 θ)}^{2}}{12 C_{μ}} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{3 θ} t^{3 θ - 3},

(90)

S_{0} (t) = \frac{h_{0} θ}{6 [θ + 2 C_{χ} (2 - 3 θ)]} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{θ} t^{θ},

(91)

χ_{0} (t) = \frac{h_{0} C_{χ} (2 - 3 θ) θ}{6 [θ + 2 C_{χ} (2 - 3 θ)]} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{θ} t^{θ - 1} .

(92)

For the four-equation model,

S_{0} (t) = \frac{2 A_{S} h_{0}}{3} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{θ} t^{θ},

(93)

χ_{0} (t) = \frac{2 A_{χ} h_{0}}{3} {(\frac{Δ v_{s} {At}^{+}}{θ h_{0}})}^{θ} t^{θ - 1}

(94)

with A_S and $A_{χ}$ given by Eqs. (A75) and (A76).

D. Production-to-dissipation/destruction and mechanical-to-scalar timescale ratios

Explicit expressions for the RM self-similar production-to-dissipation/destruction and mechanical-to-scalar timescale ratios corresponding to the two-, three-, and four-equation models are presented here as a function of the model coefficients.

1. Production-to-dissipation/destruction ratios

After “shock” passage at t = 0, the turbulent kinetic energy production-to-dissipation ratio and turbulent kinetic energy dissipation rate production-to-destruction ratio are zero. For the three-equation model, the scalar variance production-to-dissipation ratio is

\begin{matrix} \frac{P_{S}}{D_{S}} (C_{ε 2}, C_{χ}) = \frac{ν_{t} {(\partial \tilde{s} / \partial z)}^{2}}{σ_{s} χ} = 1 + \frac{θ}{2 (2 - 3 θ) C_{χ}} \\ = 1 + \frac{2 C_{ε 2} - 3}{6 C_{χ}}, \end{matrix}

(95)

and must be positive; this quantity is singular for $θ = 2 / 3$ ⁠. Increasing the value of $C_{χ}$ decreases the value of $P_{S} / D_{S}$ ⁠, which is greater than unity.

For the four-equation model, the scalar variance production-to-dissipation ratio is

\frac{P_{S}}{D_{S}} (C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{θ}{8 A_{χ}},

(96)

and the scalar variance dissipation rate production-to-destruction ratio is

\frac{P_{χ}}{D_{χ}} (C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{1}{C_{χ 2}} [C_{χ 0} \frac{ν_{t} {(\partial \tilde{s} / \partial z)}^{2}}{\sqrt{Sc} χ} + C_{χ 3}] \frac{ε S}{K χ} = \frac{1}{C_{χ 2}} (\frac{C_{χ 0} θ σ_{s}}{8 \sqrt{Sc} A_{χ}} + C_{χ 3}) (2 - 3 θ) \frac{A_{S}}{A_{χ}}

(97)

with A_S and $A_{χ}$ given by Eqs. (A75) and (A76). These quantities must be positive and are singular if the denominators vanish.

2. Mechanical-to-scalar timescale ratio

For the three-equation model, the mechanical-to-scalar timescale ratio is (15) and for the four-equation model,

R (C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{2}{2 - 3 θ} \frac{A_{χ}}{A_{S}},

(98)

which is singular for $θ = 2 / 3$ (and $θ < 2 / 3$ for R > 0), with A_S and $A_{χ}$ given by Eqs. (A75) and (A76).

V. SELF-SIMILAR SOLUTIONS FOR KELVIN–HELMHOLTZ MIXING

The solutions of the self-similar equations are summarized here; see Appendix A for details on their derivation.

A. Equations for Kelvin–Helmholtz mixing

For incompressible Kelvin–Helmholtz instability, the mean momentum equation reduces to an equation for the mean shear velocity $\bar{v}$ ⁠, and $\bar{w} = 0$ ⁠. Equations (1)–(5) reduce to (approximating $\bar{ρ} \approx ρ_{0}$ in the diffusion terms)

\frac{\partial \bar{s}}{\partial t} = \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{s}} \frac{\partial \bar{s}}{\partial z}),

(99)

\frac{\partial \bar{v}}{\partial t} = - \frac{1}{ρ_{0}} \frac{\partial τ_{y z}}{\partial z} = C_{dev} \frac{\partial}{\partial z} (ν_{t} \frac{\partial \bar{v}}{\partial z}),

(100)

\frac{\partial \bar{U}}{\partial t} = ε + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{U}} \frac{\partial \bar{U}}{\partial z}),

(101)

\frac{\partial K}{\partial t} = C_{dev} ν_{t} {(\frac{\partial \bar{v}}{\partial z})}^{2} - ε + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{K}} \frac{\partial K}{\partial z}),

(102)

\frac{\partial ε}{\partial t} = C_{ε 1} C_{dev} \frac{ε}{K} ν_{t} {(\frac{\partial \bar{v}}{\partial z})}^{2} - C_{ε 2} \frac{ε^{2}}{K} + \frac{\partial}{\partial z} (\frac{ν_{t}}{σ_{ε}} \frac{\partial ε}{\partial z}),

(103)

together with Eqs. (45) and (46) for S and χ. The shear Reynolds stress is modeled as

τ_{y z} (z, t) = - C_{dev} μ_{t} \frac{\partial \bar{v}}{\partial z}

(104)

and the mean pressure gradient vanishes by symmetry. The dimensionless coefficient C_dev allows for a different value of $C_{μ}$ in the closure of the shear production terms for KH mixing and was previously introduced in the similarity analysis of the K–2L–a model.¹⁹

The mean shear velocity is assumed to have a linear profile,

\bar{v} (η) = v_{0} + \sqrt{1 - {At}^{2}} \frac{Δ v}{2} η = v_{0} (1 + \sqrt{1 - {At}^{2}} {At}_{v} η)

(105)

with average stream velocity $v_{0} = (v_{1} + v_{2}) / 2$ ⁠, satisfying $\bar{v} (η = - 1) = v_{0} - \sqrt{1 - {At}^{2}} Δ v / 2$ and $\bar{v} (η = 1) = v_{0} + \sqrt{1 - {At}^{2}} Δ v / 2$ [and reducing to $\bar{v} (η = - 1) = v_{2}$ and $\bar{v} (η = 1) = v_{1}$ for At $= 0$ ], where the velocity Atwood number is

{At}_{v} = \frac{v_{1} - v_{2}}{v_{1} + v_{2}} = \frac{Δ v}{2 v_{0}},

(106)

so that (105) has a similar form as the mean density (35) in RT flow. The factor $\sqrt{1 - {At}^{2}}$ (unity in the single-fluid, constant density case) reduces the growth of the instability with increasing density contrast $Δ ρ$ and multiplies $Δ v$ in the expression for the linear instability growth rate.^57–60 The constant mean shear velocity gradient is

\frac{d \bar{v}}{d η} = \sqrt{1 - {At}^{2}} \frac{Δ v}{2},

(107)

and has a similar form as the mean density gradient (36). Note that $\bar{v} (z, t)$ and $\partial \bar{v} (z, t) / \partial z$ are independent of At because $η \propto 1 / h (t) \propto 1 / \sqrt{1 - {At}^{2}}$ [see Eq. (117)]. Note that there is no overall time dependence for the mean shear velocity as $f_{v} = t^{θ_{h} - 1} = 1$ ⁠.

B. Self-similar turbulent fields

Explicit expressions for the KH self-similar turbulent fields corresponding to the two-, three-, and four-equation models are presented here as a function of the model coefficients.

1. Two- and three-equation models

For the two- and three-equation model, using the solutions derived in Appendix A, Eqs. (A66)–(A69), the turbulent kinetic energy, turbulent kinetic energy dissipation rate, scalar variance, and scalar variance dissipation rate are functions of $f (\cdot)$ given in Eq. (24),

K (z, t) = \frac{(C_{ε 2} - C_{ε 1}) {(Δ v)}^{2} (1 - {At}^{2})}{8 C_{ε 2}} f (\frac{ξ z}{h (t)}),

(108)

ε (z, t) = \frac{C_{ε 1} {(Δ v)}^{2} (1 - {At}^{2})}{8 C_{ε 2}} f (\frac{ξ z}{h (t)}) \frac{1}{t},

(109)

S (z, t) = \frac{C_{ε 2} - C_{ε 1}}{4 [C_{ε 2} + (2 C_{χ} - 1) C_{ε 1}]} f (\frac{ξ z}{h (t)}),

(110)

χ (z, t) = \frac{C_{χ} C_{ε 1}}{4 [C_{ε 2} + (2 C_{χ} - 1) C_{ε 1}]} f (\frac{ξ z}{h (t)}) \frac{1}{t} .

(111)

The turbulent viscosity, turbulent lengthscale, and mechanical and scalar turbulent timescale are

ν_{t} (z, t) = \frac{C_{μ} {(C_{ε 2} - C_{ε 1})}^{2} {(Δ v)}^{2} (1 - {At}^{2})}{8 C_{ε 1} C_{ε 2}} f (\frac{ξ z}{h (t)}) t,

(112)

L (z, t) = \frac{C_{ε} {(C_{ε 2} - C_{ε 1})}^{3 / 2} \sqrt{1 - {At}^{2}}}{C_{ε 1} \sqrt{8 C_{ε 2}}} \sqrt{f (\frac{ξ z}{h (t)})} | Δ v | t,

(113)

τ_{m} (t) = \frac{C_{ε 2} - C_{ε 1}}{C_{ε 1}} t, τ_{s} (t) = \frac{1}{2 C_{χ}} \frac{C_{ε 2} - C_{ε 1}}{C_{ε 1}} t .

(114)

The turbulent timescales are independent of z.

The shear Reynolds stress (104) is

τ_{y z} (z, t) = - \frac{ρ_{0} δ {(Δ v)}^{2} (1 - {At}^{2})}{8} f (\frac{ξ z}{h (t)}),

(115)

and the centerline normalized shear stress is

\frac{τ_{y z} (0, t)}{ρ_{0} {(Δ v)}^{2}} = - \frac{δ (1 - {At}^{2})}{8},

(116)

which depends only on δ [see Eq. (119)] and At. The Atwood number dependent factors in Eqs. (108), (109), (112), (113), and (115) reduce the values of each of these quantities with increasing At.

2. Four-equation model

For the four-equation model, not assuming an algebraic model (9), the coefficients A_S and $A_{χ}$ in Eqs. (55) and (56) are given by Eqs. (A77) and (A78), respectively. The scalar turbulent timescale is given by Eq. (57) with these coefficients.

C. Self-similar parameters

Explicit expressions for the KH self-similar mixing layer growth parameter, normalized turbulent kinetic energy and turbulent kinetic energy dissipation rate, and molecular mixing parameter corresponding to the three- and four-equation models are presented here as a function of the model coefficients.

1. Mixing layer growth parameter

For incompressible Kelvin–Helmholtz instability, the self-similar mixing layer width grows as a power-law (18) with $θ_{h} = 1$ ⁠,

h (t) = \frac{ξ}{2} A_{h} t = \frac{ξ}{2} δ \sqrt{1 - {At}^{2}} | Δ v | t

(117)

with $A_{h} = δ \sqrt{1 - {At}^{2}} | Δ v |$ and $Δ v = v_{1} - v_{2}$ the difference between the velocities of the coflowing streams. The mixing layer velocity and acceleration are

v_{h} (t) = \frac{ξ}{2} δ \sqrt{1 - {At}^{2}} | Δ v |, a_{h} (t) = 0 .

(118)

Using Eq. (A40), the dimensionless mixing layer growth parameter is^7,8

\begin{matrix} δ (C_{μ}, C_{dev}, C_{ε 1}, C_{ε 2}) = \frac{1}{\sqrt{1 - {At}^{2}} | Δ v |} \sqrt{\frac{8 C_{dev} C_{μ} A_{K}^{2}}{A_{ε}}} \\ = (C_{ε 2} - C_{ε 1}) \sqrt{\frac{C_{μ} C_{dev}}{C_{ε 1} C_{ε 2}}}, \end{matrix}

(119)

and growth requires $C_{ε 2} > C_{ε 1}$ ⁠. Note that δ (like α) is degenerate, in that a given value can be obtained using different combinations of the coefficients.

With the density ratio $s \equiv ρ_{1} / ρ_{2} > 1$ ⁠, At $(s) = (s - 1) / (s + 1)$ and

\frac{δ (s)}{δ (1)} = \sqrt{1 - {At}^{2}} = \frac{2 \sqrt{s}}{s + 1},

(120)

which very well describes the reduction in mixing layer growth obtained from the experiments⁶¹ and DNSs of variable-density Kelvin–Helmholtz mixing with density ratios s = 1, 2, 4, and 8.⁶²

As expected from dimensional analysis, $K \propto {(d h / d t)}^{2} \propto {(Δ v)}^{2}, ε \propto {(d h / d t)}^{2} / t \propto {(Δ v)}^{2} / t$ ⁠, S is dimensionless (for a mass fraction variance), $χ \propto 1 / t, ν_{t} \propto hdh / d t \propto {(Δ v)}^{2} t, L \propto h \propto | Δ v | t$ ⁠, and $τ_{m}, τ_{s} \propto t$ ⁠. In particular, ${K, ε, ν_{t}, L} ↓ 0$ as $Δ v ↓ 0$ ⁠, as required in the limit of a stable interface. Using Eqs. (113), (117), and (33), the turbulent lengthscale-to-mixing layer width ratio is

\frac{L (t)}{h (t)} = \frac{C_{ε} π}{4 ξ} \sqrt{\frac{C_{ε 2} - C_{ε 1}}{2 C_{μ} C_{dev} C_{ε 1}}} .

(121)

2. Normalized turbulent kinetic energy and normalized turbulent kinetic energy dissipation rate

The dimensionless centerline turbulent kinetic energy normalized by ${(Δ v)}^{2}$ is

\frac{K (0, t)}{{(Δ v)}^{2}} (C_{ε 1}, C_{ε 2}) = \frac{1}{8} (1 - \frac{C_{ε 1}}{C_{ε 2}}) (1 - {At}^{2}) .

(122)

Integrating Eq. (108) over the mixing layer and using Eqs. (117) and (32), it follows that the integrated turbulent kinetic energy normalized by $| Δ v |^{3}$ is

\begin{matrix} \frac{K_{0}}{| Δ v |^{3}} (C_{μ}, C_{dev}, C_{ε 1}, C_{ε 2}) = \frac{{(C_{ε 2} - C_{ε 1})}^{2}}{12 C_{ε 2}} \sqrt{\frac{C_{μ} C_{dev}}{C_{ε 1} C_{ε 2}}} {(1 - {At}^{2})}^{3 / 2} t \\ = \frac{C_{ε 2} - C_{ε 1}}{12 C_{ε 2}} δ {(1 - {At}^{2})}^{3 / 2} t, \end{matrix}

(123)

which evolves linearly with time in the self-similar state.

Integrating Eq. (109) over the mixing layer and using Eqs. (117) and (32), it follows that the integrated turbulent kinetic energy dissipation rate normalized by ${| Δ v |}^{3}$ is

\begin{matrix} \frac{ε_{0}}{{| Δ v |}^{3}} (C_{μ}, C_{dev}, C_{ε 1}, C_{ε 2}) = \frac{1}{12} (1 - \frac{C_{ε 1}}{C_{ε 2}}) \\ \times \sqrt{\frac{C_{μ} C_{dev} C_{ε 1}}{C_{ε 2}}} {(1 - {At}^{2})}^{3 / 2} . \end{matrix}

(124)

The values of each of these ratios decrease with increasing Atwood number.

Equations (119), (122), (123), and (124) depend on the difference $C_{ε 2} - C_{ε 1}$ ⁠, the magnitude of which is related to the competition between turbulent kinetic energy dissipation rate shear production (i.e., $C_{ε 1}$ ⁠) and destruction (i.e., $C_{ε 2}$ ⁠). Equations (122)–(124) must be positive, again requiring $C_{ε 2} > C_{ε 1}$ ⁠.

3. Molecular mixing parameter

Using Eqs. (65) and (A68), the molecular mixing parameter is

θ_{m} (C_{ε 1}, C_{ε 2}, C_{χ}) = \frac{2 C_{ε 1} C_{χ}}{C_{ε 2} - C_{ε 1} + 2 C_{ε 1} C_{χ}}

(125)

for the three-equation model and must be positive. For specific, fixed values of $C_{ε 1}$ and $C_{ε 2}$ ⁠, θ_m is determined by the value of $C_{χ}$ ⁠.

For the four-equation model,

θ_{m} (C_{ε 1}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = 1 - 4 A_{S}

(126)

with A_S given by Eq. (A77). Expression (125) is considerably simpler than (126) and depends only on three coefficients.

D. Production-to-dissipation/destruction and mechanical-to-scalar timescale ratios

Explicit expressions for the KH self-similar production-to-dissipation/destruction and mechanical-to-scalar timescale ratios corresponding to the two-, three-, and four-equation models are presented here as a function of the model coefficients.

1. Production-to-dissipation/destruction ratios

The turbulent kinetic energy production-to-dissipation ratio is

\frac{P_{K}}{D_{K}} (C_{ε 1}, C_{ε 2}) = \frac{C_{dev} ν_{t} {(\partial \bar{v} / \partial z)}^{2}}{ε} = \frac{C_{ε 2}}{C_{ε 1}},

(127)

and the turbulent kinetic energy dissipation rate production-to-destruction ratio is

\frac{P_{ε}}{D_{ε}} = \frac{C_{ε 1}}{C_{ε 2}} \frac{P_{K}}{D_{K}} = 1,

(128)

so that production is exactly balanced by destruction of the turbulent kinetic energy dissipation rate. It follows from $C_{ε 2} > C_{ε 1}$ that $P_{K} / D_{K} > 1$ ⁠.

For the three-equation model, the scalar variance production-to-dissipation ratio is

\frac{P_{S}}{D_{S}} (C_{ε 1}, C_{ε 2}, C_{χ}) = \frac{ν_{t} {(\partial \tilde{s} / \partial z)}^{2}}{σ_{s} χ} = 1 + \frac{C_{ε 2} - C_{ε 1}}{2 C_{ε 1} C_{χ}} .

(129)

Increasing the value of $C_{χ}$ decreases the value of $P_{S} / D_{S}$ ⁠, which is greater than unity. For the four-equation model,

\frac{P_{S}}{D_{S}} (C_{ε 1}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{1}{8 A_{χ}}

(130)

and the scalar variance dissipation rate production-to-destruction ratio is

\frac{P_{χ}}{D_{χ}} (C_{ε 1}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{1}{C_{χ 2}} [C_{χ 0} \frac{ν_{t} {(\partial \tilde{s} / \partial z)}^{2}}{\sqrt{Sc} χ} + C_{χ 3}] \frac{ε S}{K χ} = \frac{1}{C_{χ 2}} (\frac{C_{χ 0} σ_{s}}{8 \sqrt{Sc} A_{χ}} + C_{χ 3}) \frac{C_{ε 1}}{C_{ε 2} - C_{ε 1}} \frac{A_{S}}{A_{χ}} = 1,

(131)

so that production of the scalar variance dissipation rate is exactly balanced by destruction, with A_S and $A_{χ}$ given by Eqs. (A77) and (A78). These quantities must be positive and are singular if the denominators vanish.

2. Mechanical-to-scalar timescale ratio

For the three-equation model, the mechanical-to-scalar timescale ratio is (15) and for the four-equation model,

R (C_{ε 1}, C_{ε 2}, C_{χ 0}, C_{χ 2}, C_{χ 3}, σ_{s}) = \frac{2 (C_{ε 2} - C_{ε 1})}{C_{ε 1}} \frac{A_{χ}}{A_{S}},

(132)

with A_S and $A_{χ}$ given by Eqs. (A77) and (A78).

VI. CALIBRATION OF MODELS USING PHYSICAL OBSERVABLES

The relationships between the Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing layer parameters and the model coefficients derived in Secs. III–V are used here to calibrate the two-, three-, and four-equation models. It may be useful to calibrate turbulence models to be able to predict the specific values of observables (e.g., in applications where the Reynolds number is very large). The models are calibrated by solving for the coefficient values that are consistent with the values of some chosen set of self-similar physical observables, as discussed below. Similar calibrations have been performed in previous similarity analyses.^8,14,18,19 The number of constraints used must equal the number of model coefficients being calibrated: the calibration process is not unique, as there is freedom to choose which coefficient values are fixed a priori and which constraints will be applied.

A. Physical considerations

Different considerations related to the intended application of a model often dictate which constraints are most relevant and what type of model should be used. Assuming consistent calibrations for each type of model, two-, three-, and four-equation models describe mechanical turbulence the same way and with identical profiles in space and time. If only the bulk mechanical turbulence is to be described, a two-equation model will suffice for which sufficient data are available for coefficient calibration for all three instabilities. If it is also important to describe molecular mixing, then a three-equation model is minimally required: the calibration of the additional coefficient is straightforward but the value is generally different for each instability case.

A four-equation model can predict all of the quantities predicted by two- and three-equation models as well as some set of additional scalar turbulence statistics that are used to calibrate the scalar dissipation rate equation coefficients ${C_{χ q}}$ ⁠. As discussed in Sec. VI E, calibration of a four-equation model for RT, RM, and KH mixing is currently challenging, as little scalar turbulence data have been obtained from experimental and numerical simulation data (that ideally should be self-similar). The scalar self-similar observable expressions are also highly nonlinear in the model coefficients, making the calibration much more difficult than for a three-equation model. The three- and four-equation models are less universal than the two-equation model, as different scalar equation coefficient values must be used for each instability case. In numerical (rather than analytical) applications, solving the least number of turbulence model equations is often desirable for robustness across a wider range of flow conditions.

B. Quantities used for calibration

The sixteen (eleven independent allowing $σ_{s} \neq σ_{ρ}$ and $σ_{s} \neq 1 / C_{dev}$ ⁠) coefficients to be calibrated are

\begin{matrix} σ_{s} = σ_{U} = σ_{K} = σ_{ε}, C_{dev}, σ_{ρ} (2 -, 3 -, 4 - equation) \\ C_{μ}, C_{ε 0}, C_{ε 1}, C_{ε 2} (2 -, 3 -, 4 - equation) \\ σ_{S} = σ_{s}, C_{χ} (3 - equation) \\ σ_{χ} = σ_{s}, C_{χ 0}, C_{χ 2}, C_{χ 3} (4 - equation) \end{matrix}

(133)

and the constant self-similar physical observables available to determine these values are

α and $ρ_{0} K / Δ Φ$ for RT mixing;
θ for RM mixing;
δ, $K (0, t) / {(Δ v)}^{2}$ ⁠, and $ε_{0} / | Δ v |^{3}$ for KH mixing;
$P_{K} / D_{K}$ and $P_{ε} / D_{ε}$ for RT mixing and $P_{K} / D_{K}$ for KH mixing;
$P_{S} / D_{S}$ for RT, RM, KH mixing and a three- or four-equation model;
$P_{χ} / D_{χ}$ for RT, RM, and KH mixing and a four-equation model;
θ_m for RT, RM, KH mixing and a three- or four-equation model; and
R for RT, RM, KH mixing and a three- or four-equation model.

An assumed value of $ρ_{0} K / Δ Φ$ has been used in previous calibrations of turbulent lengthscale models.^14,18,19 The production-to-dissipation/destruction ratios and mechanical-to-scalar timescale ratios have not been previously derived for these mixing cases for potential use in model calibration. A hierarchy of increasingly more complete coefficient calibrations can be constructed by assuming successively fewer assumed coefficient values and incorporating more physical constraints. From the mixing layer growth parameters (60), (85), and (119), it is clear that RT, RM, and KH mixing are associated with the turbulent kinetic energy dissipation rate equation coefficients ${C_{ε 0}, C_{ε 2}}, C_{ε 2}$ ⁠, and ${C_{ε 1}, C_{ε 2}}$ ⁠, respectively. From Eq. (C4), $C_{ε 3} = 2$ ⁠, which is common to all of the calibrations and is relevant to compressible applications. For the four-equation model, the scalar dissipation rate equation model coefficients $C_{χ 0}, C_{χ 2}$ ⁠, and $C_{χ 3}$ are associated with θ_m, R, and the scalar production-to-dissipation/destruction ratios, i.e., three self-similar observable values are needed to uniquely determine the values of these three coefficients.

The steps in the calibration can be summarized as follows to solve the constraint:

equation for θ for $C_{ε 2}$ ⁠;
equations for δ and $K (0, t) / {(Δ v)}^{2}$ for $C_{μ} C_{dev}$ and $C_{ε 1}$ ⁠;
equations for α and $ρ_{0} K / Δ Φ$ for $σ_{ρ}$ and $C_{ε 0}$ ⁠;
equation for θ_m for $C_{χ}$ (three-equation); and
equations for $P_{S} / D_{S}, P_{χ} / D_{χ}$ ⁠, θ_m or R for $C_{χ 0}, C_{χ 2}$ ⁠, and $C_{χ 3}$ ⁠.

For the purposes of model calibration and illustration of the self-similar solutions derived here, the following five physical observable values pertaining to each instability case will be used (universality of these values is not implied):

α = 0.12, \frac{ρ_{0} K}{Δ Φ} = 0.5, θ = 0.3, δ = 0.08, \frac{K (0, t)}{{(Δ v)}^{2}} = 0.035 .

(134)

A wide range of values of α has been reported from experiments and numerical simulations,¹ depending on the specific details of the initial conditions and perhaps other factors. The value $α \approx 2 α_{b} \approx 0.12$ represents a value typical of small Atwood number RT DNSs (At $= 0.01$ ⁠, Sc = 1)⁶³ and gas channel RT experiments using air and helium (At $\approx 0.04$ ⁠, Sc $\approx 0.7$ ⁠).⁶⁴ The value $ρ_{0} K / Δ Φ \approx 0.5$ is typical of values obtained from RT water and gas channel experiments^51,52 and numerical simulations^{23,49,50,65,66} The exponent $θ \approx 0.3$ is consistent with shock tube experimental RM data prior to reshock,^67,68 numerical simulations using narrowband initial perturbations,^55,56 and modal modeling^54,69 The values $δ \approx 0.08$ (often called the “spreading parameter”¹⁰) and $K (0, t) / {(Δ v)}^{2} \approx 0.035$ are consistent with incompressible shear layer (KH) experiments in air^10,70 and DNSs and LESs of these experiments.^62,71–73

Evidence from both experiments and numerical simulations indicates that the asymptotic values of α, θ, and δ depend on the spectral characteristics of the multimode initial conditions (e.g., broadband vs narrowband, laminar vs turbulent, and inflow).^{61,63,73–80} RT mixing experiments using salt/fresh water with Sc $\sim 10^{3}$ showed that Schmidt number effects also influence the observable values.⁸¹ DNSs of RM mixing also showed that θ and other statistics, such as θ_m depend on the initial Reynolds number.⁸² This implies that the mixing layer parameters are also sensitive to the initial conditions (and this is currently poorly understood). It is implicitly assumed in the similarity analysis that mixing layer parameters are functions only of the constant model coefficients and are independent of initial conditions. Consequently, applications of models calibrated using similarity to cases with differing initial conditions require recalibration. Note that a model with highly constrained coefficient values is not necessarily a better model in practice than a less constrained model. Adjusting the model coefficient values does not change the shapes of the assumed spatial mean and turbulent field profiles.

Several possible hierarchical calibrations of the K–ε model equations are described below. The choice of which calibration should be used depends on the intended use of the model, which is briefly discussed for each option. As in previous studies, a value of $C_{μ}$ will be assumed.

C. Two-equation model

Various self-similar and nonself-similar calibration options for the mechanical turbulence equations are developed and discussed here.

1. Calibration of $C_{ε 0}, C_{ε 1}$ , and $C_{ε 2}$ for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing

Consider a model intended to describe prescribed mixing layer growth for all three instabilities. With the reasonable assumption that the values of α, θ, and δ are independent and can, therefore, be used to calibrate these coefficients, Eq. (85) gives

C_{ε 2} (θ) = \frac{3 (1 - θ)}{2 - 3 θ},

(135)

Eq. (119) gives

\begin{matrix} C_{ε 1} (C_{μ} C_{dev}, θ, δ) = \frac{3 (1 - θ)}{2 C_{μ} C_{dev} (2 - 3 θ)} \\ \times (2 C_{μ} C_{dev} + δ^{2} - δ \sqrt{4 C_{μ} C_{dev} + δ^{2}}) \end{matrix},

(136)

and Eq. (60) gives

\begin{matrix} C_{ε 0} (C_{μ}, α, θ, σ_{ρ}) = \frac{3 [α σ_{ρ} (2 - θ) + C_{μ} (1 - θ)]}{C_{μ} (2 - 3 θ)} \\ - \frac{3 \sqrt{2 α σ_{ρ} (C_{μ} + 2 α σ_{ρ})} (2 - θ)}{2 C_{μ} (2 - 3 θ)} . \end{matrix}

(137)

From Eq. (B22), the relationship between the ε equation destruction coefficient and isotropic turbulence decay exponent is

C_{ε 2} (n) = 1 + \frac{1}{n},

(138)

and equating Eqs. (135) and (138) gives the linear relationships between n and θ,

n (θ) = 2 - 3 θ, θ (n) = \frac{2 - n}{3} .

(139)

The value $θ = 0.3$ implies n = 1.10, which is essentially the same value used in previous self-similar analyses of L-based models (which used the slightly smaller value $θ = 0.25$ ⁠).^14,18,19

This calibration is similar to that used for a buoyancy–drag–shear model describing all three instabilities, in which the analogous coefficients were calculated using specific values of α, θ, and δ.⁸³ The nonlinearity of these expressions implies that the coefficient values are quite sensitive to changes in α, θ, or δ. Note the $C_{μ}$ always appears multiplied by C_dev in the expression for $C_{ε 1}$ ⁠: the value of $C_{ε 1}$ is unchanged if the values of both $C_{μ}$ and C_dev are adjusted keeping their product fixed.

The simplest calibration that incorporates the three instabilities considered here is as follows. Assume that

σ_{ρ} = σ_{ϕ_{α}} = \frac{1}{C_{dev}} = 1

(140)

as required by strict self-similarity, taking the standard value C_dev = 1 for the K–ε model. Such a choice is motivated by retaining the standard value of $C_{μ}$ and $σ_{ϕ_{α}} \sim O (1)$ ⁠, which would also allow the model to give good predictions for boundary layer and other canonical turbulent free-shear flows (i.e., mixing layers, jets, and wakes) for which the K–ε model was originally developed and calibrated.^10,84–87 The standard K–ε model coefficient values are $C_{μ} = 0.09, C_{ε 1} = 1.44, C_{ε 2} = 1.92, σ_{K} = 1.00$ ⁠, and $σ_{ε} = 1.30 \neq σ_{K}$ (these standard values of σ_K and $σ_{ε}$ were used in previous K–ε models for RT, RM, and KH mixing.^6,8) The values $σ_{ϕ_{α}} = 1$ were also assumed in the calibration of the K–L model for RT and RM mixing.¹⁴

With α, θ, and δ from (134), it follows that

C_{μ} = 0.087, C_{ε 0} = 0.840, C_{ε 1} = 1.457, C_{ε 2} = 1.909 .

(141)

These values of $C_{ε 1}$ and $C_{ε 2}$ are $\approx 1.2 %$ larger and $\approx 0.6 %$ smaller than their standard values. As this calibration uses three physical observables corresponding to the growth parameters and exponent (and assumed values of $C_{μ}$ and C_dev), it yields a model that predicts the self-similar mixing layer growths for each of the three instabilities, but other quantities are unconstrained, such as $ρ_{0} K / Δ Φ$ and $K (0, t) / {(Δ v)}^{2}$ ⁠: using the values in Eq. (141) gives $ρ_{0} K / Δ Φ = 0.922 σ_{s} / σ_{ρ}$ and $K (0, t) / {(Δ v)}^{2} = 0.03 (1 - {At}^{2})$ ⁠. Calibrations that also predict specific values of these two quantities are described below.

2. Calibration of $C_{ε 0}, C_{ε 1}, C_{ε 2}$ , and $σ_{ρ}$ for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing

Let $Ξ = ρ_{0} K / Δ Φ, ϒ = K (0, t) / {(Δ v)}^{2}$ ⁠, and $σ_{ϕ_{α}} = σ_{s}$ ⁠. The value of $C_{ε 2}$ is given by Eq. (135). The values of $C_{ε 1}$ and $C_{μ} C_{dev}$ are determined by solving Eqs. (119) and (122):

C_{ε 1} (θ, ϒ) = C_{ε 2} (θ) \frac{1 - 8 ϒ - {At}^{2}}{1 - {At}^{2}} = \frac{3 (1 - θ)}{2 - 3 θ} \frac{1 - 8 ϒ - {At}^{2}}{1 - {At}^{2}},

(142)

C_{μ} C_{dev} (δ, ϒ) = \frac{(1 - 8 ϒ - {At}^{2}) δ^{2} (1 - {At}^{2})}{64 ϒ^{2}},

(143)

which depend on the Atwood number (this dependence propagates to other coefficients, as shown below). Note that the K and ε shear production terms $P_{K}^{s} = C_{dev} ν_{t} {(\partial \bar{v} / \partial z)}^{2}$ and $P_{ε}^{s} = C_{ε 1} (ε / K) P_{K}^{s}$ in Eqs. (102) and (103) are proportional to $C_{μ} C_{dev} \propto (1 - 8 ϒ - {At}^{2}) (1 - {At}^{2})$ and $C_{μ} C_{dev} C_{ε 1} \propto {(1 - 8 ϒ - {At}^{2})}^{2}$ ⁠, respectively, which decrease with the increasing Atwood number, thereby reducing instability and turbulence growth. The observables Ξ and $ϒ$ are also very likely Atwood number dependent. Note that for KH mixing, Eqs. (A48) and (A50) require $σ_{ϕ_{α}} = 1 / C_{dev}$ for strict self-similarity. However, it also is possible to relax this constraint such that the resulting model is not strictly self-similar.

The values of $C_{ε 0}$ and $σ_{ρ}$ are determined by solving Eqs. (60) and (64),

\begin{matrix} C_{ε 0} (α, θ, Ξ, σ_{ϕ_{α}}) = C_{ε 2} (θ) (1 - \frac{σ_{ρ}}{σ_{ϕ_{α}}} Ξ) + \frac{3 σ_{ρ} Ξ}{4 σ_{ϕ_{α}}} \\ = \frac{3 (1 - θ)}{2 - 3 θ} (1 - \frac{σ_{ρ}}{σ_{ϕ_{α}}} Ξ) + \frac{3 σ_{ρ} Ξ}{4 σ_{ϕ_{α}}} \\ = \frac{3 (1 - θ)}{2 - 3 θ} [1 - \frac{C_{μ} Ξ^{3}}{8 σ_{ϕ_{α}} α (1 - Ξ)}] \\ + \frac{3 C_{μ} Ξ^{3}}{32 σ_{ϕ_{α}} α (1 - Ξ)}, \end{matrix}

(144)

\begin{matrix} σ_{ρ} (α, Ξ) = \frac{2 C_{μ}}{α} \frac{{(C_{ε 2} - C_{ε 0})}^{2}}{(4 C_{ε 0} - 3) (4 C_{ε 2} - 3)} = \frac{C_{μ} Ξ^{2}}{8 α (1 - Ξ)} . \end{matrix}

(145)

Equations (143) and (145) both involve the ratio $C_{μ} / σ_{ρ} = C_{μ} C_{dev}$ and cannot be satisfied simultaneously, showing that the strictly self-similar assumption $σ_{ρ} = σ_{ϕ_{α}} = 1 / C_{dev}$ is inconsistent for RT and RM mixing. Therefore, assume a value $C_{μ} = 0.087$ in Eq. (143).

Assume that

σ_{ρ} = σ_{ϕ_{α}} \neq \frac{1}{C_{dev}}

(146)

(breaking strict self-similarity). With the values in (134), it follows that

\begin{matrix} C_{μ} = 0.087, C_{ε 0} = 1.330, C_{ε 1} (At) = 1.909 \frac{0.72 - {At}^{2}}{1 - {At}^{2}}, \\ C_{ε 2} = 1.909, σ_{ρ} = σ_{ϕ_{α}} = 0.045, σ_{ϕ_{α}} = 1.475 (KH), \\ C_{dev} (At) = 0.938 (1 - {At}^{2}) (0.72 - {At}^{2}) . \end{matrix}

(147)

For At $= 0, C_{ε 1} = 1.375$ ⁠, and C_dev = 0.678, the values of $C_{ε 1}$ and $C_{ε 2}$ are $\approx 4.5 %$ smaller, and $\approx 0.6 %$ smaller than their standard values. The value of C_dev is much smaller than the value required by the K–2L–a model.¹⁹ This calibration indicates that fully self-similar RT and RM mixing with typical values of α and θ leads to turbulent Schmidt numbers more than a factor $\approx 20$ smaller than those for incompressible shear flows. Note that $C_{ε 0}$ does not depend on At, and positivity of C_dev requires At $< \sqrt{0.72} = 0.849$ (this critical Atwood number is set by the value $1 - 8 ϒ = 0.72$ ⁠, which corresponds to At $= 0$ experimental data).

Instead, assume that

σ_{ρ} \neq σ_{ϕ_{α}} = \frac{1}{C_{dev}}

(148)

[again breaking strict self-similarity, but consistent with Eq. (A50)]. In this case, Eqs. (60) and (64) give

\begin{matrix} C_{ε 0} = \frac{C_{μ} C_{ε 2} + (4 C_{ε 2} - 3) α σ_{ρ}}{C_{μ}} \\ - \frac{(4 C_{ε 2} - 3) \sqrt{2 α σ_{ρ} (C_{μ} + 2 α σ_{ρ})}}{2 C_{μ}}, \end{matrix}

(149)

σ_{ρ} = \frac{4 σ_{ϕ_{α}} (C_{ε 2} - C_{ε 0})}{(4 C_{ε 2} - 3) Ξ}

(150)

together with $C_{ε 1}$ and $C_{μ}$ in Eqs. (142) and (143). Simultaneously solving Eqs. (149) and (150) gives

C_{ε 0} (α, θ, δ, Ξ, ϒ) = \frac{3 [2 α (2 - 3 θ) σ_{ϕ_{α}} + C_{μ} (1 - θ) Ξ]}{(2 - 3 θ) (8 α σ_{ϕ_{α}} + C_{μ} Ξ)} = \frac{384 α (2 - 3 θ) ϒ^{2} + 3 (1 - 8 ϒ - {At}^{2}) δ^{2} (1 - {At}^{2}) (1 - θ) Ξ}{(2 - 3 θ) [512 α ϒ^{2} + δ^{2} Ξ (1 - 8 ϒ - {At}^{2}) (1 - {At}^{2})]},

(151)

σ_{ρ} (α, δ, Ξ, ϒ, σ_{ϕ_{α}}) = \frac{8 α σ_{ϕ_{α}}^{2}}{(8 α σ_{ϕ_{α}} + C_{μ} Ξ) Ξ} = \frac{512 α σ_{ϕ_{α}} ϒ^{2}}{[512 α ϒ^{2} + δ^{2} Ξ (1 - 8 ϒ - {At}^{2}) (1 - {At}^{2})] Ξ} .

(152)

With the values in (134),

C_{μ} = 0.087, C_{ε 0} (At) = \frac{0.909 [0.06209 + 0.00672 (1 - {At}^{2}) (0.72 - {At}^{2})]}{0.07526 + 0.0032 (1 - {At}^{2}) (0.72 - {At}^{2})}, C_{ε 1} (At) = 1.909 \frac{0.72 - {At}^{2}}{1 - {At}^{2}}, C_{ε 2} = 1.909, σ_{ρ} (At) = \frac{0.222}{0.07526 + 0.0032 (1 - {At}^{2}) (0.72 - {At}^{2})}, σ_{ϕ_{α}} = \frac{1}{C_{dev}} = 1.475 .

(153)

For At $= 0, C_{ε 0} = 0.784$ and $σ_{ρ} = 2.86$ ⁠. Note that $1 / σ_{ρ} \propto a + b (1 - {At}^{2}) (0.72 - {At}^{2})$ and $C_{ε 0} / σ_{ρ} \propto c [d + e (1 - {At}^{2}) (0.72 - {At}^{2})]$ (a, b, c, d, and e are constants), so that the K and ε equation buoyancy production terms in Eqs. (43) and (44) $P_{K}^{b} = [g_{0} ν_{t} / (σ_{ρ} ρ_{0})] \partial \bar{ρ} / \partial z$ and $P_{ε}^{b} = C_{ε 0} (ε / K) P_{K}^{b}$ are proportional to $At [a + b (1 - {At}^{2}) (0.72 - {At}^{2})]$ and $At c [d + e (1 - {At}^{2}) (0.72 - {At}^{2})]$ ⁠, respectively, and grow nearly linearly with Atwood number.

Finally, assume that

σ_{ρ} \neq σ_{ϕ_{α}} = 1 \neq \frac{1}{C_{dev}}

(154)

(again breaking strict self-similarity). In this case, Eqs. (60) and (64) together with the values for $C_{ε 1}$ and $C_{μ}$ in Eqs. (142) and (143) give Eq. (149) and

σ_{ρ} (α, δ, Ξ, ϒ) = \frac{8 α}{(8 α + C_{μ} Ξ) Ξ} = \frac{512 C_{dev} α ϒ^{2}}{[512 C_{dev} α ϒ^{2} + δ^{2} Ξ (1 - 8 ϒ - {At}^{2}) (1 - {At}^{2})] Ξ}

(155)

so that

C_{ε 0} (α, θ, δ, Ξ, ϒ) = \frac{3 [2 α (2 - 3 θ) + C_{μ} (1 - θ) Ξ]}{(2 - 3 θ) (8 α + C_{μ} Ξ)} = \frac{384 C_{dev} α (2 - 3 θ) ϒ^{2} + 3 (1 - 8 ϒ - {At}^{2}) (1 - {At}^{2}) δ^{2} (1 - θ) Ξ}{(2 - 3 θ) [512 C_{dev} α ϒ^{2} + δ^{2} Ξ (1 - 8 ϒ - {At}^{2}) (1 - {At}^{2})]} .

(156)

With the values in (134),

\begin{matrix} C_{μ} = 0.087, C_{ε 0} (At) = \frac{0.909 [0.04211 + 0.00672 (1 - {At}^{2}) (0.72 - {At}^{2})]}{0.05104 + 0.0032 (1 - {At}^{2}) (0.72 - {At}^{2})}, C_{ε 1} (At) = 1.909 \frac{0.72 - {At}^{2}}{1 - {At}^{2}}, C_{ε 2} = 1.909, σ_{ρ} (At) = \frac{0.96}{0.48 + 0.03 (1 - {At}^{2}) (0.72 - {At}^{2})}, σ_{ϕ_{α}} = 1, σ_{ϕ_{α}} = 1.475 (KH), C_{dev} = 0.678 . \end{matrix}

(157)

The K and ε equation buoyancy production terms also grow nearly linearly with Atwood number. For At $= 0, C_{ε 0} = 0.800$ and $σ_{ρ} = 1.914$ ⁠: these values are comparable to the values 0.85 and 2.00, respectively, recommended by Llor.⁸ As also shown by Llor,⁸ in the broader context of self-similar variable acceleration RT flows (SSVARTs), in which $g_{0} (t) \propto t^{n}$ (n = 0 corresponds to a constant acceleration), model stability requires that $C_{ε 0} > 1$ ⁠: this again favors (147).

The calibration options above lead to considerably different values of $C_{ε 0}$ and $σ_{ρ}$ ⁠. The calibration option in (147) gives small values of the turbulent Schmidt numbers similar to those obtained in the self-similar calibration of the K–L–a and K–2L–a models, $σ_{ϕ_{α}} = 0.06$ ⁠.^18,19 If it is desired to apply the models solely to RT, RM, and KH mixing, then all four calibration options can be used in principle. If it is desirable for the models to give reasonable predictions for broader applications, such as free-shear flows, then the calibration option (141) should be used. If broader applicability of the model is not a strong consideration, then the option (147) is recommended. Calibration (153) is the only option that utilizes the same values of $σ_{ϕ_{α}}$ for all three instabilities. The At-dependence of several model coefficients arising in the calibration options can be regarded as small Atwood number corrections to the $At = 0$ coefficient expressions, and are not expected to be valid for intermediate or large Atwood number flows. All of the calibrations satisfy the inequalities identified in Secs. III–V that must be satisfied by the coefficients.

D. Three-equation model

For a three-equation model, the calibrations given for a two-equation model above are augmented by Eq. (66), which gives for RT mixing,

C_{χ} (θ_{m, R T}, C_{ε 0}, C_{ε 2}) = \frac{θ_{m, R T}}{1 - θ_{m, R T}} \frac{C_{ε 2} - C_{ε 0}}{4 C_{ε 0} - 3},

(158)

by Eq. (87) which gives for RM mixing

C_{χ} (θ_{m, R M}, C_{ε 2}) = \frac{θ_{m, R M}}{1 - θ_{m, R M}} \frac{2 C_{ε 2} - 3}{6},

(159)

and by Eq. (125) which gives for KH mixing

C_{χ} (θ_{m, K H}, C_{ε 1}, C_{ε 2}) = \frac{θ_{m, K H}}{1 - θ_{m, K H}} \frac{C_{ε 2} - C_{ε 1}}{2 C_{ε 1}}

(160)

with the expressions above for $C_{ε 0}, C_{ε 1}$ ⁠, and $C_{ε 2}$ pertaining to the selected calibration choice. Each expression for $C_{χ} \propto θ_{m} / (1 - θ_{m})$ and depends only on the turbulent kinetic dissipation rate equation coefficients ${C_{ε q}}$ ⁠.

At late times, measurements from experiments and DNSs and ILESs (that are not necessarily in a self-similar state) indicate that $θ_{m} \approx 0.8$ for incompressible RT,^{23,52,63,88,89} RM,^55,90,91 and incompressible At $= 0$ KH^73,92 mixing. In comparison with RT and RM mixing, there is considerably less available data on higher-order scalar turbulence statistics in incompressible, variable-density KH mixing. Most of the experimental studies have focused on the measurements of concentration profiles, chemical product formation, and probability density functions, while most of the numerical simulation studies have focused on the mechanical turbulence (e.g., mixing layer widths, components of the Reynolds stress tensor, energy spectra, scaling of structure functions) and Atwood number effects. A DNS study of passive scalar statistics in a self-similar turbulent shear layer reported⁹³ a peak value of passive scalar fluctuation variance $\bar{ϕ^{' 2}} / {(Δ Φ)}^{2} \approx 0.03$ ⁠, corresponding to $θ_{m, K H} = 1 - 4 \bar{ϕ^{' 2}} / {(Δ Φ)}^{2} \approx 0.88$ ⁠. Later, a DNS and LES study of passive scalar mixing in shear layers reported a late-time scalar variance $\approx 0.055$ (see Figs. 4 and 5 of Ref. 73), corresponding to a value $θ_{m, K H} = 1 - 4 \times 0.055 \approx 0.78$ ⁠. More research is needed to better establish the values of the molecular mixing parameter in KH mixing.

The value of $C_{χ}$ is not universal as it has a different value for each instability case, as determined below. With the molecular mixing parameter value

θ_{m} = θ_{m, R T} = θ_{m, R M} = θ_{m, K H} = 0.80

(161)

for all three instability cases, using the coefficient values from the two-equation model calibrations (147), (153), and (157) gives for RT mixing (158),

C_{χ} = 1.000, C_{χ} = 32.55, C_{χ} = 22.07,

(162)

for RM mixing (159),

C_{χ} = 0.545, C_{χ} = 0.545, C_{χ} = 0.545,

(163)

and for KH mixing (160),

C_{χ} = 0.778, C_{χ} = 0.776, C_{χ} = 0.778,

(164)

respectively.

The corresponding scalar variance decay exponents (B23) predicted using the $C_{χ}$ values above are

n_{s} = 2.20, n_{s} = 71.61, n_{s} = 48.55,

(165)

n_{s} = 1.20, n_{s} = 1.20, n_{s} = 1.20,

(166)

n_{s} = 1.71, n_{s} = 1.71, n_{s} = 1.71,

(167)

respectively. As briefly discussed in Appendix B 4, there is a range of values $n_{s} \approx 1.20$ –2.60, indicating that the first option is the only reasonable option, which gives a range n_s = 1.20–2.20 including all three instabilities: this option also corresponds to values $C_{χ} = 0.545$ –1.000 for all three instabilities.

The implications of the following two approximations are of interest. First, assuming that $θ_{m} = θ_{m, R T} = θ_{m, R M} = θ_{m, K H}$ ⁠, equating the three expressions in Eqs. (158)–(160), and taking $C_{ε 2} = 1.909$ gives

\frac{C_{ε 2} - C_{ε 0}}{4 C_{ε 0} - 3} = \frac{2 C_{ε 2} - 3}{6} = \frac{C_{ε 2} - C_{ε 1}}{2 C_{ε 1}} \approx 0.136,

which then implies that $C_{ε 0} = C_{ε 1} = 1.500$ ⁠, which are both larger than the calibrated coefficient values. Second, using the coefficient values pertaining to the first calibration and taking the average value of $C_{χ}$ for all three instabilities $⟨ C_{χ} ⟩ = (1.000 + 0.545 + 0.778) / 3 = 0.774$ ⁠, it follows that Eqs. (158)–(160) give

θ_{m, R T} \approx 0.76, θ_{m, R M} \approx 0.85, θ_{m, K H} \approx 0.80,

(168)

i.e., the resulting RT and KH values are close to the observed values while the RM value is somewhat larger than the observed values $\approx 0.8$ ⁠. Therefore, assuming the same value

C_{χ} = 0.774

(169)

for all three instabilities may provide a reasonable compromise for universal modeling of all three instabilities.

Finally, several general relationships between the scalar variance production-over-dissipation ratio, scalar turbulent timescale, mechanical turbulent timescale, and molecular mixing parameter can be established as follows. First, the expressions for the molecular mixing parameters (66), (87), and (125) and the expressions for the scalar variance production-to-dissipation ratios (70), (95), and (129) for RT, RM, and KH mixing, respectively, show that

\frac{P_{S}}{D_{S}} ({C_{ε q}}, C_{χ}) = 1 + \frac{θ_{h} τ_{s}}{t} = \frac{1}{θ_{m}}

(170)

for the three-equation models, and the expressions for the scalar variance production-to-dissipation ratios (71), (96), and (130) for RT, RM, and KH mixing, respectively, show that

\frac{P_{S}}{D_{S}} ({C_{ε q}}, {C_{χ q}}, σ_{s}) = \frac{θ_{h}}{8 A_{χ}}

(171)

for the four-equation models. If the algebraic model (9) is used for χ in Eq. (171), then

\begin{matrix} \frac{P_{S}}{D_{S}} ({C_{ε q}}, {C_{χ q}}, σ_{s}) \to \frac{P_{S}}{D_{S}} ({C_{ε q}}, C_{χ}) = \frac{θ_{h}}{8 C_{χ}} \frac{A_{K}}{A_{ε} A_{S}} \\ = \frac{θ_{h}}{2 C_{χ}} \frac{τ_{m}}{(1 - θ_{m}) t} = \frac{θ_{h} τ_{s}}{(1 - θ_{m}) t} . \end{matrix}

(172)

Equation (170) shows that larger values of θ_m corresponds to smaller values of $P_{S} / D_{S}$ and because

θ_{m} = \frac{1}{1 + θ_{h} τ_{s} / t} \leq 1,

(173)

it follows that $P_{S} / D_{S} \geq 1$ ⁠, i.e., the scalar variance production must exceed the dissipation if there is molecular mixing. Therefore, for the three-equation models Eq. (173) gives

\frac{τ_{s}}{t} = \frac{1}{θ_{h}} (\frac{1}{θ_{m}} - 1),

(174)

which reduces to $τ_{s} / t = 1 / (4 θ_{h})$ with $θ_{m} = 0.8$ for all three instabilities, or $τ_{s, R T} / t = 1 / 8, τ_{s, R M} / t = 1 / (4 θ) = 10 / 12$ ⁠, and $τ_{s, K H} / t = 1 / 4$ ⁠. For RM mixing, larger values of θ imply smaller values of $τ_{s, R M} / t$ ⁠. The inequality in Eq. (173) obviously requires that $θ_{h} τ_{s} / t > 0$ or $θ_{h} > 0$ ⁠. With (161), it follows that $P_{S} / D_{S} = 1.25$ for all three instabilities.

E. Four-equation model

A model with more transport equations introduces more coefficients, requiring more constraints to specify them, so a four-equation model can predict more quantities than a three-equation model. Of course, the three- and four-equation models can be cross-calibrated to give identical predictions of the same set of common physical observables for both models. A four-equation model is expected to provide better predictions of transition to turbulence of the scalar fields than a three-equation model. The scalar variance dissipation rate equation coefficients must be recalibrated for other values of Sc.

For the four-equation model $R = 2 A_{K} A_{χ} / (A_{ε} A_{S})$ and $θ_{m} = 1 - 4 A_{S}$ ⁠, so that both quantities are functions of the mechanical and scalar model coefficients. The mechanical-to-scalar timescale ratio can be written as $R = 8 A_{K} A_{χ} / [A_{ε} (1 - θ_{m})]$ ⁠, which depends on all of the scalar model coefficients. It follows that the turbulent scalar field prefactors can be expressed directly in terms of θ_m and R,

A_{S} = \frac{1 - θ_{m}}{4}, A_{χ} = \frac{A_{ε}}{8 A_{K}} (1 - θ_{m}) R,

(175)

so that

\frac{A_{χ}}{A_{S}} = \frac{A_{ε}}{2 A_{K}} R .

(176)

However, these expressions do not provide the dependence of A_S and $A_{χ}$ on the model coefficients.

The four-equation model subtracts the coefficient $C_{χ}$ in the three-equation model and adds the coefficients $C_{χ 0}, C_{χ 2}$ ⁠, and $C_{χ 3}$ ⁠, requiring three constraints for their calibration. The molecular mixing parameters (67), (88), and (126), scalar variance production-to-dissipation ratios (71), (96), and (130), scalar variance dissipation rate production-to-destruction ratios (72), (97), and (131), and mechanical-to-scalar timescale ratios (73), (98), and (132) can provide constraints for the three instability cases. However, these quantities have a highly nonlinear algebraic dependence on the model coefficients and Schmidt number and cannot be solved analytically for the coefficients ${C_{χ q}}$ ⁠. Moreover, measurements from either experiments or simulations of the scalar observables discussed here (except for θ_m for RT and RM mixing) are largely unavailable, particularly in a self-similar state—future research should address this gap.

Consider the two-equation calibrations (147), (153), and (157), and let Sc $= 0.7$ (mixing of gases). In general, the calibration of $C_{χ 0}, C_{χ 2}$ ⁠, and $C_{χ 3}$ requires solving the constraints given by the expressions for the self-similar scalar observables for the coefficients ${C_{χ q}}$ ⁠. Given the values of these observables, together with the values of the previously calibrated mechanical turbulence coefficients, the expressions for these observables can then be solved simultaneously (numerically) for the values of the scalar coefficients. The exponents for decaying scalar turbulence derived in Appendix B 4 are not required for this calibration.

Recognizing the limitations noted above and that there is currently limited data available for these observables corresponding to RT and KH mixing, the numerical calibration procedure will be simplified as follows. The values of some of the scalar coefficients will be assumed and either one or two observable constraints will be used to determine the remaining coefficient(s). For all three instabilities, θ_m will constitute one of the physical observables, and for RT mixing one additional observable will be used (⁠ $P_{S} / D_{S}$ ⁠). Even in this simplified case, the solutions to the equations are not unique (i.e., there are generally multiple roots): solutions for which any coefficient value is negative are inadmissible. In the future, when more scalar mixing data becomes available, it will not be necessary to assume values of any of the coefficients.

In order for a candidate model to also predict isotropic decay of scalar variance and its dissipation rate in homogeneous turbulence, the expression for the decay exponent $n_{s} (C_{ε 2}, C_{χ 2}, C_{χ 3})$ in Appendix B 4 requires that $R_{0} / (C_{ε 2} - 1) > 0$ and $R_{0} / (C_{ε 2} - 1) + 1 > 0$ ⁠, respectively [see Eq. (B13)]. However, it is difficult to obtain numerical solutions to the all of the constraint equations for RT, RM, and KH mixing that satisfy these inequalities and give positive values of ${C_{χ q}}$ that are not unreasonably large.

Comparing the expressions (72), (97), and (131) for the scalar variance production-to-destruction ratios shows that in general

\frac{P_{χ}}{D_{χ}} ({C_{ε q}}, {C_{χ q}}, σ_{s}) = \frac{2}{C_{χ 2} R} (\frac{C_{χ 0} σ_{s} θ_{h}}{8 \sqrt{Sc} A_{χ}} + C_{χ 3}),

(177)

which shows that this quantity is proportional to $1 / R$ ⁠.

1. Rayleigh–Taylor mixing

To simplify the numerical calibration procedure assume a value of $C_{χ 3} = 0.1$ or 1.5. In this case, the coefficients $C_{χ 0}$ and $C_{χ 2}$ are determined by two quantities chosen from ${θ_{m}, R, P_{S} / D_{S}, P_{χ} / D_{χ}}$ ⁠. The observable values that can be used are

\frac{P_{S}}{D_{S}} = 1.24, \frac{P_{χ}}{D_{χ}} = 1.01, R = 1.10

(178)

which are the latest-time values obtained from DNS data^31,34 closely corresponding to a water channel RT experiment.³⁵ The value $P_{S} / D_{S}$ above is also similar to those obtained in Ref. 63. This experimental data were consistent with a slightly larger late-time value of $α \approx 0.14$ ⁠: there is presently no other data available for these values.

Taking θ_m and $P_{S} / D_{S}$ [(161) and (178)] as the scalar observables for RT mixing gives for the calibrations (147), (153), and (157),

\begin{matrix} C_{χ 0} = 12.30, C_{χ 2} = 0.640, C_{χ 3} = 0.100, \\ C_{χ 0} = 16.05, C_{χ 2} = 0.805, C_{χ 3} = 1.500, \\ C_{χ 0} = 18.55, C_{χ 2} = 0.950, C_{χ 3} = 1.500, \end{matrix}

(179)

respectively. Using other values of $C_{χ 3}$ or other pairs of observables, such as ${θ_{m}, R}$ or ${θ_{m}, P_{χ} / D_{χ}}$ ⁠, yields other coefficient values.

2. Richtmyer–Meshkov mixing

Taking θ_m (161) as the scalar observable for RM mixing and taking $C_{χ 3} = 0$ or 1.5 give for the calibrations (147), (153), and (157),

\begin{matrix} C_{χ 0} = 0.499, C_{χ 2} = 4.000, C_{χ 3} = 1.500, \\ C_{χ 0} = 0.030, C_{χ 2} = 1.250, C_{χ 3} = 0, \\ C_{χ 0} = 0.304, C_{χ 2} = 2.000, C_{χ 3} = 0, \end{matrix}

(180)

respectively.

3. Kelvin–Helmholtz mixing

Taking θ_m (161) as the scalar observable for KH mixing and taking $C_{χ 3} = 0$ or 1.5 give for the calibrations (147), (153), and (157),

\begin{matrix} C_{χ 0} = 6.570, C_{χ 2} = 2.500, C_{χ 3} = 1.500, \\ C_{χ 0} = 0.353, C_{χ 2} = 1.000, C_{χ 3} = 0, \\ C_{χ 0} = 0.520, C_{χ 2} = 1.000, C_{χ 3} = 0, \end{matrix}

(181)

respectively.

F. Self-similar parameters predicted by the calibrated models

The predicted values of twenty-two self-similar parameters derived in Secs. III–V can be compared, providing additional insight into the relative merits of the two-, three-, and four-equation models using the calibration values in Table I. For each instability $L (t) / h (t) \propto [C_{ε} π / (4 ξ)] \sqrt{σ_{s} / C_{μ}}$ ⁠. Tables II–IV give the predicted values of self-similar parameters corresponding to each calibrated model for RT, RM, and KH mixing together with the equations corresponding to them.

TABLE I.

Summary of values of constant self-similar quantities used to calibrate the coefficients in the two-, three-, and four-equation models and their corresponding constraining expressions. The first five values are used to calibrate the two-equation mechanical turbulence model, the sixth, seventh, and eighth values are additionally used to calibrate the three- and four-equation models, and the ninth value is used in the calibration of the four-equation model for RT mixing.

	α	$\frac{ρ_{0} K}{Δ Φ}$	θ	δ	$\frac{K (0, t)}{{(Δ v)}^{2}}$	$θ_{m, R T}$	$θ_{m, R M}$	$θ_{m, K H}$	${(\frac{P_{S}}{D_{S}})}_{R T}$
Value	0.12	0.50	0.30	0.08	0.035	0.80	0.80	0.80	1.24
Equation	(60)	(64)	(85)	(119)	(122)	(66) or (67)	(87) or (88)	(125) or (126)	(71)

	α	$\frac{ρ_{0} K}{Δ Φ}$	θ	δ	$\frac{K (0, t)}{{(Δ v)}^{2}}$	$θ_{m, R T}$	$θ_{m, R M}$	$θ_{m, K H}$	${(\frac{P_{S}}{D_{S}})}_{R T}$
Value	0.12	0.50	0.30	0.08	0.035	0.80	0.80	0.80	1.24
Equation	(60)	(64)	(85)	(119)	(122)	(66) or (67)	(87) or (88)	(125) or (126)	(71)

TABLE II.

Summary of values of constant self-similar quantities predicted by the calibrated two-, three-, and four-equation models for RT mixing. The values of α, $ρ_{0} K / Δ Φ$ ⁠, and θ_m are the same for each model and are not shown. For the four-equation model, values are shown for each of the three calibrations. The total mixing layer width is h (i.e., ξ = 2) and $C_{ϵ} = 0.5$ (see Sec. VII A).

	${(\frac{P_{K}}{D_{K}})}_{R T}$	${(\frac{P_{ϵ}}{D_{ϵ}})}_{R T}$	${(\frac{P_{S}}{D_{S}})}_{R T}$	${(\frac{P_{χ}}{D_{χ}})}_{R T}$	R_RT	$\frac{τ_{m, R T}}{t}$	$\frac{τ_{s, R T}}{t}$	${(\frac{L}{h})}_{R T}$
Two-equation	2.00	1.39	⋯	⋯	⋯	0.25	⋯	0.07
Equation	(68)	(69)	⋯	⋯	⋯	(54)	⋯	(61)
Three-equation	2.00	1.39	1.24	⋯	2.00	0.25	0.13	0.07
Equation	(68)	(69)	(70)	⋯	(15)	(54)	(54)	(61)
Four-equation (I)	2.00	1.39	1.24	1.37	2.09	0.25	0.12	0.07
Four-equation (II)	34.09	14.00	1.24	1.29	70.1	8.27	0.12	2.33
Four-equation (III)	23.18	9.71	1.24	1.24	48.8	5.55	0.12	1.57
Equation	(68)	(69)	(71)	(72)	(73)	(54)	(57)	(61)

	${(\frac{P_{K}}{D_{K}})}_{R T}$	${(\frac{P_{ϵ}}{D_{ϵ}})}_{R T}$	${(\frac{P_{S}}{D_{S}})}_{R T}$	${(\frac{P_{χ}}{D_{χ}})}_{R T}$	R_RT	$\frac{τ_{m, R T}}{t}$	$\frac{τ_{s, R T}}{t}$	${(\frac{L}{h})}_{R T}$
Two-equation	2.00	1.39	⋯	⋯	⋯	0.25	⋯	0.07
Equation	(68)	(69)	⋯	⋯	⋯	(54)	⋯	(61)
Three-equation	2.00	1.39	1.24	⋯	2.00	0.25	0.13	0.07
Equation	(68)	(69)	(70)	⋯	(15)	(54)	(54)	(61)
Four-equation (I)	2.00	1.39	1.24	1.37	2.09	0.25	0.12	0.07
Four-equation (II)	34.09	14.00	1.24	1.29	70.1	8.27	0.12	2.33
Four-equation (III)	23.18	9.71	1.24	1.24	48.8	5.55	0.12	1.57
Equation	(68)	(69)	(71)	(72)	(73)	(54)	(57)	(61)

TABLE III.

Summary of values of constant self-similar quantities predicted by the calibrated two-, three-, and four-equation models for RM mixing. The values of θ and θ_m are the same for each model and are not shown. For the four-equation model, values are shown for each of the three calibrations. The total mixing layer width is h (i.e., ξ = 2) and $C_{ϵ} = 1.54$ (see Sec. VII B).

	${(\frac{P_{S}}{D_{S}})}_{R M}$	${(\frac{P_{χ}}{D_{χ}})}_{R M}$	R_RM	$\frac{τ_{m, R M}}{t}$	$\frac{τ_{s, R M}}{t}$	${(\frac{L}{h})}_{R M}$
Two-equation	⋯	⋯	⋯	0.91	⋯	0.16
Equation	⋯	⋯	⋯	(82)	⋯	(86)
Three-equation	1.25	⋯	1.09	0.91	0.83	0.16
Equation	(95)	⋯	(98)	(82)	(82)	(86)
Four-equation (I)	1.25	0.71	1.09	0.91	0.83	0.16
Four-equation (II)	1.24	0.09	1.12	0.91	0.81	0.92
Four-equation (III)	1.25	0.42	1.09	0.91	0.81	0.76
Equation	(96)	(97)	(98)	(82)	(57)	(86)

	${(\frac{P_{S}}{D_{S}})}_{R M}$	${(\frac{P_{χ}}{D_{χ}})}_{R M}$	R_RM	$\frac{τ_{m, R M}}{t}$	$\frac{τ_{s, R M}}{t}$	${(\frac{L}{h})}_{R M}$
Two-equation	⋯	⋯	⋯	0.91	⋯	0.16
Equation	⋯	⋯	⋯	(82)	⋯	(86)
Three-equation	1.25	⋯	1.09	0.91	0.83	0.16
Equation	(95)	⋯	(98)	(82)	(82)	(86)
Four-equation (I)	1.25	0.71	1.09	0.91	0.83	0.16
Four-equation (II)	1.24	0.09	1.12	0.91	0.81	0.92
Four-equation (III)	1.25	0.42	1.09	0.91	0.81	0.76
Equation	(96)	(97)	(98)	(82)	(57)	(86)

TABLE IV.

Summary of values of constant self-similar quantities predicted by the calibrated two-, three-, and four-equation models for KH mixing. The values of δ, $K (0, t) / {(Δ v)}^{2}$ ⁠, and θ_m are the same for each model and are not shown. For the four-equation model, values are shown for each of the three calibrations. The total mixing layer width is h (i.e., ξ = 2) and $C_{ϵ} = 1$ (see Sec. VII C).

	$\frac{K_{0}}{\| Δ v \|^{3} t}$	$\frac{ϵ_{0}}{\| Δ v \|^{3}}$	${(\frac{P_{K}}{D_{K}})}_{K H}$	${(\frac{P_{S}}{D_{S}})}_{K H}$	R_KH	$\frac{τ_{m, K H}}{t}$	$\frac{τ_{s, K H}}{t}$	${(\frac{L}{h})}_{K H}$
Two-equation	0.002	0.005	1.39	–	–	0.39	–	0.71
Equation	(123)	(124)	(127)	(129)	–	(114)	–	(121)
Three-equation	0.002	0.005	1.39	1.25	1.56	0.39	0.25	0.71
Equation	(123)	(124)	(127)	(129)	(15)	(114)	(114)	(121)
Four-equation (I)	0.002	0.005	1.39	1.49	1.55	0.39	0.25	0.71
Four-equation (II)	0.002	0.005	1.39	1.49	1.56	0.39	0.25	0.71
Four-equation (III)	0.002	0.005	1.39	1.49	1.55	0.39	0.25	0.71
Equation	(123)	(124)	(127)	(130)	(132)	(114)	(57)	(121)

	$\frac{K_{0}}{\| Δ v \|^{3} t}$	$\frac{ϵ_{0}}{\| Δ v \|^{3}}$	${(\frac{P_{K}}{D_{K}})}_{K H}$	${(\frac{P_{S}}{D_{S}})}_{K H}$	R_KH	$\frac{τ_{m, K H}}{t}$	$\frac{τ_{s, K H}}{t}$	${(\frac{L}{h})}_{K H}$
Two-equation	0.002	0.005	1.39	–	–	0.39	–	0.71
Equation	(123)	(124)	(127)	(129)	–	(114)	–	(121)
Three-equation	0.002	0.005	1.39	1.25	1.56	0.39	0.25	0.71
Equation	(123)	(124)	(127)	(129)	(15)	(114)	(114)	(121)
Four-equation (I)	0.002	0.005	1.39	1.49	1.55	0.39	0.25	0.71
Four-equation (II)	0.002	0.005	1.39	1.49	1.56	0.39	0.25	0.71
Four-equation (III)	0.002	0.005	1.39	1.49	1.55	0.39	0.25	0.71
Equation	(123)	(124)	(127)	(130)	(132)	(114)	(57)	(121)

1. Rayleigh–Taylor mixing

All of the models predict the calibration values of α, $ρ_{0} K / Δ Φ$ ⁠, θ_m, and $P_{S} / D_{S}$ ⁠. The four-equation model calibrations give significantly different values of $P_{K} / D_{K}, P_{ε} / D_{ε}, P_{χ} / D_{χ}$ ⁠, R, $τ_{m} / t$ ⁠, and L/h, but the same value of $τ_{s} / t$ ⁠: the variability results from the significantly different values of $C_{ε 0}$ and $σ_{ρ}$ ⁠. The first four-equation model calibration predicts a larger value of $P_{χ} / D_{χ}$ than the second and third calibrations [which are larger than that in (178)]. The second and third four-equation calibrations gives unphysically large values of $P_{K} / D_{K}, P_{ε} / D_{ε}$ ⁠, R, and $τ_{m} / t$ ⁠. Future measurements of these quantities can guide which of the three calibrations is best for RT mixing and indicate whether these models can simultaneously predict such a large number of self-similar observables.

DNS data corresponding to RT mixing with At = 0.01 and Sc = 1 exploring self-similar scalings of turbulence statistics and fields⁶³ and with At = $7.5 \times 10^{- 4}$ and Sc = 7 modeling a water channel experiment and exploring turbulent transport mechanisms³⁴ confirmed that the spatial profiles of the molecular mixing parameter, production-to-dissipation/destruction ratios, mechanical and scalar turbulent timescales, and mechanical-to-scalar timescale ratio are approximately constant across the mixing layer (varying most near the edges) at late evolution times approaching a self-similar state.^34,63,94 The profiles of $P_{K} / D_{K} (z, t)$ and $P_{S} / D_{S} (z, t)$ indicated average values $\approx 3.5$ (Fig. 33 of Ref. 63) and $\approx 1.25$ (Fig. 32 of Ref. 63), respectively. Data from a DNS of At $= 0.5$ RT mixing indicated $P_{K} / D_{K} (z, t) \approx P_{K} / D_{K} (t) \approx 2.0$ at the latest time (Fig. 4 of Ref. 94). The latest time values reported in Ref. 34 were $P_{K} / D_{K} (t) \approx 2.00, P_{ε} / D_{ε} (t) \approx 1.10, P_{S} / D_{S} (t) \approx 1.24, P_{χ} / D_{χ} (t) \approx 1.01$ ⁠. This study also indicated that at the latest time $ρ_{0} K / Δ Φ \approx 0.57$ and $θ_{m} \approx 0.6$ (and had not attained asymptotic values), somewhat larger and smaller than the calibration values used here. Note that the value $P_{S} / D_{S} = 1.24$ is nearly identical to the self-similar value 1.25 for $θ_{m} = 0.8$ ⁠. The mechanical-to-scalar timescale ratio (15) had late-time values $R \approx 1.5$ –1.8 (Fig. 13 of Ref. 63) and $R \approx 0.9$ –1.3 (Fig. 20 of Ref. 31) in reasonable agreement with the prediction of the first four-equation model. Ristorcelli and Clark computed the late-time value ${(L / h)}_{R C} = [4 / (π C_{ε})] {(L / h)}_{R T} \approx 0.35$ (Fig. 11 of Ref. 63), indicating a value ${(L / h)}_{R C} \approx 0.18$ predicted by the first four-equation model. These (and other self-similar) values surely depend on the initial conditions and other simulation details (including Schmidt number effects), but clearly favor the first four-equation model calibration.

2. Richtmyer–Meshkov mixing

All of the models predict the calibration values of θ, θ_m, and $P_{S} / D_{S}$ (nearly the same as the RT calibration value) and nearly the same value of R, $τ_{m} / t$ ⁠, and $τ_{s} / t$ ⁠. As in the RT mixing case, the four-equation model calibrations predict significantly different values of $P_{χ} / D_{χ}$ ⁠, all of which are less than one, indicating that destruction exceeds production. The two- and three-equation and the first calibration of the four-equation model predict the same value of L/h, while the other four-equation calibrations predict a value that is approximately five to six times larger. Additional data would again be useful for determining which model and calibration is best for RM mixing.

To date, there has been considerably less analysis of self-similarity in RM mixing than for RT and KH mixing. In particular, there have been no experimental or simulation studies of impulsively or shock-accelerated RM mixing that have considered the quantities discussed here (for either nonself-similar or self-similar conditions). Specifically, measurements of the quantities in Table III are not available. DNS and ILES data for RM mixing with finite initial Reynolds numbers⁸⁰ indicate decay exponents $3 θ - 2 = - 1.25$ and −1.41, respectively, for the turbulent kinetic energy integrated over the mixing layers (89), which are $\approx 14 %$ and $\approx 28 %$ smaller than the model value –1.1.

3. Kelvin–Helmholtz mixing

All of the models predict the calibration values of δ, $K (0, t) / {(Δ v)}^{2}$ ⁠, and θ_m. All of the models predict the same values of $P_{K} / D_{K}, P_{S} / D_{S}, τ_{m} / t, τ_{s} / t$ ⁠, and L/h. All of the models predict $K_{0} / (| Δ v |^{3} t) = 0.002$ and $ε_{0} / | Δ v |^{3} = 0.005$ ⁠. The three- and four-equation models predict the same value of R. Additional scalar data are required to guide model selection for KH mixing.

Thus far, there has been more consideration of self-similar evolution of KH (shear) mixing layers than for RT flows, with both At $= 0$ and At $\neq 0$ ⁠. DNS data closely corresponding to the Bell–Mehta shear layer experiment^62,71,73,75 at late time predicted $K (0, t) / {(Δ v)}^{2} \approx 0.035$ and $τ_{y z} / ρ_{0} \approx - 0.01$ in excellent agreement with the current model predictions. DNS data also indicated that at late times $ε_{0} / | Δ v |^{3} \approx 0.006$ ⁠,^62,72,95 also in very good agreement with the model predictions. The evolution of the normalized integrated turbulent kinetic energy (123) was considered using DNS: estimating the slope of the pertinent curve in Fig. 3 of Ref. 72 gives $\approx (1.8 - 1.2) / (460 - 200) \approx 0.002$ ⁠, in agreement with the model predictions. Estimating the late-time ratio of the pertinent production and dissipation profiles in Fig. 21(e) of Ref. 71 and Fig. 4 of Ref. 72 gives $P_{K} / D_{K} \approx 0.0065 / 0.004 \approx 1.39$ ⁠, which is in excellent agreement with the model predictions.

VII. MODEL APPLICATIONS

The previously calibrated and recommended four-equation model will be applied here to miscible Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz instability-induced mixing cases with Sc $= 0.7$ in order to illustrate the solutions of the model. The space–time evolution of fields is presented in unscaled coordinates to facilitate comparison of their broadening as the layer width grows in time quadratically for RT mixing, as a power-law for RM mixing, and linearly for KH mixing as well as their relative rates of growth or decay in time. In addition to the four primary turbulent fields K, ε, S, and χ, the evolution of a_z, ν_t, and L are also shown: a_z enters the buoyancy production term and is a principal field in K–L–a-type models, ν_t drives the production and diffusion terms through gradient-diffusion closures, and L is also a principal field in turbulent lengthscale-based models. The requirements and challenges involved in comparing self-similar profiles to available RT, RM, and KH experimental and simulation data are discussed in Sec. VIII F.

The scalar variance will correspond to the heavy fluid mass fraction variance, and the model coefficients used are given by Eqs. (147), (179), (180), and (181) (taking At $= 0$ ⁠),

σ_{ρ} = σ_{ϕ_{α}} = 0.045, σ_{ϕ_{α}} = 1.475 (KH), C_{dev} = 0.678, C_{μ} = 0.087, C_{ε 0} = 1.330, C_{ε 1} = 1.375, C_{ε 2} = 1.909, C_{χ 0} = 12.30, C_{χ 2} = 0.640, C_{χ 3} = 0.100 (RT), C_{χ 0} = 0.499, C_{χ 2} = 4.000, C_{χ 3} = 1.500 (RM), C_{χ 0} = 6.570, C_{χ 2} = 2.500, C_{χ 3} = 1.500 (KH),

(182)

corresponding to “Four-equation (I)” in Tables II–IV. It is apparent from the expressions for the turbulent fields in Secs. III–V that they are not valid at t = 0 (the fields are either zero or singular). These expressions will be extrapolated to t = 0 here by introducing a time offset t₀ (an initial characteristic timescale) related to the initial mixing layer width, $h (0) = h_{0} = (ξ / 2) A_{h} t_{0}^{θ_{h}}$ or

t_{0} = {(\frac{2 h_{0}}{ξ A_{h}})}^{1 / θ_{h}},

(183)

and letting $t \to t + t_{0}$ in the expressions for the mixing layer widths and in the power-laws $t^{p_{ϕ_{α}}}$ of the turbulent fields. This allows the fields to be defined at t = 0 and approach their asymptotic self-similar forms for $t ≫ t_{0}$ (with larger values of h₀ giving larger values of t₀). The budgets of the four turbulent transport equations (38), with the terms on the right sides reconstructed using the analytic expressions for the fields corresponding to each instability case given in Secs. III–V, are also shown at the latest evolution time to illustrate the relative magnitudes of the production, dissipation/destruction, and turbulent diffusion terms.

As will be seen later, the magnitudes of the turbulent fields are large within the mixing layer core

z \in [\frac{2 h_{s} (t)}{3 ξ}, \frac{2 h_{b} (t)}{3 ξ}] .

It is assumed that ξ = 2, i.e., h(t) is the total mixing layer width. Common to all three instabilities, the turbulent lengthscale grows proportionally with h(t), with a profile $\sqrt{1 - 4 z^{2} / h {(t)}^{2}}$ ⁠. For each budget, the production-to-dissipation/destruction ratios $P_{ϕ_{α}} / D_{ϕ_{α}}$ (evaluated, for example, at the centerline z = 0) are consistent with the values in Tables II–IV corresponding to “Four-equation (I).” The production and dissipation/destruction terms also have parabolic profiles. The turbulent diffusion terms $T_{ϕ_{α}} = (\partial / \partial z) [(ν_{t} / σ_{ϕ_{α}}) \partial ϕ_{α} / \partial z]$ are positive near the layer boundaries and negative in the core, indicating outward transport from the core toward the layer edges (and integrating to zero). The shear production terms in the mechanical turbulence equations $P_{K}^{s} = - (τ_{z z} / ρ_{0}) \partial \tilde{w} / \partial z$ and $P_{ε}^{s} = C_{ε 1} (ε / K) P_{K}^{s}$ are close to zero at very small Atwood number. The advection terms in all of the equations $A_{ϕ_{α}} = \tilde{w} \partial ϕ_{α} / \partial z$ are also nearly zero at very small Atwood number. The magnitudes of the shear production and advection terms, which are positive and negative over the layer, increase with increasing Atwood number.

A. Rayleigh–Taylor mixing

To illustrate the spatiotemporal evolution of self-similar Rayleigh–Taylor turbulence and mixing, the following parameters are used:

ρ_{0} = 1.00 {g / cm}^{3}, At = 0.05, g_{0} = 981 {cm / s}^{2}

(184)

with a final time t_f = 5 s. These parameters are similar to those used in the simulations by Morgan et al.²³ investigating the approach to self-similarity of RT mixing as a function of mode coupling generations. In the numerical evaluation of the model, it is assumed that

h (t) = α At g_{0} {(t + t_{0})}^{2}

(185)

with $h (0) = h_{0} = 0.01$ cm and $t_{0} = \sqrt{h_{0} / (α At g_{0})} = 0.04$ s. The resulting mixing layer width has a constant, linear-in-time, and quadratic-in-time contribution. The very small value of t₀ allows a very rapid transition to a t² self-similar growth. The growth of the mixing layer width with $α = 0.12$ is shown in Fig. 1.

FIG. 1.

View large Download slide

Time-evolution of the mixing layer width (185) for RT mixing.

1. Evolution of fields

The space–time evolution of the self-similar profiles of $\bar{ρ}$ ⁠, K, ε, S, χ, a_z, ν_t, and L (with $C_{ε} = 0.5$ the approximate average value obtained from RT simulation data⁹⁶) given by Eqs. (35) and (47)–(53) is shown in Figs. 2–9. As time evolves, the linear profile of the mean density across the mixing layer widens with a decreasing slope (asymptotically approaching the fully homogeneously mixed density ρ₀), as $\partial \bar{ρ} / \partial z \propto 1 / h (t)$ decreases with time. The mean heavy fluid mass fraction evolves the same way, but is bounded by $[0, 1]$ instead of by $[ρ_{2}, ρ_{1}]$ as for the mean density. The peak values (at z = 0) of the turbulent field profiles evolve in time consistently with the power-law exponents 2, 1, $0$ ⁠, and –1, respectively. The early-time profiles are narrow, widening in space rapidly as the mixing layer width increases as t². The peak value of S(z, t) remains the same, while that of $χ (z, t)$ decreases with time rapidly at early times (asymptotically decreasing to zero when complete mixing has occurred). The mass flux velocity is negative so that the turbulent kinetic energy buoyancy production term is $- g_{0} a_{z} > 0$ ⁠. The profiles of $ε (z, t)$ and $a_{z} (z, t)$ both evolve linearly in time. The turbulent viscosity profile increases rapidly as t³.

FIG. 2.

View large Download slide

Spatiotemporal evolution of the self-similar mean density $\bar{ρ} (z, t) = ρ_{0} [1 + 2 At z / h (t)]$ for RT mixing.

FIG. 3.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent kinetic energy (47) for RT mixing.

FIG. 4.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent kinetic energy dissipation rate (48) for RT mixing.

FIG. 5.

$FIG. 5. Spatiotemporal evolution of the self-similar mass fraction variance (49) for RT mixing.$

View large Download slide

Spatiotemporal evolution of the self-similar mass fraction variance (49) for RT mixing.

FIG. 6.

$FIG. 6. Spatiotemporal evolution of the self-similar mass fraction variance dissipation rate (50) for RT mixing.$

View large Download slide

Spatiotemporal evolution of the self-similar mass fraction variance dissipation rate (50) for RT mixing.

FIG. 7.

View large Download slide

Spatiotemporal evolution of the self-similar mass flux velocity (51) for RT mixing.

FIG. 8.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent viscosity (52) for RT mixing.

FIG. 9.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent lengthscale (53) for RT mixing.

The linear and inverted parabolic self-similar profiles used for the mean and turbulent fields in the K–L,¹⁴ BHR,¹⁷ K–L–a,¹⁸ and K–2L–a¹⁹ model studies are the same as assumed here. The model is calibrated using the same observable values as in Ref. 18 and therefore predicts the same values a priori. Simulation data consistent with this calibration indicate that the heavy fluid mass fraction is nearly linear and the turbulent kinetic energy normalized by its maximum value is well-approximated by an inverted parabolic profile across the layer (see Fig. 1 of Ref. 18, Fig. 15 of Ref. 23, and Fig. 2 of Ref. 19). Figure 26 of Ref. 63 shows that the scalar variance normalized by its maximum value also has an approximately inverse parabolic profile, and Fig. 36 shows that the turbulent lengthscale [see Eq. (12) but with $C_{ε} = 1$ ] normalized by its maximum value has a profile that is broader than an inverted parabola, more consistent with a $\sqrt{1 - 4 z^{2} / h {(t)}^{2}}$ profile. Figure 5 also shows that the molecular mixing parameter is consistent with $θ_{m} = 1 - 4 S (0, t) = 0.8$ ⁠. As also discussed in Ref. 17, the profiles of K and a_z obtained from gas channel experimental data⁵² are approximately inverse parabolic.

The $t^{p_{ϕ_{α}}}$ time-dependence of the turbulent fields complicates the comparison of the profiles to data, which is the reason why self-normalized (to their maximum values) fields have been compared previously, i.e., comparing $ϕ_{α} (z, t) / (A_{ϕ_{α}} t^{p_{ϕ_{α}}})$ obtained from a model and data [see Eq. (23)]. Also, the coordinate is scaled by h(t) so that the width of the self-similar profile automatically matches data consistent with the mixing layer growth. Such comparisons remove both the dependence on the model coefficients (⁠ $A_{ϕ_{α}}$ ⁠) and on the time (⁠ $t^{p_{ϕ_{α}}}$ ⁠), and thus become comparisons of the profile shapes only.

2. Evolution of transport equation budgets

The budget of each turbulent transport equation (43)–(46) is shown at the final time t_f = 5 s in Figs. 10–13. The buoyancy production terms $P_{K}^{b} = - g_{0} a_{z}$ and $P_{ε}^{b} = C_{ε 0} (ε / K) P_{K}^{b}$ are positive, significantly exceeding the dissipation $D_{K} = ε$ and destruction $D_{ε} = C_{ε 2} ε^{2} / K$ terms in magnitude, consistent with the rapid growth of K and ε. The turbulent diffusion terms in these equations are smaller in magnitude compared to the production and dissipation/destruction terms. The mass fraction variance production term $P_{S} = (2 ν_{t} / σ_{m}) {(\partial {\tilde{m}}_{1} / \partial z)}^{2}$ is balanced by the sum of the mass fraction variance dissipation $D_{S} = 2 χ$ and the turbulent diffusion term such that S(0, t) remains constant. The mass fraction variance dissipation rate production $P_{χ} = (C_{χ 0} ν_{t} / \sqrt{Sc}) (ε / K) {(\partial {\tilde{m}}_{1} / \partial z)}^{2} + C_{χ 3} ε χ / K$ is nearly balanced by the sum of the destruction term $D_{χ} = C_{χ 2} χ^{2} / S$ and the turbulent diffusion term, consistent with the slow decay of χ at late times.

FIG. 10.

View large Download slide

Budget of the self-similar turbulent kinetic energy equation (43) at the final time for RT mixing.

FIG. 11.

View large Download slide

Budget of the self-similar turbulent kinetic energy dissipation rate equation (44) at the final time for RT mixing.

FIG. 12.

$FIG. 12. Budget of the self-similar mass fraction variance equation (45) at the final time for RT mixing.$

View large Download slide

Budget of the self-similar mass fraction variance equation (45) at the final time for RT mixing.

FIG. 13.

$FIG. 13. Budget of the self-similar mass fraction variance dissipation rate equation (46) at the final time for RT mixing.$

View large Download slide

Budget of the self-similar mass fraction variance dissipation rate equation (46) at the final time for RT mixing.

Despite extensive research on RT mixing, there is exceptionally little data available pertaining to turbulent transport equation budgets.^34,94 The relative importance of the terms in the transport equation budget for K in Ref. 94 for At $= 0.5$ and for K, ε, S, and χ in Ref. 34 for RT mixing in the Boussinesq approximation is in generally good agreement with the model predictions shown here.

B. Richtmyer–Meshkov mixing

To illustrate the spatiotemporal evolution of self-similar incompressible impulsively driven Richtmyer–Meshkov turbulence and mixing, the following parameters are used:

ρ_{0} = 0.001 {g / cm}^{3}, {At}^{+} = 0.05, Δ v_{s} = 29160 cm / s,

(186)

with a final time t_f = 5 s. The impulsive acceleration

g_{0} (t) = \frac{Δ v_{s}}{10 t_{*}} e^{- {(t / t_{*})}^{2}}

(187)

(with $t_{*} = 0.1$ s) in the K and ε equations decreases to zero very rapidly after t = 0 and is similar to the acceleration used in a buoyancy–shear–drag turbulence model.⁸³ The value of $Δ v_{s}$ is that used in the multimode simulation of a Ma $= 1.8439$ RM instability,⁵⁵ but with Atwood number 0.05 instead of 0.5. Soulard et al.⁵⁶ investigated the dependence of various RM mixing parameters on the initial large-scale energy spectrum in a self-similar state for At $= 0.05$ ⁠. In the numerical evaluation of the model, it is assumed that

h (t) = h_{0} {(1 + t / t_{0})}^{θ}

(188)

with $h (0) = h_{0} = 0.01$ cm and $t_{0} = θ h_{0} / (Δ v_{s} {At}^{+}) = 2.06 \times 10^{- 6}$ s. The very small value of t₀ allows a very rapid transition to a $t^{0.3}$ self-similar growth. The slow power-law growth in time with exponent $θ = 0.3$ of the mixing layer width is shown in Fig. 14 (see Fig. 4 in the compressible study of Ref. 55 for comparison).

FIG. 14.

View large Download slide

Time-evolution of the mixing layer width (188) for RM mixing.

1. Evolution of fields

The space–time evolution of the self-similar profiles of $\bar{ρ}$ ⁠, K, ε, S, χ, a_z, ν_t, and L (with $C_{ε} = 1.54$ the approximate value obtained from RM DNS data⁸²) given by Eqs. (35) (but with At⁺) and Eqs. (75)–(81) is shown in Figs. 15–22. As time evolves, the linear profile of the mean density across the mixing layer widens with a reduced slope (asymptotically approaching the fully homogeneously mixed density ρ₀), as $\partial \bar{ρ} / \partial z \propto 1 / h (t)$ decreases with time as in RT mixing, but at a significantly reduced rate. The heavy fluid mass fraction evolves the same way, but is bounded by $[0, 1]$ instead of by $[ρ_{2}, ρ_{1}]$ as for the density. The peak values (at z = 0) of the turbulent field profiles evolve in time consistently with the power-law exponents –1.4, –2.4, 0, and –1, respectively. The early-time profiles are narrow, widening in space slowly as the width increases as $t^{0.3}$ ⁠. The peak value of S(z, t) remains the same, while that of $χ (z, t)$ decreases with time (asymptotically decreasing to zero when complete mixing has occurred) as in RT mixing. Figure 18 also shows that the molecular mixing parameter is consistent with $θ_{m} = 1 - 4 S (0, t) = 0.8$ ⁠. Except for the mass fraction variance and turbulent lengthscale, the other turbulent fields decay rapidly at early times and more slowly at late times. The mass flux velocity is very small in magnitude and also negative as in RT mixing. The profile of $a_{z} (z, t)$ evolves as $t^{- 0.7}$ and the profile of $ν_{t} (z, t)$ evolves as $t^{- 0.4}$ ⁠. Both K(z, t) and $ε (z, t)$ decrease after the initial impulse at t = 0 after which the fields decay due to dissipation and destruction.

FIG. 15.

View large Download slide

Spatiotemporal evolution of the self-similar mean density $\bar{ρ} (z, t) = ρ_{0} [1 + 2 {At}^{+} z / h (t)]$ for RM mixing.

FIG. 16.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent kinetic energy (75) for RM mixing.

FIG. 17.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent kinetic energy dissipation rate (76) for RM mixing.

FIG. 18.

$FIG. 18. Spatiotemporal evolution of the self-similar mass fraction variance (77) for RM mixing.$

View large Download slide

Spatiotemporal evolution of the self-similar mass fraction variance (77) for RM mixing.

FIG. 19.

$FIG. 19. Spatiotemporal evolution of the self-similar mass fraction variance dissipation rate (78) for RM mixing.$

View large Download slide

Spatiotemporal evolution of the self-similar mass fraction variance dissipation rate (78) for RM mixing.

FIG. 20.

View large Download slide

Spatiotemporal evolution of the self-similar mass flux velocity (79) for RM mixing.

FIG. 21.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent viscosity (80) for RM mixing.

FIG. 22.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent lengthscale (81) for RM mixing.

Groom and Thornber⁸² also evaluated the normalized scalar dissipation rate (proportional to the mechanical-to-scalar timescale ratio),

\frac{L}{\sqrt{2 K / 3}} \frac{χ}{S} = \sqrt{\frac{3}{2}} C_{ε} \frac{K}{ε} \frac{χ}{S} = \sqrt{\frac{3}{2}} C_{ε} \frac{R}{2} \approx 0.66,

(189)

from their DNS data (which is limited to small Re and is not self-similar), or $R \approx 0.70$ ⁠, which is smaller than the value 1.09 in Table III. Figures 11 and 12 in Ref. 97 show that late-time profiles of the heavy fluid mass fraction are approximately linear, except near the mixing layer edges.

2. Evolution of transport equation budgets

The budget of each turbulent transport equation (43)–(46) is shown at the final time t_f = 5 s in Figs. 23–26. The buoyancy production terms $P_{K}^{b}$ and $P_{ε}^{b}$ are nearly zero except at t = 0. The dissipation D_K and destruction $D_{ε}$ terms are large, consistent with the decay of K and ε across the inhomogeneous mixing layer. The turbulent diffusion terms in these equations are smaller in magnitude than the dissipation and destruction terms. The mass fraction variance production term P_S is balanced by the sum of the dissipation term D_S and the diffusion term such that S(0, t) remains constant. The mass fraction variance dissipation rate destruction term $D_{χ}$ exceeds the production term $P_{χ}$ such that χ decreases with time across the layer as in RT mixing.

FIG. 23.

View large Download slide

Budget of the self-similar turbulent kinetic energy equation (43) with impulsive acceleration (187) at the final time for RM mixing.

FIG. 24.

View large Download slide

Budget of the self-similar turbulent kinetic energy dissipation rate equation (44) with impulsive acceleration (187) at the final time for RM mixing.

FIG. 25.

$FIG. 25. Budget of the self-similar mass fraction variance equation (45) at the final time for RM mixing.$

View large Download slide

Budget of the self-similar mass fraction variance equation (45) at the final time for RM mixing.

FIG. 26.

$FIG. 26. Budget of the self-similar mass fraction variance dissipation rate equation (46) at the final time for RM mixing.$

View large Download slide

Budget of the self-similar mass fraction variance dissipation rate equation (46) at the final time for RM mixing.

Despite extensive experimental and computational research on RM mixing, there is exceptionally little data available pertaining to turbulent transport equation budgets.⁹⁷ The relative importance of the terms in the K transport equation budget shown in Fig. 23 in Ref. 97 for RM mixing with At $= 0.5$ (and computed using ILES) is in generally good qualitative agreement with the model prediction shown here: at the latest time, the evolution of K is dominated by turbulent transport and the pressure–dilatation correlation (which is not modeled here), and by numerical dissipation. The S transport equation budget at the latest time in Fig. 17 in Ref. 97 shows that production dominates and the numerical dissipation is of similar magnitude, which qualitatively agrees with the model prediction shown here.

C. Kelvin–Helmholtz mixing

To illustrate the spatiotemporal evolution of self-similar Kelvin–Helmholtz turbulence and mixing, the following parameters are used:

v_{1} = 1500, v_{2} = 900 cm / s,

(190)

which are the stream velocities in the incompressible Bell–Mehta experiment,⁷⁰ and with a final time t_f = 4.5 s. Thus, $Δ v = 600$ cm/s, At $_{v} = 0.25$ ⁠, and At $= 0$ ⁠. In the numerical evaluation of the model, it is assumed that

h (t) = δ | Δ v | (t + t_{0})

(191)

with $h (0) = h_{0} = 0.01$ cm and $t_{0} = h_{0} / (δ | Δ v |) = 2 \times 10^{- 4}$ s. The resulting mixing layer width has a constant and linear-in-time contribution. The very small value of t₀ allows a very rapid transition to a t self-similar growth. The growth of the mixing layer width with $δ = 0.08$ is shown in Fig. 27.

FIG. 27.

View large Download slide

Time-evolution of the mixing layer width (191) for KH mixing.

1. Evolution of fields

The space–time evolution of the self-similar profiles of $\bar{v}$ ⁠, K, ε, S, χ, ν_t, and L (with $C_{ε} = 1$ assumed in the absence of data) given by Eqs. (105) and (108)–(113) is shown in Figs. 28–34. As time evolves, the linear profile of the mean shear velocity across the mixing layer widens with a reduced slope (asymptotically approaching the average stream velocity v₀), as $\partial \bar{v} / \partial z \propto 1 / h (t)$ decreases with time, analogously with the mean density in RT and RM mixing. The peak values (at z = 0) of the turbulent fields evolve in time consistently with the power-law exponents 0, –1, 0, and –1. The early-time profiles are narrow, widening in space as the width increases as t. The value of S(0, t) remains the same, while $χ (z, t)$ decreases with time (asymptotically decreasing to zero). Figure 31 also shows that the molecular mixing parameter is consistent with $θ_{m} = 1 - 4 S (0, t) = 0.8$ ⁠. The profiles of $ε (z, t)$ and $χ (z, t)$ both decrease at the same rate in time. The turbulent viscosity profile increases as t.

FIG. 28.

View large Download slide

Spatiotemporal evolution of the self-similar mean shear velocity $\bar{v} (z, t) = v_{0} [1 + 2 {At}_{v} z / h (t)]$ for KH mixing.

FIG. 29.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent kinetic energy (108) for KH mixing.

FIG. 30.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent kinetic energy dissipation rate (109) for KH mixing.

FIG. 31.

$FIG. 31. Spatiotemporal evolution of the self-similar mass fraction variance (110) for KH mixing.$

View large Download slide

Spatiotemporal evolution of the self-similar mass fraction variance (110) for KH mixing.

FIG. 32.

$FIG. 32. Spatiotemporal evolution of the self-similar mass fraction variance dissipation rate (111) for KH mixing.$

View large Download slide

Spatiotemporal evolution of the self-similar mass fraction variance dissipation rate (111) for KH mixing.

FIG. 33.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent viscosity (112) for KH mixing.

FIG. 34.

View large Download slide

Spatiotemporal evolution of the self-similar turbulent lengthscale (113) for KH mixing.

The linear and inverted parabolic self-similar profiles used for the mean and turbulent fields in the K–2L–a model¹⁹ studies are the same as assumed here. The model is calibrated using the same observable values as in Ref. 19 and, therefore, predicts the same values a priori. Simulation data consistent with this calibration indicates that the mean shear velocity is nearly linear and the turbulent kinetic energy is well-approximated by an inverted parabolic profile across the layer (see Figs. 5 and 6 of Ref. 19). Experimental data (see Fig. 2 of Ref. 70) and numerical simulation data (see Fig. 13 of Ref. 71, Fig. 5 of Ref. 62, Fig. 7 of Ref. 73, and Fig. 2 of Ref. 95) are consistent with a nearly linear mean shear velocity profile within the mixing layer core. It follows from Fig. 29 that $K (0, t) / {(Δ v)}^{2} = 0.035$ ⁠, consistent with the model calibration. The shear Reynolds stress (115) has a negative parabolic profile, and the centerline normalized shear stress (116) is $τ_{y z} (0, t) / [ρ_{0} {(Δ v)}^{2}] = - 0.01$ ⁠, which is in very good agreement with the peak value in Fig. 4(d) of Ref. 70 (but with the opposite sign convention) and in Fig. 4(d) of Ref. 71.

2. Evolution of transport equation budgets

The budget of each turbulent transport equation (102)–(103) and (45)–(46) is shown at the final time t_f = 4.5 s in Figs. 35–38. The turbulent kinetic energy shear production term $P_{K}^{s} = C_{dev} ν_{t} {(\partial \bar{v} / \partial z)}^{2}$ is larger than the dissipation term D_K. The turbulent kinetic energy shear production term $P_{K}^{s}$ is balanced by the sum of the dissipation $D_{K}$ and the smaller in magnitude turbulent diffusion term. The turbulent kinetic energy dissipation rate shear production term $P_{ε}^{s} = C_{ε 1} (ε / K) P_{K}^{s}$ is balanced by the destruction term $D_{ε}$ ⁠, consistent with Eq. (128) and with the very slow late-time decay of ε. The scalar variance production term P_S is balanced by the sum of the dissipation term D_S and the smaller in magnitude turbulent diffusion term, and the scalar variance dissipation rate production term $P_{χ}$ is exactly balanced by the destruction term $D_{χ}$ such that χ decreases very slowly with time across the layer as in RT and RM mixing. The advection and shear production terms associated with the vertical velocity are close to zero.

FIG. 35.

View large Download slide

Budget of the self-similar turbulent kinetic energy equation (102) at the final time for KH mixing.

FIG. 36.

Self-similar Reynolds-averaged mechanical–scalar turbulence models for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz instability-induced mixing in the small Atwood number limit

I. INTRODUCTION

II. MECHANICAL–SCALAR TURBULENCE MODEL FORMULATIONS AND SIMILARITY ANALYSIS

A. Transport equations

B. Self-similarity assumptions

C. Self-similarity analysis

1. General mixing layer width

2. Self-similar turbulent field profiles

3. Self-similar mean and turbulent fields

III. SELF-SIMILAR SOLUTIONS FOR RAYLEIGH–TAYLOR MIXING

A. Equations for Rayleigh–Taylor mixing

B. Self-similar turbulent fields

1. Two- and three-equation models

2. Four-equation model

C. Self-similar parameters

1. Mixing layer growth parameter

2. Ratio of generated turbulent kinetic energy to released potential energy

3. Molecular mixing parameter

D. Production-to-dissipation/destruction and mechanical-to-scalar timescale ratios

1. Production-to-dissipation/destruction ratios

2. Mechanical-to-scalar timescale ratio

IV. SELF-SIMILAR SOLUTIONS FOR RICHTMYER–MESHKOV MIXING

A. Equations for Richtmyer–Meshkov mixing

B. Self-similar turbulent fields

1. Two- and three-equation models

2. Four-equation model

C. Self-similar parameters

1. Mixing layer growth exponent

2. Molecular mixing parameter

3. Decay exponents of mechanical and scalar turbulent fields

D. Production-to-dissipation/destruction and mechanical-to-scalar timescale ratios

1. Production-to-dissipation/destruction ratios

2. Mechanical-to-scalar timescale ratio

V. SELF-SIMILAR SOLUTIONS FOR KELVIN–HELMHOLTZ MIXING

A. Equations for Kelvin–Helmholtz mixing

B. Self-similar turbulent fields

1. Two- and three-equation models

2. Four-equation model

C. Self-similar parameters

1. Mixing layer growth parameter

2. Normalized turbulent kinetic energy and normalized turbulent kinetic energy dissipation rate

3. Molecular mixing parameter

D. Production-to-dissipation/destruction and mechanical-to-scalar timescale ratios

1. Production-to-dissipation/destruction ratios

2. Mechanical-to-scalar timescale ratio

VI. CALIBRATION OF MODELS USING PHYSICAL OBSERVABLES

A. Physical considerations

B. Quantities used for calibration

C. Two-equation model

1. Calibration of Cε0, Cε1, and Cε2 for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing

2. Calibration of Cε0, Cε1, Cε2, and σρ for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing

D. Three-equation model

E. Four-equation model

1. Rayleigh–Taylor mixing

2. Richtmyer–Meshkov mixing

3. Kelvin–Helmholtz mixing

F. Self-similar parameters predicted by the calibrated models

1. Rayleigh–Taylor mixing

2. Richtmyer–Meshkov mixing

3. Kelvin–Helmholtz mixing

VII. MODEL APPLICATIONS

A. Rayleigh–Taylor mixing

1. Evolution of fields

2. Evolution of transport equation budgets

B. Richtmyer–Meshkov mixing

1. Evolution of fields

2. Evolution of transport equation budgets

C. Kelvin–Helmholtz mixing

1. Evolution of fields

2. Evolution of transport equation budgets

1. Calibration of $C_{ε 0}, C_{ε 1}$ , and $C_{ε 2}$ for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing

2. Calibration of $C_{ε 0}, C_{ε 1}, C_{ε 2}$ , and $σ_{ρ}$ for Rayleigh–Taylor, Richtmyer–Meshkov, and Kelvin–Helmholtz mixing