2.1. Governing equations

To study incompressible mixing layers involving two fluid streams with strongly differing densities, the governing equations are formed by considering the full compressible flow equations for a miscible binary fluid mixture and then obtaining the infinite speed of sound incompressible limit (Livescu Reference Livescu2013). Gravity is not included here, but otherwise the governing equations are identical to those describing variable-density (non-Boussinesq) RT flow, as simulated by Cook & Dimotakis (Reference Cook and Dimotakis2001), Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007) and Wei & Livescu (Reference Wei and Livescu2012). To our knowledge, the present study is the first application of these equations to purely shear-driven variable-density fully turbulent mixing layers.

The equations for the instantaneous variables (with partial derivatives denoted following the comma in the subscript, namely $t$ representing the time variable $t$ and an index $i$ representing the relevant spatial direction $x_i$) are

(2.1)\begin{gather} \rho_{,t} + \left(\rho u_j\right)_{,j} = 0, \end{gather}

(2.2)\begin{gather}\left( \rho u_i\right)_{,t} + \left( \rho u_i u_j \right)_{,j} = -p_{,i} + \tau_{ij,j}, \end{gather}

(2.3)\begin{gather}u_{j,j} = -\mathcal{D} \left(\ln \rho \right)_{,jj}, \end{gather}

where the viscous stress, assumed to be Newtonian, is

(2.4)\begin{equation} \tau_{ij} = \mu \left[u_{i,j} + u_{j,i} - \tfrac{2}{3} u_{k,k} \delta_{ij} \right]. \end{equation}

The governing equations are supplemented by slip boundary conditions in the $y$ direction and periodic boundary conditions in $x$ and $z$ directions.

Equation (2.3) represents the non-zero divergence of velocity that occurs due to the change in volume during mixing (while the flow is incompressible). The Fickian form with diffusion coefficient $\mathcal {D}$ represents the infinite sound speed limit of the full multicomponent diffusion operator (Livescu Reference Livescu2013). Equation (2.3) can be derived from the mixture rule $\rho =1/(Y_1/\rho _1 + Y_2/\rho _2)$ (where $Y_1$ and $Y_2$ are species mass fractions of pure fluids with constant densities $\rho _1$ and $\rho _2$, respectively) and species mass fraction transport equations for each species $(\rho Y_m)_{,t} + (\rho Y_m u_j)_{,j} = \mathcal {D}(\rho Y_{m,j})_{,j}$ (Sandoval Reference Sandoval1995; Cook & Dimotakis Reference Cook and Dimotakis2001; Livescu & Ristorcelli Reference Livescu and Ristorcelli2007). The mixture rule can also be connected to the infinite speed of sound limit of the ideal gas mixture equation of state. Alternately, the same divergence relation can be derived as the infinite sound speed limit of the energy transport equation, which demonstrates the consistency of the variable-density governing equations (Livescu Reference Livescu2013). The dynamic viscosity of mixed fluid obeys a relation analogous to the density: $\mu =1/(Y_1/\mu _1 + Y_2/\mu _2)$, where $\mu _1$ is the viscosity of the pure fluid with density $\rho _1$ and $\mu _2$ is the viscosity of the pure fluid with density $\rho _2$, which ensures a uniform Schmidt number, $Sc=\mu /(\rho D)$, throughout the mixture.

2.2. Notations

Many of the statistics are based on averages, which are indicated by the symbol $\langle \rangle$. Generically, the Reynolds mean of a quantity $q$ is denoted by $\langle q \rangle$ and Reynolds fluctuation is $q' = q - \langle q \rangle$. For simple expressions, the Reynolds mean will also be indicated by an overbar, i.e. $\bar {q}$, which is equal to $\langle q \rangle$. As is typical for compressible flows, Favre averaging is employed for the mean governing equations to account for density variations. The Favre mean of a velocity component, $u_i$, is denoted by $\skew2\tilde {U}_i = \langle \rho u_i \rangle / \langle \rho \rangle$ and the Favre fluctuation is $u_i'' = u_i - \skew2\tilde {U}_i$, in contrast to the Reynolds mean, $\bar {U}_i=\langle u_i \rangle$, and fluctuation, $u_i' = u_i - \bar {U}_i$.

Numerical quantities presented in sections below are obtained from averages computed based on homogeneities present within the flows. Since the flow is periodic and homogeneous in the streamwise and spanwise coordinates $x$ and $z$, area averages are computed across $y$-normal planes. Self-similar statistics will also be considered in which profiles should not change with time (except for noise due to lack of statistical convergence) when the $y$ coordinate is scaled by an appropriate length scale. For these statistics, time averaging is also performed over the self-similar growth duration to improve statistical convergence (§ 5.2). The averages computed to obtain Reynolds ($\bar {U}_i$) and Favre ($\skew2\tilde {U}_i$) averages are $x-z$ area averages only when the statistic is a function of time or not in self-similar coordinates, but time averages are taken of the area averages when self-similar statistics are presented and the same set of notations is used for the averaged quantities.

2.3. Numerical approach

The governing equations (2.1)–(2.4) are solved numerically using a pseudo-spectral scheme for spatial discretization in the periodic (streamwise and spanwise) directions and a compact difference scheme for the inhomogeneous (cross-stream) direction of the flow. The algorithm and code are slightly modified from those employed and described by Wei & Livescu (Reference Wei and Livescu2012), Livescu et al. (Reference Livescu, Ristorcelli, Petersen and Gore2010) and Livescu, Wei & Petersen (Reference Livescu, Wei and Petersen2011) for variable-density RT simulations; the equations solved are the same except non-zero mean streamwise velocity is present in the mixing layer.

The cross-stream (normal) velocities at the lower and upper slip wall boundaries are maintained at zero, and this is consistent with the governing equations for this temporal mixing layer. Averaging the divergence equation (2.3) with diffusivity $\mathcal {D}$ assumed constant, and then omitting the terms of the summed indices that vanish due to the homogeneities present in the flow results in

(2.5)\begin{equation} \langle u_2 \rangle_{,2} = - \mathcal{D} \langle \ln \rho \rangle_{,22}. \end{equation}

Integrating across the $y$ domain, this expression becomes $u_2(y_{max})-u_2(y_{min}) = -\mathcal {D}\{[\ln {\rho }]_{,2}(y_{max}) - [\ln {\rho }]_{,2}(y_{min})\}$. Since density remains constant at the free streams existing at the upper and lower walls, it follows that $u_2(y_{max})-u_2(y_{min}) = 0$. Thus, the variable-density equations are consistent with the boundary conditions $u_2(y_{min}) = u_2(y_{max}) = 0$. This argument also holds for thermally induced single-fluid variable-density mixing layers (for which the governing equations are summarized and contrasted with the present equations in Baltzer & Livescu (Reference Baltzer, Livescu, Livescu, Battaglia and Givi2020) and Livescu (Reference Livescu2020)) if the heat conduction coefficient is constant. More complicated cases such as heat release with chemical reaction necessitates non-zero normal velocity at the boundaries, e.g. Higuera & Moser (Reference Higuera and Moser1994). Spatially developing mixing layers also include streamwise gradients in the streamwise velocity, leading to another term remaining in the left-hand side of the divergence equation, which leads to cross-stream velocities at the upper and lower domain velocities associated with entrainment in even the single-density case.

The third-order accurate variable time stepping Adams–Bashforth–Moulton method is used for time integration, coupled with the usual fractional step method. This is adapted for the pressure equation with variable coefficients due to non-zero velocity divergence associated with the variable-density equations. Fourier representations in the periodic coordinate directions allow the variable coefficient Poisson equation for pressure to reduce to an ordinary differential equation in the inhomogeneous direction. Taking advantage of the structure of the compact derivative, direct solvers can be employed for constant coefficient Poisson equations. The algorithm was initially devised for triply periodic buoyant turbulence simulations by Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007) to provide an exact divergence of momentum and thus avoid degrading the overall order of accuracy. This was an advancement from the algorithm used by Sandoval (Reference Sandoval1995) that required an extrapolation of velocity in time in order to determine the divergence of momentum but could degrade the overall temporal order of accuracy from second order.

The variable coefficient Poisson equation for pressure is decomposed into the form $\boldsymbol {\nabla } p/\rho ^{(n+1)}=\boldsymbol {\nabla } q + \boldsymbol {\nabla } \times {\boldsymbol{A}} + \langle {\boldsymbol{L}} \rangle$, which results in a constant coefficient equation corresponding to the dilatational (curl-free) component, $\boldsymbol {\nabla } q$, and implicit equations for the curl (divergence free), $\boldsymbol {\nabla } \times \boldsymbol{A}$, and mean components. The implicit equations are solved iteratively, using the direct Poisson solvers at each step. Due to the periodic boundary conditions, the mean term $\langle \boldsymbol{L} \rangle$ is non-zero only in $y$ direction. Differences compared to the RT algorithm appear in the mean term for the mixing layer because of the mean flow in the streamwise direction. For the RT case, the mean velocity is zero in both (periodic) horizontal directions, while for the mixing layer case, it is zero only in the (periodic) spanwise direction.

This algorithm avoids introducing additional errors that could affect mass conservation or degrade the accuracy from the time stepping method. The dilatational component of $\boldsymbol {\nabla } p/\rho$ is related to mass conservation, which is enforced to machine precision due to the direct solvers involved. The curl component, $\boldsymbol {\nabla } \times \boldsymbol{A}$, is related to the baroclinic production of vorticity. The iterative procedure is performed until the maximum $x$–$z$ planar-average squared change in $\boldsymbol {\nabla } \times \boldsymbol{A}$ relative to the previous iteration value reduces to 0.01 times the squared value of $\boldsymbol {\nabla } \times \boldsymbol{A}$ averaged within the plane, for each component $\alpha$

(2.6)\begin{align} \max_{\substack{j \in \left\{1, \ldots, N_y\right\}\\ \alpha \in \left\{1, 2, 3\right\}}}{\frac{\displaystyle\sum_{k=1}^{N_z} \sum_{i=1}^{N_x}\left[(\boldsymbol{\nabla} \times \boldsymbol{A})_\alpha^{(n)}(x_i,y_j,z_k)-(\boldsymbol{\nabla} \times \boldsymbol{A})_\alpha^{(n-1)}(x_i,y_j,z_k)\right]^2}{\displaystyle\sum_{k=1}^{N_z} \sum_{i=1}^{N_x}\left[(\boldsymbol{\nabla} \times \boldsymbol{A})_\alpha^{(n)}(x_i,y_j,z_k)\right]^2}}<0.01, \end{align}

where $n$ denotes the iteration number. This tolerance ensures small differences compared to convergence to machine precision. Note that each step of the iterative procedure is based on a direct Poisson solver.

No filtering was used in the simulations, so that the small scales are not affected by numerical artefacts. The spatial resolutions were determined by the requirement that the Kolmogorov scale is well resolved and a series of lower resolution, early time mesh convergence studies. The higher Atwood number cases have more stringent spatial resolution requirements, but for consistency, the same resolution was used for all simulations with Atwood number of $0.75$ or below. Therefore, the lowest Atwood number simulations are over-resolved but should yield very high-quality vorticity and velocity gradient statistics. As described below in the discussion of self-similarity, at late times the peak local dissipation decays linearly with time, so the simulations require the finest resolution during the initial growth stage.

Moin & Mahesh (Reference Moin and Mahesh1998) note that the Kolmogorov length scale is often cited as the smallest scale that needs to be resolved, but suggest that this requirement is more stringent than necessary for reliable first- and second-order statistics. For spectral methods, resolution is often expressed as $k_{max} \eta$, where $\eta$ is the average Kolmogorov length scale $(\nu ^3/\epsilon )^{1/4}$ and $k_{max}=a (2{\rm \pi} /L)$ for a spectral representation of $N$ grid points in a domain of length $L$. The leading coefficient of the $k_{max}$ definition depends on the dealiasing employed, up to a maximum of $N/2$ if no truncation is used. The present simulations calculate the advective terms in skew-symmetric form to reduce the aliasing errors for cubic terms (Blaisdell, Mansour & Reynolds Reference Blaisdell, Mansour and Reynolds1991). In DNS intended to maximize Reynolds number, typical values are $1 \le k_{max}\eta \le 2$ (Gotoh & Yeung Reference Gotoh, Yeung, Davidson, Kaneda and Sreenivasan2013), with $k_{max}\eta \approx 1.5$ typical for adequately resolved DNS of isotropic turbulence (Pope Reference Pope2000; Petersen & Livescu Reference Petersen and Livescu2010). Greater resolution may be required when special attention is focused on certain features, such as fine scale structure associated with stretched spiral vortices in isotropic turbulence that requires $k_{max}\eta \gtrsim 4$ (Horiuti & Fujisawa Reference Horiuti and Fujisawa2008) or the alignment of strain rate and vorticity (Hamlington, Schumacher & Dahm Reference Hamlington, Schumacher and Dahm2008).

