Variable-density effects in incompressible non-buoyant shear-driven turbulent mixing layers

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, 2013). 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 (2001), Livescu & Ristorcelli (2007) and Wei & Livescu (2012). 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

where the viscous stress, assumed to be Newtonian, is

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

Equation (4) represents the nonzero 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, 2013). Equation (4) can be derived from the mixture rule $\rho=1/\left(Y_{1}/\rho_{1}+Y_{2}/\rho_{2}\right)$ (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 $\left(\rho Y_{m}\right)_{,t}+\left(\rho Y_{m}u_{j}\right)_{,j}=\mathcal{D}\left(\rho Y_{m,j}\right)_{,j}$ (Sandoval, 1995; Cook & Dimotakis, 2001; Livescu & Ristorcelli, 2007). 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 VD governing equations (Livescu, 2013). The dynamic viscosity of mixed fluid obeys a relation analogous to the density: $\mu=1/\left(Y_{1}/\mu_{1}+Y_{2}/\mu_{2}\right)$, 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.

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^{\prime}=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 $\tilde{U}_{i}=\langle\rho u_{i}\rangle/\langle\rho\rangle$ and the Favre fluctuation is $u_{i}^{\prime\prime}=u_{i}-\tilde{U}_{i}$, in contrast to the Reynolds mean, $\bar{U}_{i}=\langle u_{i}\rangle$, and fluctuation, $u_{i}^{\prime}=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 ($\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.

The governing equations (2–5) 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 (2012); Livescu et al. (2010, 2011) 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 (4) 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

Integrating across the $y$ domain, this expression becomes $u_{2}(y_{\mathrm{max}})-u_{2}(y_{\mathrm{min}})=-\mathcal{D}\left\{[\ln{\rho}]_{,2}(y_{\mathrm{max}})-[\ln{\rho}]_{,2}(y_{\mathrm{min}})\right\}$. Since density remains constant at the free streams existing at the upper and lower walls, it follows that $u_{2}(y_{\mathrm{max}})-u_{2}(y_{\mathrm{min}})=0$. Thus, the variable-density equations are consistent with the boundary conditions $u_{2}(y_{\mathrm{min}})=u_{2}(y_{\mathrm{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 Livescu, 2020; Baltzer & Livescu, 2020) if the heat conduction coefficient is constant. More complicated cases such as heat release with chemical reaction necessitates nonzero normal velocity at the boundaries, e.g., Higuera & Moser (1994). 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 (2007) 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 (1995) 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 $\nabla p/\rho^{(n+1)}=\nabla q+\nabla\times\vec{A}+\langle\vec{L}\rangle$, which results in a constant coefficient equation corresponding to the dilatational (curl-free) component, $\nabla q$, and implicit equations for the curl (divergence free), $\nabla\times\vec{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\vec{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 $\nabla p/\rho$ is related to mass conservation, which is enforced to machine precision due to the direct solvers involved. The curl component, $\nabla\times\vec{A}$, is related to the baroclinic production of vorticity. The iterative procedure is performed until the maximum $x$–$z$ planar average squared change in $\nabla\times\vec{A}$ relative to the previous iteration value reduces to 0.01 times the squared value of $\nabla\times\vec{A}$ averaged within the plane, for each component $\alpha$:

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 artifacts. 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 (1998) 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_{\mathrm{max}}\eta$, where $\eta$ is the average Kolmogorov length scale $(\nu^{3}/\epsilon)^{1/4}$ and $k_{\mathrm{max}}=a(2\pi/L)$ for a spectral representation of $N$ grid points in a domain of length $L$. The leading coefficient of the $k_{\mathrm{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 et al., 1991). In DNS intended to maximize Reynolds number, typical values are $1\leq k_{\mathrm{max}}\eta\leq 2$ (Gotoh & Yeung, 2013), with $k_{\mathrm{max}}\eta\approx 1.5$ typical for adequately-resolved DNS of isotropic turbulence (Petersen & Livescu, 2010; Pope, 2000). 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_{\mathrm{max}}\eta\gtrsim 4$ (Horiuti & Fujisawa, 2008) or the alignment of strain rate and vorticity (Hamlington et al., 2008).

In the present mixing layer at negligible Atwood number, $k_{\mathrm{max}}\eta$ for the Fourier spectral representation of each homogeneous direction reaches a minimum of $\approx 1.7$ at early times at the centerline (where turbulence is most developed) and continuously increases thereafter. Pantano & Sarkar (2002) report final values of $k_{\mathrm{max}}\eta\approx 1.0$, and they rely on spatial filtering that was shown to produce a relatively small amount of nonphysical dissipation to improve stability in their simulations. Resolution can also be be quantified in terms of grid spacing relative to the average Kolmogorov scale. Almagro et al. (2017) reported horizontal grid spacing finer than $1.8\eta$ during the self-similar growth, whereas the corresponding values in Pantano & Sarkar (2002) are $3$–$4\eta$. In the present low $A$ simulation, the horizontal grid spacing ($\Delta x$ and $\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_{\mathrm{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_{\mathrm{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 6th order accurate for both the momentum and pressure equations. The uniform grid spacing is finer (reduced to a factor of 0.8: $\Delta_{y}=0.8\Delta_{x}=0.8\Delta_{z}$) in the inhomogeneous direction, in order to compensate for the lower accuracy relative to the Fourier directions. Modified wavenumber analysis for 6th order compact difference equations indicates errors in differentiating modes become larger at higher wavenumbers (Petersen & Livescu, 2010). 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_{\mathrm{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 (4); 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 (3) 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.

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 6 successive mergers; Vreman et al. (1997) 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, 1994; O’Brien et al., 2014).

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 center) and the most intense turbulence drift to the light fluid side as the flow develops (§5.3). The interface is centered 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 (2002) and 3.9 times those of Almagro et al. (2017).

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

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 et al. (1986); Pantano & Sarkar (2002); Olson et al. (2011); O’Brien et al. (2014); Almagro et al. (2017).

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

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 (2002) and Almagro et al. (2017), though 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}=\left(\rho_{1}+\rho_{2}\right)/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}-\Delta\rho/2$ and $\rho_{2}=\rho_{0}+\Delta\rho/2$, where $\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 behaviors 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 (2002). 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.

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 et al., 2008), 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 centered 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 3D energy spectrum obeying a Gaussian behavior at high wavenumbers with $k^{4}$ behavior at low wavenumbers as $E(k)=\left({k}/{k_{0}}\right)^{4}e^{-2\left({k}/{k_{0}}\right)^{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 (2002) 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 $\left\langle u_{i}^{\prime}u_{i}^{\prime}\right\rangle(y)=Ae^{-\frac{1}{2}\left({y}/{\sigma}\right)^{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. (1986) simulations based on measurements of the intensity profile in a mixing layer experiment and to the thickness used by Pantano & Sarkar (2002). The peak amplitude $A$ is specified for peak intensity $\left\langle u_{i}^{\prime}u_{i}^{\prime}\right\rangle$ of $0.03\Delta U^{2}$ by prescribing a $0.1\Delta U$ RMS 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. (1986) (further described in Riley & Metcalfe, 1979) and Pantano & Sarkar (2002) (further described in Pantano-Rubino, 2000), 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, 2000), but the present method produces an initial velocity field divergence fully consistent with the variable-density incompressible divergence condition (4). 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 $\Delta U$ in all simulations shown).

Momentum thickness Reynolds number, $Re_{m}=\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 $\Delta U$ to initialize the simulation is arbitrary and can be scaled out. In consistent units, $\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, 1995; Cook & Dimotakis, 2001; Livescu & Ristorcelli, 2007), though maintaining $Sc=0.7$ (which is typical for gases) is also common (e.g., Olson et al., 2011). 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., 2010).

Integrating the mean density conservation equation (13) over the $y$ domain indicates that $\int_{y_{1}}^{y_{2}}\bar{\rho}\;dy$ is constant with respect to time (total mass within the domain is conserved). The mean momentum equations (14)-(15), when similarly integrated over the $y$ domain, show that $\int_{y_{1}}^{y_{2}}\bar{\rho}\tilde{U}_{i}\;dy$ 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 (14) and (15) 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.

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 centered 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 (10) and Favre mean streamwise momentum equation (11) are satisfied for self-similar growth when the growth rate $dh/dt$ is constant and the mean variables are non-dimensionalized as

Analyzing 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 [$d\hat{U}_{2}/{d\eta}(\eta_{2})=0$] and $\eta_{12}$ as the point where Favre shear stress has its inflection [$d\hat{R}_{12}/{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 (2002) after developing the self-similar analysis framework while analyzing 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).

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 (11). Momentum thickness is defined as:

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_{-}=-\Delta U/2$ and $U_{+}=\Delta U/2$. An analogous thickness could also be defined on a per-mass basis to depend only upon the mean velocity profile:

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, (21) is commonly given as the definition of the momentum thickness because (20) reduces to this when density is constant, though (20) 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

as the vorticity magnitude reduces to $|d\bar{U}/dy|$ 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 specific percents (e.g., 10%) of the difference $\Delta U$ between its free-steam values $U_{-}$ and $U_{+}$ (which are associated with fluids having densities $\rho_{1}$ and $\rho_{2}$, respectively):

While momentum thickness and vorticity thickness have been the most commonly used thickness measurements in the historic mixing layer literature, Pope (2000) adopts $h_{0.1}$ in treating planar mixing layers, and it has also been recently used by Schwarzkopf et al. (2016), for example. For brevity, $h$ will be used herein to indicate $h_{0.1}$. This choice of velocity percent 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$, though 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

Note that $\delta_{\rho}$ is equivalent to the width measurement introduced by Youngs (1991, 2009) and also used by Livescu et al. (2010); it is typically written as $W=\int_{-\infty}^{\infty}F_{1}F_{2}\;dy$, defined based on the mean volume fractions of each species $F_{1}=\left(\bar{\rho}-\rho_{1}\right)/\left(\rho_{2}-\rho_{1}\right)$ and $F_{2}=\left(\rho_{2}-\bar{\rho}\right)/\left(\rho_{2}-\rho_{1}\right)$. 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 (Youngs, 2013; Livescu et al., 2010). 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 (26) 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., 2009; Zhou & Cabot, 2019) are also relevant. One such quantity, used by Cook & Dimotakis (2001); Livescu & Ristorcelli (2008); Livescu et al. (2010), is $h_{X_{\rho}}=\int_{-\infty}^{\infty}X_{P}(\bar{\rho})\;dy$, where $X_{P}$ represents the amount of product in a hypothetical fast reaction between the two species

$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.

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 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.

It should be noted that $\delta_{m,0}$, used for nondimensionalization, 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_{d\rho/dy}$ 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 4c. 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.

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 (1994), 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:

In earlier simulations emphasizing roll-ups of KH vortex structures and their subsequent mergers, Moser & Rogers (1993) showed that large values of $\mathcal{V}$ are associated with these features. Conversely, when Rogers & Moser (1994) 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 (2013) 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 behavior is compared to that of the present simulation with negligible Atwood number in Figure 5a. The present simulation produces a much weaker peak in $\mathcal{V}$ than the Attili & Bisetti (2013) simulation. Despite the weaker peak, the present simulation follows similar behavior of approaching self-similarity after the peak. This behavior contrasts with the asymptoting from below that appears to occur for the fully-turbulent initial condition of Rogers & Moser (1994). 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 behavior also occurs at increased Atwood numbers.

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

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 $\Delta U$. Since $\mathcal{E}$ has units of velocity cubed, it can be argued on dimensional grounds that $\mathcal{E}$ scales with $\Delta U$ only and therefore is constant with respect to time during self-similar growth (Rogers & Moser, 1994). 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 5b. 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 (2013) that begin from a thin disturbance. The latter applies to the simulation of Rogers & Moser (1994) that begins from a field containing fully-turbulent fluid and slowly approaches the self-similar state from below (in terms of dissipation). Attili & Bisetti (2013) 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, though this is weaker than in Attili & Bisetti (2013) 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 (2013) 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., 2009), unlike the typical budgets written for single-fluid incompressible mixing layers (e.g., Rogers & Moser, 1994). Therefore, the integrated dissipation must be divided by density to have the units of $(\Delta U)^{3}$. One option is to nondimensionalize 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:

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 ${\Delta U}^{3}$ decrease strongly with increasing Atwood number, while $\tilde{\mathcal{E}}$ scaled by ${\Delta U}^{3}$ displays a much weaker dependence.

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, 1994), comparing additional flow statistics profiles produces further useful indications. This was recognized by Vreman et al. (1997), 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 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.

Variable-density mixing layers introduce additional quantities to be considered for self-similarity, most importantly the density fluctuation intensity $\langle\rho^{\prime}\rho^{\prime}\rangle$. The integral values of this planar-mean quantity are shown for each Atwood number in Figure 7d. $\langle\rho^{\prime}\rho^{\prime}\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 behavior, 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^{\prime}\rho^{\prime}\rangle$ for $A=0.75$ (blue curve) is suggestive of reaching self-similar growth at particularly late time, with a leveling 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 behavior 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 7b.

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}^{\prime}u_{j}^{\prime}\rangle$, and density fluctuation intensity $\langle\rho^{\prime 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\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., 1996, 1997; Pantano & Sarkar, 2002). In general, questions remain about the universality of the self-similar state (e.g., Dimotakis & Brown, 1976; Rogers & Moser, 1994; Vreman et al., 1997). 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 behavior have been observed to occur at a Reynolds number threshold (i.e., the mixing transition, Dimotakis, 2000). 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 (2000), 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, 1994). 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.

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}\frac{\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 $u^{\prime}_{\mathrm{rms}}\approx\sqrt{\left(\langle u_{1}^{\prime}u_{1}^{\prime}\rangle+\langle u_{2}^{\prime}u_{2}^{\prime}\rangle+\langle u_{3}^{\prime}u_{3}^{\prime}\rangle\right)/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}\frac{1}{\tilde{\varepsilon}\nu}}$. Using $u^{\prime}_{\mathrm{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). Though similar definitions are also used for other relevant flows (e.g., Sekimoto et al., 2016), mixing layer literature often uses $\sqrt{2k}$ as the velocity scale (rather than $u^{\prime}_{\mathrm{rms}}=\sqrt{2k/3})$ to form $Re_{\lambda}=(2k)\lambda_{g}/\nu=k\sqrt{20/(\varepsilon\nu)}$ (e.g., Pantano & Sarkar, 2002; O’Brien et al., 2014; Almagro et al., 2017). Renormalized to the present convention, the $Re_{\lambda}$ range during self-similar growth for the single-density mixing layer of Rogers & Moser (1994) is $84$–$99$ and for Almagro et al. (2017) 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 (2000) 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.

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.

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 centered 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 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.

The positions of the neutral points (i.e., $\rho=\rho_{0}$ and $\tilde{U}_{1}=0$) and positions of extrema for various statistical quantities (e.g., $\min(\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 percent of the free-stream difference $\Delta U$: $y_{\left[\tilde{U}_{1}=U_{-}+0.1*\Delta U\right]}$ and $y_{\left[\tilde{U}_{1}=U_{+}-0.1*\Delta U\right]}$. These are the points whose separation define $h$ in (23). In Figure 11a, 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 $\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 centered at the initial interface. In addition, Figure 11b indicates that the mean cross-stream velocity $\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, 2002; Almagro et al., 2017). 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.

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.

Overall, the behaviors 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 (1994) 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. (2017) 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 (2012) 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 (2012) 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 (1990) and Spencer & Jones (1971) as shown in Almagro et al. (2017). Rogers & Moser (1994) 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}^{\prime}u_{2}^{\prime}\rangle$). They also noted the perspective of Dimotakis & Brown (1976) that persisting influence of the initial conditions may be responsible.

When Atwood number is increased, the behavior 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}\tilde{U}_{1,2}$ that dominates the budget for $\bar{\rho}\tilde{k}$.) Velocity variances $\langle u_{i}^{\prime}u_{j}^{\prime}\rangle$ and Reynolds stresses $\tilde{R}_{ij}$ [Figure 12(a–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-cross-stream ($i=1$, $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}^{\prime}u_{2}^{\prime}\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 $\Delta U$ and $dh/dt$ as suggested by the self-similar analysis, rather than by $\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. $\left\langle u_{i}u_{j}\right\rangle$ and $\tilde{R}_{ij}=\left\langle\rho u_{i}^{\prime\prime}u_{j}^{\prime\prime}\right\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 $\left\langle\rho u_{i}^{\prime\prime}u_{j}^{\prime\prime}\right\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).

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 $d\delta_{m}/dt/\Delta U=0.012$ agrees well with the value of 0.014 reported by Rogers & Moser (1994). Almagro et al. (2017) 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 $dh/dt/\Delta U$ of 0.069 is consistent with the 0.062 value for the simulation of Rogers & Moser (1994), though both of these growth rates are toward the lower end of the $0.06$–$0.11$ values typically observed in experiments (Pope, 2000; Dimotakis, 1991).

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_{d\rho/dy}$), however, display less smooth growth rate reduction behavior (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. Though $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.

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. (2017) (low-speed limit) as well as Pantano & Sarkar (2002) (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, 2020). A detailed comparison between simplified INBM equations and LMNOB equations has been made at $A=0.75$ by Baltzer & Livescu (2020). Atwood number for the wide range shown in Figure 14 affects both momentum and vorticity thickness growth rates relatively similarly between density variation mechanisms (Figure 14), 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 (2002) 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).

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

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

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

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 (2002). The growth rates calculated from this equation using the self-similar averaged profiles match those directly measured in the flow to within several percent.

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 Figure 12(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 (33) 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 (33), 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 (20). 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.