In the present mixing layer at negligible Atwood number, $k_{max} \eta$ for the Fourier spectral representation of each homogeneous direction reaches a minimum of $\approx 1.7$ at early times at the centreline (where turbulence is most developed) and continuously increases thereafter. Pantano & Sarkar (Reference Pantano and Sarkar2002) report final values of $k_{max} \eta \approx 1.0$, and they rely on spatial filtering that was shown to produce a relatively small amount of non-physical dissipation to improve stability in their simulations. Resolution can also be quantified in terms of grid spacing relative to the average Kolmogorov scale. Almagro et al. (Reference Almagro, García-Villalba and Flores2017) reported horizontal grid spacing finer than $1.8\eta$ during the self-similar growth, whereas the corresponding values in Pantano & Sarkar (Reference Pantano and Sarkar2002) are $3$–$4\eta$. In the present low $A$ simulation, the horizontal grid spacing (${\rm \Delta} x$ and ${\rm \Delta} z$) peaks at $1.8\eta$ during the early-time transition and reduces to $1.0\eta$ during self-similar growth. Since the mixing layer is inhomogeneous and the Kolmogorov microscales shown above calculated from the dissipation at the peak $y$ position does not account for inhomogeneities in the flow scales, these values merely represent a guideline.

For the present high Atwood number simulations, resolutions can be similarly estimated using the isotropic turbulence formula for $\eta$ that does not address how scales may vary with local density variations. For the present $A=0.75$ simulation, which has the same grid spacing as the $A=0.001$ simulation, $k_{max} \eta$ attains a minimum value of 1.8 at early times and is 3.2 to 3.7 during the self-similar growth (which is similar to the values attained in the $A=0.001$ case). For $A=0.75$, the horizontal grid spacing corresponds to a maximum of $1.8\eta$ at early time and decreases to $1.0\eta$ by the end of self-similar growth. For $A=0.87$, the simulation requires a greater number of grid points for the same physical domain size to maintain numerical stability. The calculated $k_{max} \eta$ reaches a minimum value of 2.7 at early times but remains between 4.4 and 5.3 during the identified self-similar growth interval. The horizontal grid spacing corresponds to a maximum of $1.2\eta$ at early time and decreases to $0.6\eta$ by the end of self-similar growth for $A=0.87$. Nonetheless, these values based on isotropic turbulence $\eta$ are not sensitive to localized steep velocity and density gradients at increased Atwood number that are hypothesized to necessitate greater resolution for numerical stability.

The compact finite difference scheme used for the cross-stream ($y$) direction is sixth-order accurate for both the momentum and pressure equations. The uniform grid spacing is finer (reduced to a factor of 0.8: ${\rm \Delta} _y=0.8{\rm \Delta} _x=0.8{\rm \Delta} _z$) in the inhomogeneous direction, in order to compensate for the lower accuracy relative to the Fourier directions. Modified wavenumber analysis for sixth-order compact difference equations indicates errors in differentiating modes become larger at higher wavenumbers (Petersen & Livescu Reference Petersen and Livescu2010). Since differentiation with the Fourier method is exact up to its highest resolved wavenumber, the Fourier method has no error until the Nyquist frequency. This corresponds to a grid spacing of $2 \eta$ if $k_{max} \eta =1.5$. Requiring the compact difference method to produce less than 25 % error in differentiating a mode with this same wavelength dictates that the grid spacing must be refined relative to that of the spectral method by a factor of 0.8. Note that the vast majority of the energy in the flow is at longer wavelengths that have negligible error, according to the modified wavenumber analysis: the lowest $3/4$ of the wavenumbers have errors of less than 3.5 %.

The pressure determined by the fractional step method restores the velocity field divergence to be consistent with (2.3); however, it represents the average pressure over the time step. To recover the instantaneous pressure for calculating budgets and other statistics, the Poisson equation resulting from obtaining the divergence of (2.2) is computed as a post-processing step after the flow has been advanced in time by the fractional step method. The numerical algorithm has been verified to accurately satisfy the governing equations by comparing the time derivatives calculated for various quantity budgets that appear throughout this paper with the appropriate budget right-hand sides.

2.4. Domain size

The domain lengths in the homogeneous streamwise and spanwise directions $L_x$ and $L_z$ are directly related to the convergence of statistical quantities obtained by planar averaging. In addition, these dimensions potentially affect the sizes of structures that grow within the domain. Convergence can be improved either by enlarging the domain size or by using an ensemble of smaller domain simulations. However, a sufficiently large domain is necessary to achieve correct structure growth and interactions.

Several domain sizes were tested and the final dimensions used were found to have minimal evidence of structure growth restriction compared to smaller sizes. From the perspective of initial KH rollup structures with an assumed streamwise wavelength of the most unstable linear instability mode $\lambda _{ls}$, the present mixing layer domain accommodates $64\lambda _{ls}$ in the streamwise direction. This corresponds to six successive mergers; Vreman, Geurts & Kuerten (Reference Vreman, Geurts and Kuerten1997) found that lengths of $8\lambda _{ls}$ (i.e. three successive mergers) were required to reach reasonable self-similarity. In shear flows, the longest scales are oriented along the streamwise direction. The domain therefore has a $L_x/L_z$ ratio of 4, which was adopted by a number of previous temporal mixing layer simulations (e.g. Rogers & Moser Reference Rogers and Moser1994; O'Brien et al. Reference O'Brien, Urzay, Ihme, Moin and Saghafian2014).

The cross-stream domain size, $L_y$, must also be sufficiently large that the mixing layer evolves freely without the slip walls at the $y$ domain boundaries influencing the growth. A series of simulations with different thicknesses has been performed to ensure the statistics are not influenced by the walls for the self-similar time of interest. The initial interface is positioned so that it is nearer the heavy fluid wall than the light fluid wall in proportion to the Atwood number, since the mean velocity neutral point (interface centre) and the most intense turbulence drift to the light fluid side as the flow develops (§ 5.3). The interface is centred within the domain for the $A=0.001$ case (as this effect is negligible at low density ratios). The domain sizes are summarized in table 1. Although initial momentum thickness $\delta _{m,0}$ (defined below) is somewhat ill-defined for making comparisons, comparing $L_x/\delta _{m,0}$ suggests that the domain lengths are approximately 10 times those of Pantano & Sarkar (Reference Pantano and Sarkar2002) and 3.9 times those of Almagro et al. (Reference Almagro, García-Villalba and Flores2017).

Table 1. Summary of simulation domain parameters. The initial interfaces of both velocity and density are each positioned at $y=0$.

3.1. Mean velocity and density profiles

The initial mean velocity profile that approaches the free-stream velocities of $\pm {\rm \Delta} U/2$ at the $y$ boundaries is specified as

(3.1)\begin{equation} \bar{U}_1(y) = \frac{{\rm \Delta} U}{2} \tanh\left(\frac{y}{2 \delta_{m,0}}\right), \end{equation}

where the momentum thickness $\delta _{m,0}$ specifies the initial thickness of the interface. The hyperbolic tangent profile is commonly used in a wide range of mixing layer simulations, such as Riley, Metcalfe & Orszag (Reference Riley, Metcalfe and Orszag1986), Pantano & Sarkar (Reference Pantano and Sarkar2002), Olson et al. (Reference Olson, Larsson, Lele and Cook2011), O'Brien et al. (Reference O'Brien, Urzay, Ihme, Moin and Saghafian2014) and Almagro et al. (Reference Almagro, García-Villalba and Flores2017).

An initial density profile is prescribed to specify the differing compositions (and thus densities) of the fluid streams. The simulations focus on the simplest case of two separate streams of different velocities and densities meeting at a thin interface, so the initial density profiles are aligned with and of the same thickness as the velocity profiles. Thus, the initial density profile is

(3.2)\begin{equation} \bar{\rho}(y) = \rho_0 + \frac{{\rm \Delta} \rho}{2}\tanh \left(\frac{y}{2\delta_{\rho,0}}\right) \end{equation}

with density profile thickness $\delta _{\rho,0}$ chosen to equal $\delta _{m,0}$. This specification of aligned tanh profiles of density and velocity is similar to the approach of Pantano & Sarkar (Reference Pantano and Sarkar2002) and Almagro et al. (Reference Almagro, García-Villalba and Flores2017), although their density variations were attained by varying the thermodynamic properties for a single fluid. In either approach, the mean density of the lower and upper streams of fluid $\rho _0 = (\rho _1 + \rho _2)/2$ is matched between all of the simulations within the set. The desired Atwood numbers $A$ are then attained by specifying free-stream densities $\rho _1 = \rho _0 - {\rm \Delta} \rho /2$ and $\rho _2 = \rho _0 + {\rm \Delta} \rho /2$, where ${\rm \Delta} \rho = \rho _2 - \rho _1 = 2A\rho _0$. Symmetries present in the temporal mixing layer (but not the spatially developing case) result in the flow behaviours being equivalent whether the negative mean streamwise velocity is associated with the light fluid and the positive velocity is associated with the heavy fluid or vice versa, as also noted by Pantano & Sarkar (Reference Pantano and Sarkar2002). Thus, results from a different profile convention can be compared by selecting coordinates to match density profiles and then changing the sign of the mean streamwise velocity to also match.

3.2. Initial disturbance

Only the velocity field is perturbed relative to the mean profile given above to induce the transition to turbulence. This is appropriate because the velocity field drives the instability and turbulence, as observed in the single-density case; this approach also allows the disturbance to be consistent between Atwood numbers. Different velocity disturbances can produce significantly different growth rates at early times in mixing layers (Fathali etal. Reference Fathali, Meyers, Rubio, Smirnov and Baelmans2008), but the present goal is to quickly establish self-similar growth and minimize long-lived large-scale structures that are uniquely associated with initial disturbances. To roughly resemble physical experiments, the velocity perturbation is confined to a thin (in $y$) region centred at the mean velocity profile interface.

In the present simulations, this is accomplished by generating a random field (filling the full domain) that is divergence free and has a 3-D energy spectrum obeying a Gaussian behaviour at high wavenumbers with $k^4$ behaviour at low wavenumbers as $E(k) = ({k}/{k_0})^4 \,\textrm {e}^{-2({k}/{k_0})^2}$. Here, $k=\sqrt {k_1^2+k_2^2+k_3^2}$ is wavenumber and $k_0$ is the prescribed peak wavenumber; $k_0$ is selected to be $\lambda _{ls}/4$, where $\lambda _{ls}$ is the streamwise wavelength of the least stable mode calculated from temporal linear stability analysis for the base velocity profile ($\lambda _{ls}=28\delta _{m,0}$ for the present set-up). This places much of the energy at small scales to quickly establish turbulent motions. The disturbance spectrum is that used by Pantano & Sarkar (Reference Pantano and Sarkar2002) and the positions of the peak wavelength (relative to the least stable wavelength) are similar. The field is then tapered to a thin interface region by multiplying by the Gaussian profile in $y$ to obey $\langle u_i' u_i' \rangle (y) = A \,\textrm {e}^{-({1}/{2}) ({y}/{\sigma })^2}$, where $\sigma$ is the intensity profile thickness chosen to be $2\delta _m$. This is nearly equivalent to the thickness used in Riley et al. (Reference Riley, Metcalfe and Orszag1986) simulations based on measurements of the intensity profile in a mixing layer experiment and to the thickness used by Pantano & Sarkar (Reference Pantano and Sarkar2002). The peak amplitude $A$ is specified for peak intensity $\langle u_i' u_i' \rangle$ of $0.03{\rm \Delta} U^2$ by prescribing a $0.1{\rm \Delta} U$ root-mean-square fluctuation for each velocity component. This relatively strong disturbance reduces the time to reach self-similar growth. The self-similar value of the streamwise turbulent velocity fluctuation intensity reaches approximately 2.5 times this initial value.

This initial velocity disturbance is similar to those used by Riley et al. (Reference Riley, Metcalfe and Orszag1986) (further described in Riley & Metcalfe Reference Riley and Metcalfe1979) and Pantano & Sarkar (Reference Pantano and Sarkar2002) (further described in Pantano-Rubino Reference Pantano-Rubino2000), but details of the implementations differ. The present approach of multiplying the field by the $y$-intensity profile produces divergence, which is corrected by applying the pressure step of the projection method to the velocity field. This step slightly weakens the intensity of the $u_2$ velocity component. Alternatives exist (e.g. applying the profile to a vorticity field, thereby producing a divergence-free velocity field as in Pantano-Rubino Reference Pantano-Rubino2000), but the present method produces an initial velocity field divergence fully consistent with the variable-density incompressible divergence condition (2.3). A small mean $u_2$ velocity is also produced by this step, which is consistent with the divergence condition (as further explained in appendix A). This mean velocity is concentrated at the interface and decays toward the $y$ boundaries; the magnitude is also very small ($<1\,\%$ of ${\rm \Delta} U$ in all simulations shown).

3.3. Viscosity and diffusivity

Momentum thickness Reynolds number, $Re_{m}={\rm \Delta} U \delta _m/\nu$, can be maximized during the self-similar stage by either growing to a large final thickness $\delta _m$ or having a small viscosity $\nu$. The initial configuration is chosen to maximize the thickness growth so that the fully turbulent state is less affected by the initial disturbance. This is achieved by selecting a relatively small initial momentum thickness and appropriate viscosity such that all scales are well resolved and the initial growth is not overly damped. The fundamental velocity scale ${\rm \Delta} U$ to initialize the simulation is arbitrary and can be scaled out. In consistent units, ${\rm \Delta} U=1$ is prescribed with initial momentum thickness of $0.5$ and viscosity of $0.00625$. This initialization results in a Reynolds number $Re_{m}$ of $80$; however, this value has limited meaning before mixing layer evolution sustains the scales of motion.

The Schmidt number $\mathrm {Sc}=\nu /\mathcal {D}$ is chosen to maintain a constant value of $1$ everywhere as the fluids mix. This is imposed by selecting the same values of kinematic viscosity $\nu =\mu /\rho$ for each of the participating fluids (i.e. $\nu _1=\nu _2$) with constant diffusivity $\mathcal {D}$. The choice of constant kinematic viscosity to maintain constant Schmidt number of $1$ is frequently used in other multi-fluid mixing studies (e.g. Sandoval Reference Sandoval1995; Cook & Dimotakis Reference Cook and Dimotakis2001; Livescu & Ristorcelli Reference Livescu and Ristorcelli2007), although maintaining $Sc=0.7$ (which is typical for gases) is also common (e.g. Olson et al. Reference Olson, Larsson, Lele and Cook2011). Note that the choice of constant $\nu$ implies that $\mu \sim \rho$, whereas with real fluids there is typically a weaker dependence on density such as $\mu \sim \sqrt {\rho }$ (Livescu et al. Reference Livescu, Ristorcelli, Petersen and Gore2010).

4.1. Conservation properties

Integrating the mean density conservation equation (4.4) over the $y$ domain indicates that $\int _{y_1}^{y_2} \bar {\rho }\,\textrm {d} y$ is constant with respect to time (total mass within the domain is conserved). The mean momentum equations (4.5) and (4.6), when similarly integrated over the $y$ domain, show that $\int _{y_1}^{y_2} \bar {\rho }\skew2\tilde {U}_i \,\textrm {d} y$ are also constant with respect to time (total momentum within the domain is conserved), when the remaining terms vanish at the boundaries. This is approximately satisfied for (4.5) and (4.6) throughout the duration of the simulation, since the velocity fluctuations remain at low values near the slip walls, and therefore the advective term and Reynolds stress are negligible at the $y$ domain boundaries, while the mean pressure gradients and viscous stresses have relatively little effect.

4.2. Self-similarity

Another property expected of mixing layers is attaining states of self-similar growth. For the temporal configuration, the statistics are functions only of time and the inhomogeneous $y$ position. Assuming self-similarity and that both mean density and velocity profiles are initially centred at $y=0$ (so that no other length scale is introduced in the problem), the time- and $y$-dependencies are eliminated by introducing a new variable $\eta =y/h$, where for the present purpose, $h$ generically represents a length scale that characterizes the $y$-thickness of the mixing layer and grows with time. Specific choices for defining this thickness are discussed below. The scaled coordinate $\eta$ defined here is separate from the Kolmogorov length scale $\eta$ of § 2.3.

As described in appendix B, the mean mass conservation equation (4.1) and Favre mean streamwise momentum equation (4.2) are satisfied for self-similar growth when the growth rate $\textrm {d} h/\textrm {d} t$ is constant and the mean variables are non-dimensionalized as

(4.7)\begin{gather} \bar{\rho}(y,t) = \rho_0 \hat{\rho} (\eta), \end{gather}

(4.8)\begin{gather}\tilde{U}_1(y,t) = ({\rm \Delta} U) \hat{U}_1 (\eta), \end{gather}

(4.9)\begin{gather}\tilde{U}_2(y,t) = \left(\textrm{d} h/\textrm{d} t \right) \hat{U}_2(\eta), \end{gather}

(4.10)\begin{gather}\tilde{R}_{12}(y,t) = ({\rm \Delta} U) \left( \textrm{d} h/\textrm{d} t\right) \hat{R}_{12}(\eta). \end{gather}

Analysing the resulting self-similar mass conservation and streamwise momentum equations (appendix B) reveals relations between the scaled $y$ positions at which features in the statistical profiles occur. Let $\eta _2$ be defined as the $\eta$ point where the Favre cross-stream velocity inflection point occurs [$\textrm {d}\hat {U}_2/{\textrm {d}\eta }(\eta _2)=0$] and $\eta _{12}$ as the point where Favre shear stress has its inflection [$\textrm {d}\hat {R}_{12}/{\textrm {d}\eta }(\eta _{12})=0$]. Then the self-similar analysis proves that $\eta _{12} < \eta _{2} < 0$. That is, the Reynolds stress peak is located further in the light fluid than the peak of mean cross-stream velocity. This analysis does not determine the position $\eta _1$ of the zero crossing of Favre streamwise velocity [$\hat {U}_1(\eta _1)=0$], but this can be empirically investigated in the simulations.

The above analysis and arguments reach similar conclusions to those presented by Pantano & Sarkar (Reference Pantano and Sarkar2002) after developing the self-similar analysis framework while analysing their variable-density flow. It should be noted that these self-similar equations and results pertain to any variable-density mixing layer that obeys the compressible mass conservation and streamwise momentum equations. Specifying particular cases of the flow (in this case, incompressible binary mixing of species, as opposed to thermodynamic variations or high-speed flow) influences the specific forms of the self-similar quantities (assuming states of self-similar growth are reached).

4.3. Thickness definitions

The thicknesses of the density and streamwise velocity profiles are among the most basic global quantities characterizing mixing layers growth. Though the density and velocity mean profiles initially coincide, they need not grow identically as the flow evolves, so various thickness measurements are defined based on both profiles.

Thickness of a mixing layer is traditionally quantified based on the mean streamwise momentum profile, which has a clear connection to the momentum equation (4.2). Momentum thickness is defined as

(4.11)\begin{align} \delta_m(t) &= \frac{1}{\rho_0 {\rm \Delta} U^2} \int_{-\infty}^{\infty} \bar{\rho} \left[ \tilde{U}_1(y,t) - U_{-} \right] \left[U_{+} - \tilde{U}_1(y,t) \right] \textrm{d} y\nonumber\\ &= \int_{-\infty}^{\infty} \frac{\bar{\rho}}{\rho_0} \left( \frac{1}{4} - \frac{\tilde{U}_1^2}{{\rm \Delta} U^2} \right) \textrm{d} y. \end{align}

As the first form emphasizes, this corresponds to the integral of the product representing deficits relative to free streams, which have streamwise velocities of $U_{-}=-{\rm \Delta} U/2$ and $U_{+}={\rm \Delta} U/2$. An analogous thickness could also be defined on a per-mass basis to depend only upon the mean velocity profile

(4.12)\begin{align} \delta_{m,pm} (t) & = \frac{1}{{\rm \Delta} U^2} \int_{-\infty}^{\infty} \left[ \bar{U}_1(y,t) - U_{-} \right] \left[ U_{+} - \bar{U}_1(y,t) \right] \textrm{d} y\nonumber\\ & = \int_{-\infty}^{\infty} \left( \frac{1}{4} - \frac{\bar{U}_1^2}{{\rm \Delta} U^2} \right) \textrm{d} y. \end{align}

This definition uses Reynolds-averaged streamwise velocity rather than Favre averaged to avoid any explicit dependence on the density field. For single-density mixing layers, (4.12) is commonly given as the definition of the momentum thickness because (4.11) reduces to this when density is constant, although (4.11) is the most formal definition.

Several other quantities also are commonly used to characterize mixing layer thickness based on the mean velocity profile, but these are generally less smooth (i.e. more sensitive to lack of statistical convergence) than the integral thicknesses defined above. These other measurements include lengths based on gradients of profiles. Vorticity thickness is obtained from gradients of the Reynolds mean velocity profiles as

(4.13)\begin{equation} \delta_\omega (t) = \frac{{\rm \Delta} U}{\max(|\textrm{d}\bar{U}_1/\textrm{d} y|)}, \end{equation}

as the vorticity magnitude reduces to $|\textrm {d}\bar {U}/\textrm {d} y|$ in the absence of a mean streamwise gradient in cross-stream velocity. This measure based on only a small portion of the mixing layer (where the mean gradient is steepest) has the potential to produce a misleading representation of the thickness of the layer when significant asymmetries are present.

The distance between positions at which the mean velocity reaches a specific per cent (e.g. 10 %) of the difference ${\rm \Delta} U$ between its free-steam values $U_{-}$ and $U_{+}$ (which are associated with fluids having densities $\rho _1$ and $\rho _2$, respectively)

(4.14)\begin{equation} h_{0.1} (t) = y_{\left[\tilde{U}_1=U_{+}-0.1*{\rm \Delta} U\right]} - y_{\left[\tilde{U}_1=U_{-}+0.1*{\rm \Delta} U\right]} = y_{\left[\tilde{U}_1=0.4*{\rm \Delta} U\right]} - y_{\left[\tilde{U}_1=-0.4*{\rm \Delta} U\right]}. \end{equation}

While momentum thickness and vorticity thickness have been the most commonly used thickness measurements in the historic mixing layer literature, Pope (Reference Pope2000) adopts $h_{0.1}$ in treating planar mixing layers, and it has also been recently used by Schwarzkopf etal. (Reference Schwarzkopf, Livescu, Baltzer, Gore and Ristorcelli2016), for example. For brevity, $h$ will be used herein to indicate $h_{0.1}$. This choice of velocity per cent produces measurements that are smoother and less sensitive to statistical fluctuations than selecting a smaller fraction (e.g. $h_{0.01}$) that would yield thicknesses based on the flow far out in the intermittently turbulent/non-turbulent interface. Favre-averaged velocity is used for $h$, although it could alternatively be based on Reynolds-averaged velocity, as could any of the other thickness quantities. For even the highest Atwood numbers, the effect of averaging type on the calculated thickness is negligible: using Reynolds averages instead of Favre averages for $A=0.87$ produces about $1\,\%$ larger values for $h$ and $5\,\%$ larger values for $\delta _m$. Favre averaging is used for all of the velocity-based thicknesses shown except for $\delta _{m,pm}$ and $\delta _\omega$.

For variable-density mixing layers, similar thicknesses may be defined based instead on the density profiles as

(4.15)\begin{gather} \delta_{\rho} (t) = \frac{1}{{\rm \Delta} \rho^2} \int_{-\infty}^{\infty} \left[ \left(\rho_0 - \frac{{\rm \Delta} \rho}{2}\right) - \bar{\rho}(y,t)\right] \left[ \bar{\rho}(y,t) - \left( \rho_0 + \frac{{\rm \Delta} \rho}{2}\right)\right] \textrm{d} y, \end{gather}

(4.16)\begin{gather}\delta_{\textrm{d}\rho/\textrm{d} y} (t) = \frac{{\rm \Delta} \rho}{\max(|\textrm{d}\rho/\textrm{d} y|)}, \end{gather}

(4.17)\begin{gather}h_{\rho,0.1} (t) = y_{\left[\bar{\rho}=\rho_{2}-0.1*{\rm \Delta} \rho\right]} - y_{\left[\bar{\rho}=\rho_{1}+0.1*{\rm \Delta} \rho\right]} = y_{\left[\bar{\rho}=\rho_{0}-0.4*{\rm \Delta} \rho\right]} - y_{\left[\bar{\rho}=\rho_{0}+0.4*{\rm \Delta} \rho\right]}. \end{gather}

Note that $\delta _\rho$ is equivalent to the width measurement introduced by Youngs (Reference Youngs1991, Reference Youngs2009) and also used by Livescu et al. (Reference Livescu, Ristorcelli, Petersen and Gore2010); it is typically written as $W = \int _{-\infty }^{\infty } F_1 F_2 \,\textrm {d} y$, defined based on the mean volume fractions of each species $F_1 = (\bar {\rho } - \rho _1)/(\rho _2 - \rho _1)$ and $F_2 = (\rho _2 - \bar {\rho })/(\rho _2 - \rho _1)$. Typically, a scaling constant $\beta$ is used with $W$ to approximate bubble height in RT flows as $h_b^*=\beta W$; $\beta$ depends on Atwood number in order to represent asymmetries that develop in RT flow structure as Atwood number increases (Livescu et al. Reference Livescu, Ristorcelli, Petersen and Gore2010; Youngs Reference Youngs2013). For brevity, $h_\rho$ will be used herein to indicate $h_{\rho,0.1}$. As shown below, the mean density profiles develop significant asymmetries at high Atwood numbers, which implies that (4.17) can only accurately represent the layer thickness at very low density ratios.

Additional width quantities commonly used for variable-density flows (particularly RT instabilities, e.g. Livescu et al. (Reference Livescu, Ristorcelli, Gore, Dean, Cabot and Cook2009), Zhou & Cabot (Reference Zhou and Cabot2019)) are also relevant. One such quantity, used by Cook & Dimotakis (Reference Cook and Dimotakis2001), Livescu & Ristorcelli (Reference Livescu and Ristorcelli2008) and Livescu et al. (Reference Livescu, Ristorcelli, Petersen and Gore2010), is $h_{X_\rho } = \int _{-\infty }^{\infty } X_{P}(\bar {\rho }) \,\textrm {d} y$, where $X_{P}$ represents the amount of product in a hypothetical fast reaction between the two species

(4.18)\begin{equation} X_{P}(\rho) = \begin{cases} 2\dfrac{\rho - \rho_1}{\rho_2 - \rho_1}, & \rho\leq \dfrac{\rho_1 + \rho_2}{2}, \\ 2\dfrac{\rho_2 - \rho}{\rho_2 - \rho_1}, & \rho\geq \dfrac{\rho_1 + \rho_2}{2}. \end{cases} \end{equation}

Here, $X_P(\bar {\rho })$ corresponds to the mole fraction of fluid fully mixed to the mean density. Physically, $h_{X_\rho }$ is the thickness of mixed fluid that would result if the two fluids were perfectly homogenized within the mixing layer.

5.1. Time evolution of mean profiles and thickness growth

Area-averaged mean profiles of streamwise velocity and density illustrate the basic properties of the mixing layers’ evolution with respect to time. These profiles are shown for two representative Atwood numbers (almost single density and strongly variable density) in figure 2. These profiles form the basis for the thickness scales defined in § 4.3.

Figure 2. Area-averaged mean profiles throughout the simulation run for $A=0.001$ (a,b) and $A=0.75$ (c,d). For Favre mean streamwise velocity $\skew2\tilde {U}_1$ (a,c) and scaled density (b,d), the cross-stream coordinate is scaled by the initial thickness $h_0$ and the profiles demonstrate the interfaces thickening with time. The lack of symmetry about $y=0$ that develops with increased Atwood number is apparent.

Figure 3 displays the time evolution of thickness by several definitions involving the above profiles. All measurements indicate that simulations for each Atwood number approach linear thickness growth with respect to time at late times. Regardless of the specific thickness definition, thickness growth is retarded with increasing Atwood number. The momentum thickness (a) indicates a strong reduction in growth with Atwood number, whereas the momentum thickness per mass (b) and $h$ (c) quantities both indicate weaker growth reduction, as does vorticity thickness (not shown). The thickness evolutions also highlight that the mixing layers grow to many times their initial thicknesses, as desired to reach self-similar growth.

Figure 3. Mixing layer widths: (a) momentum thickness $\delta _m$, (b) momentum thickness per unit mass $\delta _{m,pm}$ and (c) mean velocity thickness $h$ time evolution for each Atwood number. Lines are coloured by Atwood number: $A=0.001$ (black); $A=0.25$ (green); $A=0.50$ (red); $A=0.75$ (blue); $A=0.87$ (orange).

It should be noted that $\delta _{m,0}$, used for non-dimensionalization, is based on initially aligned profiles at $t=0$, before the shifts of mean streamwise velocity relative to mean density have developed. Thus, $\delta _m$ and $\delta _{m,pm}$ are initially essentially equal but evolve differently as the profile shifts develop. The correspondence between initial $\delta _m$ and other initial length scales is $h_0 = 4.39 \delta _{m,0}$ and $\delta _{\omega,0} = 4 \delta _{m,0}$; similar relations apply to the analogous initial density thicknesses $\delta _\rho$, $h_\rho$, and $\delta _{\textrm {d}\rho /\textrm {d} y}$ as well. However, as the profile shapes evolve in transition and turbulent flow, these relations no longer apply.

To evaluate whether constant values for self-similar temporal growth rates are reached, the time derivatives of thicknesses are shown as functions of time for each Atwood number in figure 4. Thicknesses based on integral measures produce relatively smooth growth rates that in each case asymptote to constant values at late time. Growth rates based on $h$ contain more noise than the rates based on the integral quantities, but applying a Hann filter to smooth the thickness vs. time functions produces the result shown in figure 4(c). These results are also consistent with asymptoting growth rate (though statistical fluctuations are present). For vorticity and density gradient thicknesses, calculating the gradient of a mean profile and then extracting its $y$-maximum makes these measurements more sensitive to noise associated with lack of statistical convergence. The sensitivity of the gradients to small-scale noise dictates that a small amount of spatial smoothing (via a Hann filter) first be applied to the instantaneous mean profiles to remove the finest scales of noise before calculating peak gradients.

Figure 4. Time evolution of mixing layer thickness growth rates based on (a) momentum thickness $\delta _m$, (b) momentum thickness per unit mass $\delta _{m,pm}$ and (c) mean velocity thickness $h$. These correspond to the time derivatives of the thickness evolutions shown in figure 3.

5.2. Determining the time interval of self-similar growth

In addition to constant growth rate, another consequence of self-similar growth is the statistical profiles collapsing when appropriately scaled. For example, the mean streamwise velocity and density profiles would collapse to single curves for all times during self-similar growth when $y$ is scaled by thickness (e.g. $\delta _m$ or $h$). As observed by Rogers & Moser (Reference Rogers and Moser1994), mean velocity profiles are relatively insensitive to deviations from self-similar growth. However, fluctuation intensity profiles generally continue to converge after the mean velocities reach their self-similar profiles. Statistical profiles for many quantities are expected to have constant peak values and thus linearly increasing integral values as thickness grows linearly with time. Directly evaluating the time histories of statistics’ peak values comprises a more stringent test of self-similarity, but evaluating their corresponding integral quantities instead is less sensitive to noise.

One statistic that is meaningful for evaluating self-similar growth is integral of cross-stream velocity fluctuation intensity

(5.1)\begin{equation} \mathcal{V} = \frac{1}{{\rm \Delta} U^2 \delta_m} \int_{-\infty}^{\infty} \langle u_2' u_2' \rangle \,\textrm{d} y. \end{equation}

In earlier simulations emphasizing roll-ups of KH vortex structures and their subsequent mergers, Moser & Rogers (Reference Moser and Rogers1993) showed that large values of $\mathcal {V}$ are associated with these features. Conversely, when Rogers & Moser (Reference Rogers and Moser1994) began a mixing layer simulation from a fully turbulent field, no large values were attained but instead $\mathcal {V}$ slowly increased and then asymptoted to the self-similar value. Attili & Bisetti (Reference Attili and Bisetti2013) examined $\mathcal {V}$ for their spatially developing mixing layer beginning from a thin disturbance (similar to that for the present simulations). It overshot the self-similar growth value when the vortices played an important role at early time, but decreased and asymptoted thereafter as the mixing layer reached a self-similar growth regime. This behaviour is compared to that of the present simulation with negligible Atwood number in figure 5(a). The present simulation produces a much weaker peak in $\mathcal {V}$ than the Attili & Bisetti (Reference Attili and Bisetti2013) simulation. Despite the weaker peak, the present simulation follows similar behaviour of approaching self-similarity after the peak. This behaviour contrasts with the asymptoting from below that appears to occur for the fully turbulent initial condition of Rogers & Moser (Reference Rogers and Moser1994). All of the simulations shown in figure 5 display $\mathcal {V}$ values remaining approximately constant throughout their respective self-similar growth periods, and these values are in good agreement between the simulations. In the present simulations, similar behaviour also occurs at increased Atwood numbers.

Figure 5. Evolutions of $y$-integrated planar-average (a) cross-stream component velocity fluctuation and (b) dissipation, two indicators of self-similar growth. Self-similar growth periods are indicated between the vertical lines; the horizontal axes are scaled to match the beginnings such that the tall left-most vertical line applies to all of the flows. The right vertical lines mark the end of self-similar growth for each flow individually. Note that this scaling is not intended to quantify the relative durations of self-similar growth between simulations. Of the lower horizontal axes, the upper-most corresponds to the present $A=0.001$ mixing layer ($\!$), the middle in (b) only corresponds to the density ratio $s=1$ simulation of Almagro et al. (Reference Almagro, García-Villalba and Flores2017) (- - -) and the lower-most corresponds to the simulation of Rogers & Moser (Reference Rogers and Moser1994) ($\cdots$ $\cdots$). The upper horizontal axis corresponds to the spatially developing simulation of Attili & Bisetti (Reference Attili and Bisetti2013) ($\square$).

An important indication of self-similarity employed by Rogers & Moser (Reference Rogers and Moser1994) is total dissipation of turbulent kinetic energy (TKE), which is planar-averaged dissipation $\varepsilon =- \langle \tau _{ij}' u_{i,j}' \rangle$ (from the TKE budget equation) integrated across the entire mixing layer

(5.2)\begin{equation} \mathcal{E} = \int_{-\infty}^{\infty} \varepsilon \,\textrm{d} y. \end{equation}

The rate at which TKE ultimately is dissipated is set by the large-scale motions that scale (in magnitude) with the velocity difference between streams ${\rm \Delta} U$. Since $\mathcal {E}$ has units of velocity cubed, it can be argued on dimensional grounds that $\mathcal {E}$ scales with ${\rm \Delta} U$ only and therefore is constant with respect to time during self-similar growth (Rogers & Moser Reference Rogers and Moser1994). Unlike the velocity fluctuation intensities, the dissipation peak value does not remain constant with respect to time but instead decays in magnitude proportionally with the mixing layer thickness. Thus, its integral over the increasing width as the mixing layer thickens remains constant.

For the essentially single-density case, the dissipation evolution is compared with those of other mixing layer simulations in figure 5(b). The self-similar growth durations are marked as identified in each corresponding reference. Depending on the route of transition, the peak dissipation may also correspond to an overshoot in dissipation prior to self-similar growth or to part of the self-similar growth regime. The former scenario applies to the simulation of Attili & Bisetti (Reference Attili and Bisetti2013) that begin from a thin disturbance. The latter applies to the simulation of Rogers & Moser (Reference Rogers and Moser1994) that begins from a field containing fully turbulent fluid and slowly approaches the self-similar state from below (in terms of dissipation). Attili & Bisetti (Reference Attili and Bisetti2013) discuss these differences and the role of KH structures in the transition in further detail. The present flow corresponds to the former scenario, beginning from a thin disturbance leading to structures that cause dissipation to overshoot, although this is weaker than in Attili & Bisetti (Reference Attili and Bisetti2013) likely due to the form of the disturbance and the temporally developing nature of the flow.

Compared to the close agreement of self-similar $\mathcal {V}$ value with the other simulations in the literature, there is significantly more variation among the self-similar integrated dissipation values. However, the Attili & Bisetti (Reference Attili and Bisetti2013) mixing layer appears to be asymptoting to a value near that observed in the present $A=0.001$ simulation. The self-similar time interval shown for this present simulation (for which $\mathcal {E}$ is one of the determining considerations) maintains $\mathcal {E}$ to a nearly constant value.

The dimensional argument described above for constant $\mathcal {E}$ in self-similar growth holds for the variable-density mixing layers as well. For variable-density mixing layers, the TKE budget equation terms are often defined to include density (e.g. Livescu et al. Reference Livescu, Ristorcelli, Gore, Dean, Cabot and Cook2009), unlike the typical budgets written for single-fluid incompressible mixing layers (e.g. Rogers & Moser Reference Rogers and Moser1994). Therefore, the integrated dissipation must be divided by density to have the units of $({\rm \Delta} U)^3$. One option is to non-dimensionalize by $\rho _0$, the average of the two streams. However, the most typical treatment is to divide $\varepsilon$ by the mean density $\bar {\rho }$, in analogy to Favre averaging other quantities

(5.3)\begin{equation} \tilde{\mathcal{E}} = \int_{-\infty}^{\infty} \frac{\varepsilon}{\bar{\rho}}\, \textrm{d} y. \end{equation}

Figure 6 demonstrates that $\mathcal {E}$ and $\tilde {\mathcal {E}}$ become constant in self-similar growth for each Atwood number. The values for $\mathcal {E}$ scaled by $\rho _0$ and ${{\rm \Delta} U}^3$ decrease strongly with increasing Atwood number, while $\tilde {\mathcal {E}}$ scaled by ${{\rm \Delta} U}^3$ displays a much weaker dependence.

Figure 6. The time histories of $y$-integrated planar-average dissipation (a) divided by local mean density and (b) divided by the average density of the two streams. Both quantities asymptote to constant values for every Atwood number, which is consistent with self-similar growth. The self-similar time periods are marked by vertical lines in (b) and are based on time histories of dissipation as well as other quantities reaching values that are constant within a specified tolerance. Line colours for each Atwood number are as for figure 3.

While linear growth of thickness and constant integrated dissipation are key indicators of self-similar growth (which have been long been employed, e.g. Rogers & Moser Reference Rogers and Moser1994), comparing additional flow statistics profiles produces further useful indications. This was recognized by Vreman et al. (Reference Vreman, Geurts and Kuerten1997), who determined mixing layer growth to be self-similar when ‘the development of the shear layer thickness is linear in time and profiles of normalized statistical quantities at different times coincide’. The time evolutions of profiles can be evaluated by monitoring the peak values of these statistics or examining their integrals in $y$ divided by the thickness (as with $\mathcal {V}$). This latter approach is less sensitive to statistical variability than the peaks. A number of profile quantities are considered in determining the self-similar growth time interval; integral velocity variances and Reynolds stresses are shown in figure 7(a–c), while additional profiles (e.g. cross-correlations between velocity and density) are considered but not shown for brevity. For each Atwood number, the integral turbulence intensities match very closely with the corresponding integral Favre-averaged Reynolds stresses and are nearly identical for $A=0.5$ and below. Comparing between Atwood numbers, there is a consistent trend to lower intensities with increasing $A$ during transition (when the values peak); during self-similar growth, the trend is weak and easily obscured by statistical variability. The $y$-integrated values shown may conceal some of the complexity in weakly changing profile shapes. For the cross-stream component (figure 7b), the intensity increasing at late time is hypothesized to be associated with the turbulent fluctuations reaching and accumulating near the slip walls to affect the interior of the mixing layer. This is expected to occur soonest for the lowest Atwood numbers because they experience the fastest growth. The self-similar time interval is determined to end before this phenomenon affects the flow.

Figure 7. (a–c) Time histories of $y$-integrated planar-average turbulence intensities $\langle u_i' u_i'\rangle$ ($\!$) and corresponding Favre-average Reynolds stresses (- - -), each divided by Favre mean streamwise velocity thickness $h$ for the (a) streamwise, (b) cross-stream and (c) spanwise components. (d) Time histories of $y$-integrated planar-average density fluctuation intensities. Line colours for each Atwood number are as for figure 3.

Variable-density mixing layers introduce additional quantities to be considered for self-similarity, most importantly the density fluctuation intensity $\langle \rho ' \rho ' \rangle$. The integral values of this planar-mean quantity are shown for each Atwood number in figure 7(d). $\langle \rho ' \rho ' \rangle$ can remain within a tolerance of a constant value later than other statistics and thus determine when the self-similar interval begins. These profiles are related to the mixing of the two streams, which is dependent on how fluid is transported into the cores of the mixing layers. Despite the complex mixing behaviour, the simulations indicate that the density fluctuation intensity profiles for each Atwood number approach a unique self-similar scaled profile that remains approximately constant with respect to time.

The integral $\langle \rho ' \rho ' \rangle$ for $A=0.75$ (blue curve) is suggestive of reaching self-similar growth at particularly late time, with a levelling occurring at earlier time before it again increases and levels off. It appears that the flow configuration changes during the second period of rapid increase. This behaviour is responsible for the late starting time of the self-similar period. This increase in maximum density fluctuation intensity also appears to be associated with a smaller increase in integral cross-stream component velocity fluctuation intensity, as shown in figure 7(b).

In summary, the self-similar periods are determined by seeking constant thickness growth rates, constant values of integrated dissipation, and statistical profiles that remain constant when the cross-stream coordinate is self-similarly scaled. In addition to the velocity intensity profiles, density fluctuation intensity profiles must also be considered for variable-density mixing layers. To identify self-similar growth periods in a consistent manner for all Atwood numbers, these conditions are approximated by requiring that thickness growth rates as well as integrals across the cross-stream domain of dissipation, velocity fluctuation intensity $\langle u_i' u_j' \rangle$ and density fluctuation intensity $\langle \rho '^2 \rangle$ be constant to within a specified threshold. The integrals of fluctuation intensity profiles are scaled by thickness ($h$) to attain constant values (or equivalently are integrated with respect to $y/h$), since the integrals would grow proportional to thickness if the self-similar scaled profiles remain constant. Mean profile convergence is accomplished by ensuring the more sensitive fluctuation intensity profiles are converged. This algorithm is consistently applied by determining the longest time interval that each of the quantities specified above remains within 10 % of any value and then retaining the intersection of these time intervals as the self-similar time interval. The very large simulations produce satisfactory adherence to a relatively stringent set of criteria that must be simultaneously satisfied, as indicated by the self-similar periods marked in figure 6. The self-similar periods for other simulations compared in figure 5 are taken from their respective publications. Due to the effects of differing initial momentum thicknesses (and how they relate to the disturbances), the scaled times $t{\rm \Delta} U/\delta _{m,0}$ (or scaled downstream position $x/\delta _{m,0}$ for the spatial-developing case) in this comparison cannot be meaningfully related between simulations. The significantly smaller domains that were feasible for many previous studies could contribute to the difficulties reported in reaching self-similarity (e.g. Vreman et al. Reference Vreman, Sandham and Luo1996, Reference Vreman, Geurts and Kuerten1997; Pantano & Sarkar Reference Pantano and Sarkar2002). In general, questions remain about the universality of the self-similar state (e.g. Dimotakis & Brown Reference Dimotakis and Brown1976; Rogers & Moser Reference Rogers and Moser1994; Vreman et al. Reference Vreman, Geurts and Kuerten1997). However, the thin and broadband disturbance is intended to reduce idiosyncratic large-scale vortices that persist after transition as a result of the initial condition so the present simulations reach generic self-similar states.

Another consideration relevant to the self-similar growth regime is flow Reynolds number. For the flow statistics to be representative of the fully turbulent mixing in practical applications, the Reynolds numbers must be sufficiently large throughout the averaging time duration. In general, significant changes in mixing behaviour have been observed to occur at a Reynolds number threshold (i.e. the mixing transition, Dimotakis Reference Dimotakis2000). Relevant Reynolds numbers are typically defined using the mixing layer thickness or the Taylor microscale. Both scales continuously grow as the mixing layers thicken with time. According to Dimotakis (Reference Dimotakis2000), general necessary conditions for passing the mixing transition for turbulent flows are that the outer-flow Reynolds number exceeds $Re \approx 1$–$2 \times 10^4$ and that Taylor Reynolds number exceeds $Re_\lambda \approx 100$–$140$. Dimotakis defines the former Reynolds number using a visual thickness scale $\delta _{sh}$ that is used in experiments; it has been estimated as $\delta _{sh} \approx 2 \delta _{\omega }$ for numerical simulations (e.g. Rogers & Moser Reference Rogers and Moser1994). This criterion corresponds to attaining $Re_{\omega } \approx 0.5$–$1 \times 10^4$. Table 2 confirms that this condition is satisfied for the self-similar growth statistical averaging periods. The decrease of $Re_m$ values with Atwood number is a consequence of $\delta _m$ decreasing as the velocity profiles shift into lighter density fluid. This complicates interpreting $Re_m$ in variable-density mixing layers.

Table 2. Reynolds numbers during self-similar growth for the present simulations.

Though Taylor microscale is anisotropic in its most fundamental definition, it is estimated using a relation that strictly only applies to homogeneous isotropic turbulence, $\lambda _g = \sqrt {10 \tilde {k} ({\nu }/{\tilde {\varepsilon }})}$. (Averaging the homogeneous-coordinate components of the fundamental Taylor microscale shows good agreement with this estimate for the present mixing layers.) The velocity scale is also taken as its root-mean-square (r.m.s.) value $u'_{\textit{rms}} \approx \sqrt {(\langle u_1'u_1' \rangle + \langle u_2'u_2' \rangle + \langle u_3'u_3' \rangle )/3} = \sqrt {2k/3}$. Using the turbulent kinetic energy and dissipation at the $y$ position of most intense turbulence, the estimate of Taylor microscale Reynolds number is $Re_\lambda =\tilde {k}\sqrt {\frac {20}{3}({1}/{\tilde {\varepsilon }\nu })}$. Using $u'_\textit{rms}$ produces consistency with the velocity scale used in the $\lambda _g$ definition as well as consistency between the turbulent kinetic energy and dissipation included in turbulent kinetic energy budget (in analogy to isotropic turbulence). Although similar definitions are also used for other relevant flows (e.g. Sekimoto, Dong & Jiménez Reference Sekimoto, Dong and Jiménez2016), the mixing layer literature often uses $\sqrt {2k}$ as the velocity scale (rather than $u'_\textit{rms}=\sqrt {2k/3})$ to form $Re_\lambda =(2k)\lambda _g/\nu =k \sqrt {20/(\varepsilon \nu )}$ (e.g. Pantano & Sarkar Reference Pantano and Sarkar2002; O'Brien et al. Reference O'Brien, Urzay, Ihme, Moin and Saghafian2014; Almagro et al. Reference Almagro, García-Villalba and Flores2017). Renormalized to the present convention, the $Re_\lambda$ range during self-similar growth for the single-density mixing layer of Rogers & Moser (Reference Rogers and Moser1994) is $84$–$99$ and for Almagro et al. (Reference Almagro, García-Villalba and Flores2017) is $81$–$87$, for example. The present simulations generally satisfy the $Re_\lambda \approx 100$ (with $Re_\lambda$ is defined in this way) mixing transition guideline given by Dimotakis (Reference Dimotakis2000) before their self-similar growth periods end. The consistency of the statistics within the self-similar growth periods suggests the turbulence is well developed throughout. The initial condition that produces rapid transition is expected to lead to this state more quickly than the large-scale features that persist through other mixing layers’ transitions.

5.3. Time-averaged self-similar statistical profiles

Figure 8. The times included in the plots correspond to the self-similar growth regimes, for which the determination is explained below (§ 5.2). Figure 8 demonstrates that the time series of mean streamwise velocity and density profiles collapse to single curves when the cross-stream coordinate is scaled by the thickness measurement $h$. Similar collapse is also observed when the cross-stream coordinate is instead scaled by $\delta _m$, $\delta _{m,pm}$, or $\delta _\omega$. While $\delta _m$ was used as the thickness length scale in the discussions above to allow comparison with other studies, scaling statistics in terms of the $h$ scale offers interpretive advantages in variable-density flow. For consistency, $h$ will be used as the thickness scale henceforth, except for when making certain comparisons with other studies. The collapse of mean profiles is one indication that self-similar growth is achieved. During self-similar growth, it is thus appropriate to time average the scaled profiles to improve statistical convergence. This averaging is also applied to all of the other scaled statistics presented below.

Figure 8. Mean profiles during self-similar growth for $A=0.001$ (a–d) and $A=0.75$ (e–h). For Favre mean streamwise velocity $\skew2\tilde {U}_1$ (a,b,e,f) and scaled density (c,d,g,h), scaling the cross-stream coordinate by the initial thickness $h_0$ shows only the self-similar time interval curves from figure 2 and demonstrates the growth of the mixing layer (a,c,e,g), whereas scaling instead by each curve's thickness $h$ causes this same time series to collapse to a single curve for each simulation (b,d,f,h).

Comparing the self-similar scaled profiles among Atwood numbers (figure 9) illustrates several basic changes that occur as the density difference between streams increases. For $A=0.001$, the mean streamwise velocity and mean density profiles are essentially centred at $y=0$ and symmetric about that point. A shift in the mean streamwise velocity profiles to the light fluid side (i.e. $\eta _1<0$) that increases in magnitude with increasing Atwood number is apparent. With increasing Atwood number, the shapes of these velocity profiles remain generally similar as they shift to the light fluid side, but the asymmetry in their gradients (figure 9b) reveals an additional steepening on the light fluid side and shallowing on the heavy fluid side. Conversely, the neutral points of the density profiles (where $\bar {\rho }(y)=\rho _0$) remain relatively stationary while the density profiles steepen on the heavy fluid side but become shallower on the light fluid side with increasing Atwood number.

Figure 9. Self-similar time-averaged profiles for all Atwood numbers showing (a) mean streamwise velocity, (b) its $y$ gradient and (c) scaled mean density. In (a) and (b), the solid line represents Reynolds mean, while the dashed line represents Favre mean. Lines are coloured by Atwood number as in figure 3.

Figure 10 displays the corresponding profiles for the cross-stream mean velocity component. The magnitudes are much smaller than those of the streamwise velocity. However, as the self-similar analysis indicates, the cross-stream velocity has an important relationship with mass conservation and mixing layer growth in variable-density mixing layers. In figure 10(b,c), these velocity profiles are shown with the scaling suggested by the self-similar analysis, using $h$ based on the Favre mean streamwise velocity for the thickness scaling. The Reynolds-averaged cross-stream velocity is much smaller in magnitude than the corresponding Favre-average quantity. In addition, $V$ can be shown to strongly depend on the mean density gradient (appendix A) and therefore not reach a time-constant magnitude during self-similar growth; the averages in figure 10(c) should be understood to pertain only to their particular averaging time periods. It is shown below (§ 5.6) that $\tilde {V}$ is dominated by the turbulent mass flux, which does approach a constant value during self-similar growth.

Figure 10. Self-similar time-averaged profiles of cross-stream velocity for all Atwood numbers. Solid lines represent Reynolds mean, while the dashed lines represent Favre mean. Line colours are by Atwood number as in figure 3. The velocities are scaled using (a) mean streamwise velocity and (b,c) growth rate as suggested by the self-similar analysis (§ 4.2); (c) is scaled to show the Reynolds mean, which is negligible compared to the Favre mean in (b); the Reynolds means should be understood to pertain only to their respective averaging times, since they do not become constant in time. No scaling by Atwood number is applied, so the magnitudes for the $A=0.001$ (black line) are small in comparison to the other cases and appear near $\skew2\tilde {U}_2/{\rm \Delta} U=0$ on this vertical scale.

The positions of the neutral points (i.e. $\rho =\rho _0$ and $\skew2\tilde {U}_1=0$) and positions of extrema for various statistical quantities (e.g. $\min (\skew2\tilde {U}_2)$) are important in characterizing the shape of the mixing layer during the self-similar regime. The mixing layer growth and its asymmetry can be summarized by tracking the points at which the mean streamwise velocity is equal to 10 % and 90 % of the free-stream difference ${\rm \Delta} U$: $y_{[\skew2\tilde {U}_1=U_{-}+0.1*{\rm \Delta} U]}$ and $y_{[\skew2\tilde {U}_1=U_{+}-0.1*{\rm \Delta} U]}$. These are the points whose separation define $h$ in (4.14). In figure 11(a), the linear growth of these positions (scaled by initial thickness) with respect to time, approximately extending from $y=0$ at $t=0$, is consistent with the positions collapsing to fixed self-similar scaled (e.g. $y/h$) values. The $\skew2\tilde {U}_1$-based positions also evolve linearly and likewise collapse to fixed $y/h$ values. Plotting the scaled positions of these points as a function of Atwood number (figure 11b) highlights the prominent features observed in figure 9: an increasing drift of the mean streamwise velocity profile to the light fluid side with increasing Atwood number, while the density profile remains approximately centred at the initial interface. In addition, figure 11(b) indicates that the mean cross-stream velocity $\skew2\tilde {U}_2$ peak similarly drifts to the light fluid side, as well as the peak Reynolds stress $\tilde {R}_{12}$ (§ 5.4). The relative magnitudes of the drifts confirm the predictions of the self-similar analysis (§ 4.2) and are consistent with previous simulations of other variable-density mixing layers (e.g. Pantano & Sarkar Reference Pantano and Sarkar2002; Almagro et al. Reference Almagro, García-Villalba and Flores2017). For the range of Atwood numbers simulated, $\eta _{12} < \eta _{1} < \eta _{2} < 0$: the Reynolds stress peak is located further in the light fluid than the neutral point of mean streamwise velocity, which itself is further than the peak of mean cross-stream velocity.

Figure 11. (a) Favre mean streamwise velocity profile edge position (10 %, as defined in (4.14)) evolutions for the range of Atwood numbers (coloured as in figure 3) are shown by solid lines. Their neutral positions (i.e. $\skew2\tilde {U}_1=0$) are also shown by dotted lines. (b) The self-similar positions are compared as a function of Atwood number for density neutral point (${\blacksquare }$), peak cross-stream velocity $\eta _2$ (${\blacktriangle }$), streamwise velocity neutral point $\eta _1$ (${\blacktriangledown }$) and peak Reynolds stress $\eta _{12}$ (${\blacklozenge }$).

5.4. Velocity fluctuation intensity profiles

Statistical profiles for velocity fluctuations are similarly obtained using self-similar scaling applied to the $y$ coordinate. It has also been verified that these profiles collapse over the self-similar growth time period (apart from a small amount of statistical variability) when scaled in this manner. These time-averaged profiles are compared among Atwood numbers in figure 12.

Figure 12. Self-similar profiles for all Atwood numbers (coloured as in figure 3) showing velocity variances and Reynolds stresses. In (a–d), the solid line represents the velocity variance $\langle u_i' u_j' \rangle$ while the dashed line represents Favre-averaged Reynolds stress $\tilde {R}_{ij}=\langle \rho u_i'' u_j'' \rangle /\bar {\rho }$. (e) $\tilde {R}_{12}$ as in (d) but scaled by ${\rm \Delta} U (\textrm {d} h/\textrm {d} t)$ as suggested by the self-similar analysis. (f) Compares the Favre shear stress of (d) not scaled by local average density but $R_{ij}=\langle \rho u_i'' u_j'' \rangle$ scaled by the average of the free-stream densities $\rho _0$.

Overall, the behaviours of the velocity variances for the low Atwood number case agree well with other published single-density mixing layer simulations. However, there can be significant differences in the magnitudes. The peak variance magnitudes of Rogers & Moser (Reference Rogers and Moser1994) are 23 % larger than those of the present $A=0.001$ simulation. The peak magnitudes of the density ratio $1$ simulation of Almagro et al. (Reference Almagro, García-Villalba and Flores2017) are on average 52 % larger than those of the present simulation. The magnitudes for the Reynolds stress $\tilde {R}_{12}$ peak likewise differ between the simulations by similar amounts. The spatially developing mixing layer simulations of Attili & Bisetti (Reference Attili and Bisetti2012) that reach relatively high Reynolds numbers have peak magnitudes on average 19 % greater than the present results.

One factor likely contributing to the differences of intensity magnitude is the determination of self-similar averaging time. With the present initial disturbance, an overshoot in the turbulence intensities occurs, and after a significant period of time the overshoot decays and asymptotes to the final self-similar growth state as the mixing layer thickens. Other simulations approach self-similar growth differently, and the self-similar period may be determined differently. Despite the difference of the spatial vs. temporally developing configuration, the Attili & Bisetti (Reference Attili and Bisetti2012) intensity profiles appear to agree most closely with the present simulation. Their simulation attains higher Reynolds number and greater thickness growth than the other temporal simulations cited. Differences in simulation domain sizes could potentially alter the turbulence dynamics by restricting structure growth and thereby affect fluctuation intensities. An additional factor may be persisting effects of the differing initial disturbances. Among experiments, there is significant scatter in the intensity magnitudes, e.g. the differences between Bell & Mehta (Reference Bell and Mehta1990) and Spencer & Jones (Reference Spencer and Jones1971) as shown in Almagro et al. (Reference Almagro, García-Villalba and Flores2017). Rogers & Moser (Reference Rogers and Moser1994) summarized the wide range of magnitudes for streamwise velocity variances measured in experiments (as well as the mixing layer growth rates, which are closely related to $\langle u_1' u_2' \rangle$). They also noted the perspective of Dimotakis & Brown (Reference Dimotakis and Brown1976) that persisting influence of the initial conditions may be responsible.

When Atwood number is increased, the behaviour of the intensities and Reynolds stresses remain similar to the $A=0.001$ case, except they shift to the light fluid side and generally decay slightly in magnitude. As shear moves to the light fluid side with increasing Atwood number, the turbulence intensity peak moves to the light fluid side as well. (The close relation between mean shear and the production of turbulent kinetic energy $\tilde {k} = \tilde {R}_{ii}/2$ is apparent from the shear production term $- \bar {\rho } \tilde {R}_{12} \skew2\tilde {U}_{1,2}$ that dominates the budget for $\bar {\rho } \tilde {k}$.) Velocity variances $\langle u_i' u_j' \rangle$ and Reynolds stresses $\tilde {R}_{ij}$ (figure 12a–d) both increasingly shift to the light fluid side with increasing Atwood number; this applies to the on-diagonal ($i=j$) elements as well as the streamwise and cross-stream ($i=1$ and $j=2$).

Figure 12(a–c) suggests that there is only a weak reduction in peak turbulent kinetic energy with increasing Atwood number. The reduction in peak $\tilde {R}_{12}$ Reynolds stress (or $\langle u_1' u_2' \rangle$) is as strong as that experienced by any of the on-diagonal turbulent kinetic energy contributions, yet it is reduced by no more than about 30 % from $A=0.001$ to $A=0.87$. When $\tilde {R}_{12}$ Reynolds stress is scaled by ${\rm \Delta} U$ and $\textrm {d} h/\textrm {d} t$ as suggested by the self-similar analysis, rather than by ${\rm \Delta} U^2$ as is typically reported, the peak magnitudes weakly increase with increasing Atwood number (figure 12e).

If Reynolds stress is scaled using the average density of the two free streams ($\rho _0$) rather than the local mean density, the reduction in peak value with Atwood number is enormous (figure 12f). This is further confirmation that the intense turbulent motions move to (and are sustained in) light density fluid; $\langle u_i u_j \rangle$ and $\tilde {R}_{ij}= \langle \rho u_i'' u_j'' \rangle /\bar {\rho }$ agree very closely for even the highest Atwood number throughout self-similar growth (while there are significant differences during transition with high $A$). This agreement is remarkable because these quantities do not agree well with $R_{ij}=\rho \tilde {R}_{ij}$ due to the shift of strong fluctuations to fluid on average lighter than $\rho _0$. In other words, at elevated $A$, $\bar {\rho }$ is much smaller than $\rho _0$ at the position of peak turbulence intensity, but $\langle \rho u_i'' u_j'' \rangle$ is also commensurately smaller so their ratio is nearly the same as for low $A$. Details of the local density distributions and how they correlate with velocity-based fluctuations will be further considered (§ 6).

5.5. Analysis of thickness growth rate during self-similar growth

The statistical profiles discussed above can be related to growth rate attained by each mixing layer during its self-similar growth regime. The average growth rates calculated over the self-similar growth time intervals obtained above are first summarized as a function of Atwood number. At very low Atwood number ($A=0.001$), the momentum thickness growth rate $\textrm {d}\delta _m/\textrm {d} t/{\rm \Delta} U=0.012$ agrees well with the value of 0.014 reported by Rogers & Moser (Reference Rogers and Moser1994). Almagro et al. (Reference Almagro, García-Villalba and Flores2017) reports a somewhat higher growth rate of 0.017 when density is constant. In terms of thickness measured by $h$, the present simulation's growth rate $\textrm {d} h/\textrm {d} t/{\rm \Delta} U$ of 0.069 is consistent with the 0.062 value for the simulation of Rogers & Moser (Reference Rogers and Moser1994), although both of these growth rates are toward the lower end of the $0.06$–$0.11$ values typically observed in experiments (Dimotakis Reference Dimotakis, Murthy and Curran1991; Pope Reference Pope2000).

Assessing the growth rate reductions as a function of Atwood number across the present simulations, figure 13(a) shows that the momentum thickness growth rate for $A=0.87$ is reduced by 77 % from the rate for the single-density case, while the momentum thickness per mass growth rate ($\delta _{m,pm}$) and the analogous integral growth rate for the density profile ($\delta _{\rho }$) experience lesser but nonetheless significant reductions. The stronger reduction for $\delta _m$ can largely be explained by a misalignment that develops between density and velocity profiles (§ 5.3). The reductions in growth rate based on $h$ and $h_{\rho }$ (figure 13b) are similar to those of the $\delta _{m,pm}$ and $\delta _{\rho }$ integral quantities (figure 13a); the reductions for $h_{X_p}$ are also similar (figure 13d). The thicknesses derived from mean profile $y$-gradients ($\delta _\omega$ and $\delta _{\textrm {d}\rho /\textrm {d} y}$), however, display less smooth growth rate reduction behaviour (figure 13c). This could be a consequence of greater sensitivity to noise associated with a lack of statistical convergence when the maximum gradient is calculated on the smoothed profile, but could also appear due to the greater sensitivity of gradient-based measurements to the details of the profile asymmetries that appear with increased Atwood number. Although $h_{X_p}$ is an integral measurement and therefore lacks the extreme sensitivity to local noise in the mean profile of gradient-based measurements, it appears to display less smooth reductions in growth rate than the other integral thickness quantities. This suggests that some measurements may be particularly sensitive to specific features of the profile shapes. The close correspondence between most of the growth rates (particularly for $\delta _m$, $\delta _{m,pm}$ and $h$) confirms that any of the corresponding thickness measurements would be acceptable for scaling the flow statistics profiles during self-similar growth. Generally, the growths of density thickness quantities (due to mixing of fluids) also behave similarly to the growths of velocity thickness quantities.

Figure 13. The effects of Atwood number on growth rate are displayed for a variety of thickness measurements based on the mean velocity profile (blue) and the mean density profile (dashed green). The momentum thickness calculated from the mean velocity profile and weighted by density (orange) is shown in (a). The fast reaction product thickness $h_{X_\rho }$ in (d) is based solely on the density profiles.

To compare the growth rate effects of Atwood numbers for other variable-density mixing layers, figure 14 includes single-species variable-temperature simulations of Almagro et al. (Reference Almagro, García-Villalba and Flores2017) (low-speed limit) as well as Pantano & Sarkar (Reference Pantano and Sarkar2002) (moderately compressible 0.7 convective Mach number). In contrast to the present species mixing governed by simplified INBM (incompressible non-Boussinesq mixing) equations, the latter cases with thermal variations are governed by the LMNOB (low-Mach number non-Oberbeck–Boussinesq) and fully compressible non-Oberbeck-Boussinesq (NOB) equations (Livescu Reference Livescu2020). A detailed comparison between simplified INBM equations and LMNOB equations has been made at $A=0.75$ by Baltzer & Livescu (Reference Baltzer, Livescu, Livescu, Battaglia and Givi2020). Atwood number for the wide range shown in figure 14 affects both momentum and vorticity thickness growth rates relatively similarly between density variation mechanisms (figure14), particularly given the differences between simulations in addition to species vs. thermal transport, e.g. domain sizes, initial disturbance, determination of self-similar growth period, etc. The growth rates are normalized by the $A \approx 0$ growth rate for each configuration (simulation set), and doing so conceals important physical differences and their effects between cases: the growth rate of the $A=0$ Pantano & Sarkar (Reference Pantano and Sarkar2002) mixing layer had reduced by $40\,\%$ relative to their low-speed simulation solely due to compressibility effects. Their simulations investigating the range of Atwood numbers are only available at a Mach number of 0.7. Since vorticity thickness is based on the maximum gradient of the mean streamwise velocity profile, more scatter in these growth rates is to be expected as they are sensitive both to the shape of the profile and noise in this quantity (smoothing was applied when calculating for the present simulations).

Figure 14. Atwood number effects on mixing layer growth rates measured by (a) momentum thickness and (b) vorticity thickness for the (${\blacksquare }$) present incompressible INBM simulations, (${\bullet }$) the LMNOB (thermal variation) simulations of Almagro et al. (Reference Almagro, García-Villalba and Flores2017) and (${\blacktriangle }$) 0.7 Mach number NOB (thermal variation) simulations of Pantano & Sarkar (Reference Pantano and Sarkar2002). The growth rates are scaled by the corresponding growth rate with $A\approx 0$ for each data set.

The influence of Atwood number on growth rate can be further explained based on the behaviour of the statistical profiles. One useful property of the momentum thickness definitions (4.11) and (4.12) is that they straightforwardly lead to relations between growth rates and flow statistical quantities. This was explored by Vreman et al. (Reference Vreman, Sandham and Luo1996), who showed that an informative relation can be formed based on the time derivative of (4.11) that yields

(5.4)\begin{equation} \frac{\textrm{d}\delta_m}{\textrm{d} t} = -\frac{1}{\rho_0 {\rm \Delta} U^2} \int \frac{\textrm{d}}{\textrm{d} t} \left( \bar{\rho} \tilde{U}_1 \tilde{U}_1 \right) \textrm{d} y \end{equation}

as other terms vanish using the $y$-integrated averaged continuity and momentum equations. Multiplying the Favre mean momentum equation (4.2) for $i=1$ by $\skew2\tilde {U}_1$ produces an equation relating the time derivative of mean kinetic energy to Reynolds stress as

(5.5)\begin{equation} \frac{\textrm{d}}{\textrm{d} t} \left( \bar{\rho} \tilde{U}_1 \tilde{U}_1 \right) = \left(\bar{\rho} \tilde{U}_1 \tilde{U}_1 \tilde{U}_2 + 2\bar{\tau}_{12} \tilde{U}_1 - 2\tilde{U}_1\bar{\rho}\tilde{R}_{12}\right)_{,2} + 2\tilde{U}_{1,2}\bar{\rho} \tilde{R}_{12} - 2\bar{\tau}_{12}\tilde{U}_{1,2}. \end{equation}

Many terms are cast in conservative forms, so they vanish when integrated over the $y$ domain in (5.4). Then, (5.4) becomes

(5.6)\begin{equation} \frac{\textrm{d}\delta_m}{\textrm{d} t} \approx -\frac{2}{\rho_0 {\rm \Delta} U^2} \int_{-\infty}^{\infty} \bar{\rho} \tilde{R}_{12} \tilde{U}_{1,2}\,\textrm{d} y \end{equation}

after neglecting the mean viscous term, which is consistent with the self-similar analysis (§ 4.2) and has been shown to be small by scaling arguments of Pantano & Sarkar (Reference Pantano and Sarkar2002). The growth rates calculated from this equation using the self-similar averaged profiles match those directly measured in the flow to within several per cent.

The above relation shows that the momentum thickness growth rate depends on the density, mean streamwise velocity and Reynolds shear stress profiles. As shown in figure12(d), there is relatively little change in Favre-averaged Reynolds shear stress with the density normalized by the local mean. However, when normalized by the average of the two free-stream densities $\rho _0$, the $\bar {\rho } \tilde {R}_{12}$ quantity that appears in (5.6) reduces strongly with increasing Atwood number, as shown in figure 12(f). The simulation suite was designed such that the average of the free-stream densities $\rho _0$ remains consistent across all Atwood numbers. A consequence is that, if the fluids were fully mixed in equal proportion in the core of the mixing layer, the density there would be the same for every Atwood number. Therefore, reductions as observed in the aforementioned quantity normalized by $\rho _0$ are indicative of $\tilde {R}_{12}$ having its strongest magnitude increasingly further in the light fluid. Furthermore, this density becomes smaller relative to $\rho _0$ as Atwood number increases. In (5.6), the growth rate is the product of $\bar {\rho } \tilde {R}_{12}$ with the mean streamwise velocity gradient, which weakly changes in magnitude with Atwood number (figure 9b). This density weighting reflects the dependence on density in the momentum thickness definition (4.11). Thus, the principal cause of growth rate reduction for $\delta _m$ is the turbulent motions and the associated momentum deficit shifting to the light fluid side.

The dominance of the profile shifting effect is demonstrated by artificially realigning the profiles such that the peak in $\tilde {R}_{12}$ and inflection point in $\skew2\tilde {U}_1$ are returned to the point where $\bar {\rho } = \rho _0$ (figure 15). Eliminating the shifts in this way significantly weakens the growth rate reduction effect. The magnitudes of growth rate reductions after these artificial shifts are approximately those observed for the growth rate of momentum thickness per mass (which lacks the density weighting). For other turbulent variable-density mixing layers, Pantano & Sarkar (Reference Pantano and Sarkar2002) and Almagro et al. (Reference Almagro, García-Villalba and Flores2017) have developed semi-empirical formulas that estimate momentum thickness growth rate reductions as functions of density ratio (or Atwood number) largely based on the shifts that develop between the mean streamwise and mean density profiles.

Figure 15. (a) Momentum thickness growth rate predicted from self-similar averaged statistical profiles using (5.6). The total momentum growth rate (—$\blacksquare$—) is decomposed into light fluid (- - $\blacktriangledown$ - -) and heavy fluid ( - - $\blacktriangle$ - -) contributions according to (5.7). The hypothetical growth rate predicted by (5.6) if the profiles’ drifts were artificially removed ($\cdots \blacklozenge \cdots$) highlights that much of the growth rate reduction is associated with the intense turbulence and shear shifting to the light fluid. (b) Momentum thickness per mass growth rate predicted from (5.11) (—$\bullet$—). The mean shear-Reynolds stress term (- -$\bullet$- -) has the same form as the production for TKE and dominates the other growth rate terms. However, variable-density terms also appreciably reduce the growth rate at high Atwood numbers, as shown by the distance between the curves.

Equation (5.6) also shows that the mixing layer growth rate is the integral over the entire width of the mixing layer of a term that is closely related with the production of TKE (through the shear production mechanism, $- \bar {\rho } \tilde {R}_{12} \skew2\tilde {U}_{1,2}$). It is informative to split this integral into contributions from the light fluid and heavy fluid sides of the domain

(5.7)\begin{equation} \frac{\textrm{d}\delta_m}{\textrm{d} t} \approx -\frac{2}{\rho_0 {\rm \Delta} U^2} \int_{-\infty}^{y_{\rho_0}} \bar{\rho} \tilde{R}_{12} \tilde{U}_{1,2}\,\textrm{d} y -\frac{2}{\rho_0 {\rm \Delta} U^2} \int_{y_{\rho_0}}^{\infty} \bar{\rho} \tilde{R}_{12} \tilde{U}_{1,2}\,\textrm{d} y, \end{equation}

where $y_{\rho _0}$ is the $y$ value at which $\bar {\rho }(y) = \rho _0$. Since mean density increases monotonically across the interior of the mixing layer, it follows that the first term is the light fluid contribution to the growth rate and the second term is the heavy fluid contribution (in a mean sense). The $y_{\rho _0}$ position remains at $y=0$ for $A \rightarrow 0$, thus splitting the mixing layer in half. As Atwood number increases, the $y_{\rho _0}$ position moves much more weakly compared to points on the mean velocity profiles or $\tilde {R}_{12}$ (and $y_{\rho _0}$ instead moves toward the heavy fluid side), as shown by figure 11(b).

When the density differences are very weak ($A=0.001$), the contributions are essentially equally split between the light and heavy sides (figure 15a). The heavy fluid growth rate contribution monotonically decays with Atwood number, which is consistent with the intense turbulence and momentum deficit (relative to the free streams) drifting to the light fluid side. As Atwood number increases from $A=0.001$ to $0.25$, the light fluid growth rate contribution weakly increases, but then decays for higher Atwood numbers. These growth rate changes can be interpreted in light of the time histories of mean profile positions for various Atwood numbers (figure 11). Beginning from the symmetric growth of the mean streamwise profile edge positions for $A=0.001$ in figure 11(a), the edge positions reduce in their penetration to the heavy fluid side but increase in penetrating the light fluid side for Atwood numbers up to $A=0.75$. The penetrations into the light fluid sides stagnate as Atwood number increases beyond $0.75$, and by $A=0.87$ the growth rate has decreased so much that even the growth into the light fluid has reduced below that for $A=0.75$ while growth into the heavy fluid side is negligible. This reduction in growth into the light fluid side is manifested in the sharply reducing light fluid growth contribution in figure 15(a).

While much of the reduction in momentum thickness growth rate with increasing Atwood number can be explained by the velocity neutral point moving into the light fluid, weaker although significant reductions in growth rate also occur in all other thickness measurements, despite their lack of explicit density profile dependencies. This differs from the assumption of constant temporal growth rates with Atwood number sometimes used to develop theories addressing growth rate for spatially developing variable-density mixing layers, as described in the introduction. To address these subtler reductions in growth rate with Atwood number, we derive an analogous equation relating the growth rate of momentum thickness per mass to statistical profiles of the flow. A similar derivation instead beginning from (4.12) leads to

(5.8)\begin{align} \frac{\textrm{d}\delta_{m,pm}}{\textrm{d} t} = \frac{U_{-} + U_{+}}{{\rm \Delta} U^2} \int \frac{\textrm{d} \bar{U}_1}{\textrm{d} t} \textrm{d} y - \frac{1}{{\rm \Delta} U^2} \int \frac{\textrm{d}}{\textrm{d} t}\left(\bar{U}_1 \bar{U}_1\right) \textrm{d} y = -\frac{1}{{\rm \Delta} U^2} \int \frac{\textrm{d}}{\textrm{d} t}\left(\bar{U}_1 \bar{U}_1\right) \textrm{d} y, \end{align}

since $U_{-} = -U_{+}$ for this flow configuration. Developing an expression for the time derivative in the integral based on the Reynolds mean momentum equation leads to

(5.9)\begin{align} \left(\bar{U}_1 \bar{U}_{1}\right)_{,t} &= 2 \left[ \left(-\bar{U}_1 \bar{U}_1 \bar{U}_j - \bar{U}_1 \left\langle u_1' u_j'\right\rangle \right)_{,j} + \bar{U}_{1,j} \bar{U}_1 \bar{U}_j + \bar{U}_{1,j} \left\langle u_1' u_j'\right\rangle \vphantom{\left\langle \frac{p_{,1}}{\rho} \right\rangle}\right. \nonumber\\ &\quad \left. -\left\langle \frac{p_{,1}}{\rho} \right\rangle \bar{U}_1 + \left\langle \frac{\tau_{1j,j}}{\rho} \right\rangle \bar{U}_1 + \bar{U}_{j,j} \bar{U}_1 \bar{U}_1 + \left\langle u_{j,j}' u_1'\right\rangle \bar{U}_1 \right]. \end{align}

Again using the homogeneities of this flow and noting that both $\bar {U}_2$ and $\langle u_1' u_j'\rangle$ vanish on the boundaries (so the conservative terms integrate to 0) simplifies the expression to

(5.10)\begin{align} \frac{\textrm{d}\delta_{m,pm}}{\textrm{d} t} &= -\frac{2}{{\rm \Delta} U^2} \int \left[\bar{U}_{1,2} \left\langle u_1' u_2'\right\rangle + \bar{U}_{1,2} \bar{U}_1 \bar{U}_2 + \left\langle u_{j,j}' u_1'\right\rangle \bar{U}_1 \vphantom{\left\langle \frac{p_{,1}}{\rho} \right\rangle}\right. \nonumber\\ &\quad \left. + \bar{U}_{2,2} \bar{U}_1 \bar{U}_1 - \left\langle \frac{p_{,1}}{\rho} \right\rangle \bar{U}_1 + \left\langle \frac{\tau_{1j,j}}{\rho} \right\rangle \bar{U}_1 \right] \textrm{d} y. \end{align}

It can be shown that the terms based on the Reynolds mean cross-stream velocity decay as the flow evolves (and are of small magnitude, figure 10). In contrast, the $\langle u_{j,j}' u_1'\rangle$ and $\langle {p_{,1}}/{\rho } \rangle$ terms, as well as the dominant $\bar {U}_{1,2} \langle u_1' u_2'\rangle$ term, can be shown to maintain constant magnitudes with the appropriate self-similar scaling. Retaining only these terms, the relation simplifies to

(5.11)\begin{equation} \frac{\textrm{d}\delta_{m,pm}}{\textrm{d} t} \approx -\frac{2}{{\rm \Delta} U^2} \int \left[ \bar{U}_{1,2} \left\langle u_1' u_2'\right\rangle + \left\langle u_{j,j}' u_1'\right\rangle \bar{U}_1 - \left\langle \frac{p_{,1}}{\rho} \right\rangle \bar{U}_1 \right] \textrm{d} y. \end{equation}

For this particular flow in which kinematic viscosity maintains a constant value everywhere, the mean viscous term simplifies to $\nu \bar {U}_{1,22} \bar {U}_1$. It decays with time and during the self-similar growth generates a very small effect, so it is neglected. Of the variable-density terms, $\int \langle {p_{,1}}/{\rho } \rangle \bar {U}_1 \,\textrm {d} y$ is negative and dominates in magnitude over $\int -\langle u_{j,j}' u_1'\rangle \bar {U}_1 \,\textrm {d} y$, which is positive, for all of the Atwood numbers considered. Note that in the single-density case, with $\bar {U}_2=0$ and $u_{j,j}$ everywhere zero, this equation simplifies to ${\textrm {d}\delta _{m,pm}}/{\textrm {d} t} = -2/{\rm \Delta} U^2 \int \bar {U}_{1,2} \langle u_1' u_2'\rangle \,\textrm {d} y$, which is consistent with the usual momentum thickness growth equation (5.6).

This growth rate equation indicates that variable-density effects can modify the $\delta _{m,pm}$ growth rate from its single-density value both through the additional terms (related to the density variations and corresponding divergence of the velocity field) and through changes in the mean velocity gradient $\bar {U}_{1,2}$ and/or the Reynolds stress $\langle u_1' u_2'\rangle$ that appear in the mean shear-Reynolds shear stress term. Applying this equation to each Atwood number (figure 15b) demonstrates that the mean shear-Reynolds shear stress term dominates the growth rate equation for all Atwood numbers and significantly decreases in magnitude for the highest Atwood numbers. To understand the reduction in growth rate associated with the dominant $\bar {U}_{1,2} \langle u_1' u_2'\rangle$ term, it can be seen that there is little change in peak $\bar {U}_{1,2}$ value (figure 9) while peak $\langle u_1' u_2'\rangle$ magnitude decreases moderately (figure 12) with increasing Atwood number. A combination of lower peak stress magnitude and subtle changes in the mean streamwise velocity and stress profile shapes account for the reduction in growth rate. These phenomena also contribute to the conventional (density-weighted) momentum thickness growth rate reductions but are weaker than the effect of shifting mean streamwise velocity and density profiles.

5.6. Profiles involving density fluctuations

Density-velocity correlations are contained in the normalized mass flux quantity $a_i=\langle \rho ' u_i' \rangle /\bar {\rho }$. The turbulent mass fluxes also quantify the relationship between Favre and Reynolds averages for the velocity and Reynolds stress quantities considered above. Relations include $\skew2\tilde {U}_i - \bar {U}_i = u_i' - u_i'' = a_i$, $a_i = - \langle u_i'' \rangle$ and $\bar {\rho } \tilde {R}_{ij} = \bar {\rho } \langle u_i' u_j' \rangle - \bar {\rho } a_i a_j + \langle \rho ' u_i' u_j' \rangle$ (with additional identities provided by Livescu et al. Reference Livescu, Ristorcelli, Gore, Dean, Cabot and Cook2009). As observed in figure 10, the Favre-average cross-stream velocity $\skew2\tilde {U}_2$ is much larger in magnitude than the Reynolds-average cross-stream velocity $\bar {U}_2$. According to the above relations, $\skew2\tilde {U}_2 = \bar {U}_2 + a_2$ is dominated by the $a_2$ turbulent mass flux while the Reynolds mean is relatively insignificant and is explained by the relation developed in appendix A. In single-fluid incompressible temporal mixing layers, the mean cross-stream velocity is zero, in order to satisfy the divergence-free condition. Thus, the correlations between the cross-stream velocity and density in these variable-density mixing layers dominate the Favre mean cross-stream velocity. In addition, density fluctuation properties as revealed by fluctuation intensity are relevant to the structure of the flow.

Two types of normalizations for $a_i$ and density fluctuation intensity are considered. The $\bar {\rho }$ denominator included in the $a_i$ definition removes the density dimensional dependency, but a consequence is that $a_i$ generally grows with Atwood number as density fluctuations become more pronounced. Likewise, as Atwood number increases, density fluctuations increase in magnitude. No mix, denoted by nm, corresponds to the quantities’ values in a hypothetical configuration in which the two fluids are distributed without any mixing, and results in the highest possible magnitudes for $\langle \rho '^2 \rangle$. It can be shown that $\langle \rho '^2 \rangle _{nm} = ({\rm \Delta} \rho /2)^2 = \rho _0^2 A^2$ (Livescu & Ristorcelli Reference Livescu and Ristorcelli2008). By analogy, $A {\rm \Delta} U$ is an appropriate scaling for $a_i$ that is equivalent to scaling by $[{\rm \Delta} \rho /(2 \rho _0)] {\rm \Delta} U$. Figure 16(c,d) shows the self-similar profiles of (a,b) with these scalings applied.

Figure 16. Profiles of correlations involving density scaled using $\rho _0$: (a) normalized mass flux $a_i/{\rm \Delta} U$ with $a_1$ as solid lines and $a_2$ as dashed lines and (b) density variance $\langle \rho '^2 \rangle /\rho _0^2$. Profiles of correlations involving density scaled using no-mix (nm) intensities based on Atwood number: (c) normalized mass flux $a_i/(A {\rm \Delta} U)$ with $a_1$ as solid lines and $a_2$ as dashed lines and (d) density variance $\langle \rho '^2 \rangle /\langle \rho '^2 \rangle _{nm}$. Line colours represent the different Atwood number simulations as in figure 3.

In contrast to the velocity fluctuations’ peak magnitudes shift to the light fluid side, the strongest density fluctuations position in the heavy fluid side with increasing Atwood number (figures 16b,d). Fluctuation profiles for these density-based quantities (intensity and $a_i$) are also almost completely symmetric about $y=0$ for the $A=0.001$ case as the intense turbulence remains at the initial interface position. At increased Atwood number, large-scale disturbances (e.g. corrugations) form on the heavy fluid side while the fine scales of motion producing faster mixing are concentrated on the light fluid side. This behaviour can be adduced from the mean density and velocity fluctuation intensity profiles and by density visualizations (figure 17). The smoother yet disturbed heavy fluid side interface suggests the dominance of large-scale structures while small scales concentrate at the light fluid side for high Atwood number. The large displacements of the largely unmixed fluids relative to the background density gradient leads to large density fluctuation magnitudes on the heavy fluid side. In contrast, the greater mixing on the light fluid side produces local densities on average nearer to the mean density $\bar {\rho }$, thus resulting in smaller density fluctuation intensities.

Figure 17. Surfaces of density coloured from light (blue) to heavy (red) during self-similar growth at approximately $t{\rm \Delta} U/h_0 = 730$ for (a) $A=0.001$ and (b–d) $A=0.75$. Viewing the $A=0.75$ case from below (d) makes apparent the much finer scales on the light fluid side relative to the heavy fluid side (c). This is consistent with the strongest turbulent vortices being concentrated near the light fluid side. Thickness $h$ as defined based on the velocity field is also included.

In absolute terms (scaled by $\rho _0$ in figure 16b), the density fluctuation intensity magnitudes increase up to $A=0.75$ but then stagnate for higher Atwood number. Peak intensity positions penetrate less and less deeply into the heavy fluid side for increasing Atwood numbers. These effects are likely a consequence of the less energetic turbulence sustained on the light fluid side reducing in ability to overcome the heavy fluid side's inertia and disturb it. Scaling using the mean density $\rho _0$ of the two streams, which remains constant for all Atwood numbers, reveals the relative magnitudes of the fluctuation intensities. The no-mix scaled density fluctuation intensity profiles (figure 16d) increase in magnitude up to $A=0.75$ but decay for higher Atwood number. No-mix intensity is proportional to the differences in densities between streams ${\rm \Delta} \rho$, so scaling by this quantity would be expected to scale out the effect of increasing density differences for Boussinesq mixing. Thus, magnitude differences under this scaling reveal non-Boussinesq effects in the fluctuation intensities.