2 The physical model for compressible miscible Newtonian fluids

We present in this section the governing equations and the mixing model for two miscible Newtonian fluids. The motion takes place in a three-dimensional domain $(L_{x},L_{y},L_{z})$ . The heavy fluid (referenced with the subscript $H$ ) is initially located in the upper side of the domain, for $0\leqslant z\leqslant z_{t}$ , on top of a light fluid (subscript $L$ ) occupying the region $z_{b}\leqslant z\leqslant 0$ , where $z_{b}$ and $z_{t}$ denote the coordinates of the bottom and the top of the domain. The mixing of these two compressible miscible fluids is analysed within the framework of the single fluid approximation. The partial densities, $\unicode[STIX]{x1D70C}_{H}$ and $\unicode[STIX]{x1D70C}_{L}$ are defined such that $\int _{V}\unicode[STIX]{x1D70C}_{H}\,\text{d}x\,\text{d}y\,\text{d}z=m_{H}$ and $\int _{V}\unicode[STIX]{x1D70C}_{L}\,\text{d}x\,\text{d}y\,\text{d}z=m_{L}$ , where $m_{H,L}$ are the masses of the heavy and light fluids contained in the total volume $V$ . The ‘partial pressures–partial densities’ mixing model is written

A fluid concentration, $c$ , is also defined such that the partial densities are

Let us recall the full NSEs for a binary mixture of two Newtonian miscible fluids (Hirschfelder, Curtiss & Bird Reference Hirschfelder, Curtiss and Bird1954, chap. 11), (Cook Reference Cook2009)

where $u_{i}$ , $(i=1,2,3)$ , are the three velocity components. The equation of state of each component is

The momentum equation (2.3) contains the buoyancy term $-g\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FF}_{i3}$ , where $g$ is the acceleration due to gravity. The stress tensor, $\unicode[STIX]{x1D70E}_{ij}=\unicode[STIX]{x1D707}(\unicode[STIX]{x2202}_{j}u_{i}+\unicode[STIX]{x2202}_{i}u_{j}-2/3\unicode[STIX]{x1D6FF}_{ij}\unicode[STIX]{x2202}_{\ell }u_{\ell })$ , where $\unicode[STIX]{x1D707}$ is the dynamic viscosity coefficient, is calculated within the Stokes approximation, and $\unicode[STIX]{x1D634}_{ij}=1/2(\unicode[STIX]{x2202}_{j}u_{i}+\unicode[STIX]{x2202}_{i}u_{j})/2$ is the rate-of-deformation tensor. This is valid for monoatomic gases and consistent with the value $\unicode[STIX]{x1D6FE}=5/3$ . The dissipation function in the energy equation is $\unicode[STIX]{x1D70E}_{ij}\unicode[STIX]{x1D634}_{ij}$ . The heat flux expression is $q_{i}=+h_{\unicode[STIX]{x1D6FC}}J_{\unicode[STIX]{x1D6FC}i}-\unicode[STIX]{x1D705}\unicode[STIX]{x2202}_{i}T$ , ( $\unicode[STIX]{x1D6FC}=H,L$ ) where $h_{\unicode[STIX]{x1D6FC}}$ is the enthalpy of the species $\unicode[STIX]{x1D6FC}$ and $J_{\unicode[STIX]{x1D6FC}i}=-\unicode[STIX]{x1D70C}{\mathcal{D}}\unicode[STIX]{x2202}_{i}c_{\unicode[STIX]{x1D6FC}}$ is the diffusive mass flux. The diffusion coefficient of species is ${\mathcal{D}}$ and the thermal conductivity coefficient is denoted  $\unicode[STIX]{x1D705}$ . The molar weights are ${\mathcal{M}}_{H}$ and ${\mathcal{M}}_{L}$ and the perfect gas constant is ${\mathcal{R}}$ . Finally, let us recall that this system contains the three Kovásznay modes, namely the vorticity, entropic and acoustic modes (Chu & Kovásznay Reference Chu and Kovásznay1958), see (Monin & Yaglom Reference Monin and Yaglom1962, § 1.7). System (2.3) is written in a dimensionless form with the following units. The reference of length is the width box, $L_{y}$ (or $L_{x}$ ). The unit of time is $(L_{y}/g)^{1/2}$ . The unit of mass, denoted $\unicode[STIX]{x1D70C}_{r}$ , is given by the half-sum of densities on each side of the pseudo-interface and the temperature reference is given by a uniform temperature, denoted  $T_{r}$ . As a result, the complete dimensionless NSEs for a binary mixing of Newtonian miscible fluids are written

The specific heats at constant volume and pressure of each fluid and their ratios are denoted $C_{v,p;H,L}$ and $\unicode[STIX]{x1D6FE}_{H,L}$ . The mixing adiabatic index is $\unicode[STIX]{x1D6FE}_{m}(c)=C_{p,m}/C_{v,m}$ , where the expression of the mixture-specific heats read $C_{v,p;m}(c)=cC_{v,p;H}+(1-c)C_{v,p;L}$ . The reference concentration is chosen to be $c_{ref}=(1-A_{t})/2$ , where $A_{t}=(\unicode[STIX]{x1D70C}_{H}-\unicode[STIX]{x1D70C}_{L})/(\unicode[STIX]{x1D70C}_{H}+\unicode[STIX]{x1D70C}_{L})$ is the Atwood number. We also use the reference value $\unicode[STIX]{x1D6FE}_{r}=\unicode[STIX]{x1D6FE}_{m}(c_{ref})$ . The derivative with respect to the concentration is denoted $d_{c}$ and the dimensionless difference of specific heats at constant pressure is $\unicode[STIX]{x1D6E5}_{H,L}^{\star }$  (see appendix A). The expression of the stratification parameter and the Reynolds, Schmidt and Prandtl numbers are

where $2/{\mathcal{M}}_{r}=1/{\mathcal{M}}_{H}+1/{\mathcal{M}}_{L}$ and $C_{v,ref}=C_{v;m}(c_{ref})$ is the reference value. System of PDEs (2.3) has to be equipped with appropriate initial and boundary conditions. In physical situations, e.g. in astrophysical configurations or in a small pellet in inertial confinement fusion, the boundaries of the computational domain used here are actually contact surfaces between the fluids. As a result, acoustic waves are partially reflected or transmitted through these contact surfaces. A limit case is studied here where acoustic waves are reflected and remain inside the domain. The opposite limit case where acoustic waves are fully transmitted is also possible (Reckinger et al. Reference Reckinger, Livescu and Vasilyev2016). In any case, the number of boundary conditions that has to be specified is given by the analysis of Strikwerda (Reference Strikwerda1977). The numerical implementation may be achieved by using the approach proposed by Thompson (Reference Thompson1990) and it has been applied in a spectral framework by Boudesocque-Dubois et al. (Reference Boudesocque-Dubois, Clarisse and Gauthier2003), among others. As a result, we chose periodic boundary conditions in the horizontal $(x,y)$ -directions for all physical quantities with stress-free conditions along the top and bottom boundaries for the velocity. There is no flow and therefore no mass flux through the top and bottom boundaries. The temperature is kept fixed at the bottom boundary and there is no heat flux at the top. These boundary conditions are written

There is no boundary condition on the density and pressure.

3 Numerical simulations

4 Numerical results

The RTI phenomenology where both fluid layers are stably stratified is the following. This scenario will be deepened and clarified in the next sections. The linear regime is classically defined by $a(t)k\leqslant 1$ , where $a(t)=a_{o}\cosh (\unicode[STIX]{x1D70E}t)$ . For the unstable scale $k\approx 182.2$ , $a_{o}=10^{-3}$ and $\unicode[STIX]{x1D70E}=4.3$ , the linear regime ends very early at $t\approx 0.43$ . The nonlinear regime of the instability is developing and transition to turbulence occurs, then a turbulent mixing layer begins to develop. This regime is qualitatively close to a classical RT regime. However, quantitatively, this regime is affected by the stratification, since the mixing layer starts to smooth the density jump so that physical quantities start decaying, but at different times and at different rates. The root-mean-square (r.m.s.) density $\overline{\unicode[STIX]{x1D70C}^{\prime 2}}^{1/2}$ maximum occurs at $t\approx 3.25$ and the baroclinic production reaches its maximum at $t\approx 3.15$ , while the global $A_{LS}$ Atwood number decays monotonically from the beginning and vanishes at $t\approx 3.50$ . The first vorticity maximum occurs at $t\approx 4.35$ and the second at $t\approx 6.85$ . The vertical Taylor–Reynolds number reaches its maximum ( $Re_{Tz}\approx 76$ ) at $t\approx 8.30$ and the r.m.s. temperature maximum occurs at $t\approx 8.40$ . The mixing layer thickness saturates around $t\approx 9$ and from this instant the concentration flattens. The simulation is stopped at $t\approx 18.55$ where the vertical Taylor–Reynolds number is only 20. The density-gradient length scale reaches very high values at the beginning of the process. It then decreases and reaches a quasi-uniform value within the mixing layer during the freely decaying turbulent regime (hereafter denoted FD regime).

Figure 2. (a) Zoom of the mean density profiles, $\overline{\unicode[STIX]{x1D70C}}(z,t)$ , at six selected times. (b) The large-scale Atwood number, $A_{LS}(t)$ , as a function of time. The initial value $A_{LS}(t=0)$ , is slightly below the value of the simulation $A_{t}=0.25$ , due to the regularization of the initial density jump (see equation (3.1)) with the small-scale Atwood number, $A_{SS}(t)$ , defined by the root-mean-square density.

5 Anisotropy

Since the flow and mixing anisotropies are a current concern in RT-induced turbulent mixing layers, this point is considered in this section. We use the Reynolds stress anisotropy tensor for the velocity, where no distinction between the various scales is made. We also use a spectral anisotropy indicator for the velocity and the scalars, to have access to the anisotropy of a given scale. The horizontal, $\unicode[STIX]{x1D609}_{ij}(z,t)$ and volume average, $\unicode[STIX]{x1D623}_{ij}(t)$ , of the Reynolds stress anisotropy tensor are defined as

Vertical profiles of the $\unicode[STIX]{x1D609}_{zz}$ -component are displayed in figure 16(a) at six different times. Notice that expressions (5.1) are defined for all values of the $z$ -coordinate, although outside of the turbulent mixing layer the three contributions of the turbulent kinetic energy are very small. At $t=3.25$ , turbulence anisotropy reaches 0.55, which means that $88\,\%$ of the turbulent kinetic energy is in the third direction, while at $t=7.45$ , $64\,\%$ of the turbulent kinetic energy is in the third direction. The time evolution of the three components of the mean value $\unicode[STIX]{x1D623}_{ij}$ are represented in figure 16(b). It also shows a strong anisotropy during the RT regime, where the kinetic energy corresponding to the vertical direction is about $90\,\%$ of the total turbulent kinetic energy. After this maximum, the anisotropy abruptly decreases to $\unicode[STIX]{x1D623}_{zz}\approx 0.40$ between $t\approx 3$ and 4. Afterwards, the isotropy is slowly re-established meanwhile turbulence is decaying. The horizontal diagonal Reynolds stress anisotropy tensor components, $\unicode[STIX]{x1D623}_{xx}$ and $\unicode[STIX]{x1D623}_{yy}$ , are very closed to each other (not included), all along the process. The off-diagonal term, $\unicode[STIX]{x1D623}_{yz}$ (or $\unicode[STIX]{x1D623}_{xz}$ ), is negligible at all times, which shows that turbulence production is essentially, or totally, of baroclinic type, no shear production occurs here. On the other hand, it has already been shown that intermediate scales are isotropic while small scales are not, at low and large Atwood numbers, within the Sandoval model (Livescu & Ristorcelli Reference Livescu and Ristorcelli2008; Chung & Pullin Reference Chung and Pullin2010; Livescu et al. Reference Livescu, Ristorcelli, Petersen and Gore2010). Similar conclusions have also been drawn from a RT Boussinesq large-scale simulation (Schneider & Gauthier Reference Schneider and Gauthier2016a ). This anisotropy is stronger in the middle of the mixing layer where turbulence is stronger. In addition, in this compressible case, acoustic waves may also contribute to this anisotropy.

Figure 17. Velocity spectral anisotropy, as defined by (5.2), at four different times, $t=4.35$ (first vorticity maximum), $t=6.85$ (second vorticity maximum), $t=7.45$ (kinetic energy maximum) and $t=11$ (within the FD regime).

Spectral anisotropy has been investigated by several authors (Livescu & Ristorcelli Reference Livescu and Ristorcelli2008; Chung & Pullin Reference Chung and Pullin2010; Livescu et al. Reference Livescu, Ristorcelli, Petersen and Gore2010; Cabot & Zhou Reference Cabot and Zhou2013). The following indicator for the velocity anisotropy (Chung & Pullin Reference Chung and Pullin2010; Schneider & Gauthier Reference Schneider and Gauthier2016c ) is used

where $\widehat{u}$ denotes the Fourier transform of the vector variable $u$ . Results for the velocity anisotropy, $S_{u}(k_{x},t)$ , are shown in figure 17, at four different times, $t=4.35$ (first vorticity maximum), $t=6.85$ (second vorticity maximum), $t=7.45$ (kinetic energy maximum) and $t=11$ (within the FD regime). At any time, the intermediate scales are isotropic, as opposed to the RT Boussinesq turbulence (Schneider & Gauthier Reference Schneider and Gauthier2016c ). However, from time $t=7.45$ , intermediate scales are more isotropic than the small scales. Within the FD regime, all scales are slightly anisotropic. The scalar anisotropy is evaluated with the following indicator

similar to equation (5.2). The anisotropy variations with respect to the horizontal wavenumber $k_{x}$ are shown in figures 18 and 19 for the concentration and temperature. Isotropy clearly occurs at intermediate scales only, for the four times displayed and for the two scalars, concentration and temperature. We notice a very strong resemblance between the concentration and temperature spectra at a given time. Very similar results have been obtained for the Boussinesq simulation, $Sr0\text{-}Re3\times 10^{4}$ , and the compressible simulation, $Sr6\text{-}Re3\times 10^{4}$ .

Figure 18. Concentration-gradient spectral anisotropy, as defined by (5.3), at four different times, $t=4.35$ , $6.85$ , $7.45$ and $11$ .

Figure 19. Temperature-gradient spectral anisotropy, as defined by (5.3), at four different times, $t=4.35$ , $6.85$ , $7.45$ and $11$ .

6 Averaged turbulent equations

In this section, the numerical data are further analysed by using the Favre-averaged equations. Such an analysis is also relevant for second-order modelling (Grégoire, Souffland & Gauthier Reference Grégoire, Souffland and Gauthier2005; Livescu et al. Reference Livescu, Ristorcelli, Petersen and Gore2010). The equation for the turbulent kinetic energy is first recalled. The right-hand side of this equation contains the three families of terms, i.e. production, diffusion and dissipation. From equations (2.5) and definition (2.10), one obtains the exact averaged equation for the turbulent kinetic energy (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002, § 5.7)

where the expression of the dissipation rate is

The first two terms of the right-hand side of equation (6.1) are producing turbulence, one from the velocity shear, the second from the baroclinic term. The third contribution is the diffusion term, then the pressure velocity–divergence correlation, which is mainly a dissipation term, and the viscous dissipation. The evolution of the turbulent kinetic energy $\langle k\rangle _{\unicode[STIX]{x1D6FD}}$ of the three simulations, $Sr6\text{-}Re6\times 10^{4}$ , $Sr6\text{-}Re3\times 10^{4}$ and $Sr0\text{-}Re3\times 10^{4}$ , are given in figure 20. The kinetic energy of the $Sr6\text{-}Re6\times 10^{4}$ -simulation grows during the RT regime and reaches a first maximum at $t\approx 3.9$ . A second peak is reached at $t\approx 6.55$ , essentially due to acoustic waves. Similar evolution is obtained from the compressible run, $Sr6\text{-}Re3\times 10^{4}$ . The same amplitude is observed in these two simulations with a time shift of $\unicode[STIX]{x0394}t\approx 1.40$ . These evolutions are different to the one obtained from the Boussinesq approximation where the turbulent kinetic energy grows as $t^{2}$ in the turbulent regime. Dissipation rates of the two compressible runs have the same behaviour, i.e. a strong increase in the RT regime followed by two maxima and the FD regime. They also differ from a time shift of $\unicode[STIX]{x0394}t\approx 1.40$ with a factor 1.08 in amplitude. The dissipation rate behaviour obtained from the Boussinesq approximation grows as $t$ at a much larger level. The turnover time scale is estimated as the ratio $\langle \widetilde{k}/\widetilde{\unicode[STIX]{x1D700}}\rangle _{\unicode[STIX]{x1D6FD}}$ and its initial value is often used to measure a flow duration. If we assume that the flow is turbulent when the Taylor–Reynolds number reaches $Re_{Tz}\approx 30$ ( $t=2.55$ ), the initial turnover time scale is $\unicode[STIX]{x1D70F}_{turn-over}=0.43$ . The RT regime thus undergoes several turnover time scales. At later times, the turnover time scale grows continuously and is equal to $\unicode[STIX]{x1D70F}_{turnover}=2.69$ at time $t=3.25$ , $\unicode[STIX]{x1D70F}_{turnover}=4.45$ ( $t=4.15$ ) and $\unicode[STIX]{x1D70F}_{turnover}=13.0$ ( $t=7.45$ ). It is still growing during the FD regime, up to $t\approx 15$ , where it slowly decays. This remark about the turnover time scales also applies to the $Sr6\text{-}Re3\times 10^{4}$ simulation.

Figure 20. (a) Turbulent kinetic energy $\overline{\unicode[STIX]{x1D70C}}\widetilde{k}=\overline{\unicode[STIX]{x1D70C}u_{i}^{\prime \prime }u_{i}^{\prime \prime }}/2$ . Comparison between results obtained from simulations $Sr6\text{-}Re6\times 10^{4}$ (blue), $Sr6\text{-}Re3\times 10^{4}$ (green) and $Sr0\text{-}Re3\times 10^{4}$ (red). (b) Dissipation rate. Comparison between results obtained from simulations $Sr6\text{-}Re6\times 10^{4}$ (blue), $Sr6\text{-}Re3\times 10^{4}$ (green) and $Sr0\text{-}Re3\times 10^{4}$ (red).

Figure 21. Mean profiles related to (6.1). (a) Turbulent kinetic energy $\overline{\unicode[STIX]{x1D70C}}\widetilde{k}$ , dissipation rate $\overline{\unicode[STIX]{x1D70C}}\widetilde{\unicode[STIX]{x1D700}}$ , pressure source term $-1/Sr\overline{u_{z}^{\prime }}\unicode[STIX]{x2202}_{z}\overline{p}$ and pressure–velocity divergence correlation $1/Sr\overline{p^{\prime }\unicode[STIX]{x2202}_{i}u_{i}^{\prime }}$ . The mean effective Atwood number $3\times 10^{-3}A_{t}(t)/A_{t}(0)$ versus time is represented by the black dashed curve. (b) Disequilibrium as defined by (6.3) for two values of $\unicode[STIX]{x1D6FD}=0.8$ ; (blue) and $\unicode[STIX]{x1D6FD}=0.2$ (green). The black dashed line is located at the ordinate 1, the equilibrium value.

Mean value evolutions of the source terms of equation (6.1), i.e. turbulent kinetic energy, dissipation rate and production term $\overline{u_{z}^{\prime }}\unicode[STIX]{x2202}_{z}\overline{p}$ are given in figure 21(a). The $A_{LS}$ Atwood number computed from the $(x,y)$ -averaged density profile is also plotted. The turbulent kinetic energy exhibits two maxima, the first one is associated with the RT instability, the second with acoustic production. During the RT regime, kinetic energy follows closely the production term, $\overline{u_{z}^{\prime }}\unicode[STIX]{x2202}_{z}\overline{p}$ , while dissipation lags behind. However, from the end of the RT regime, where the large-scale $A_{LS}$ Atwood number vanishes, this production term and the dissipation decay while the kinetic energy is still fed by acoustic waves. Turbulent kinetic energy only decays after its second maximum. Turbulence production from velocity shear is negligible, its maximum is $10^{-4}$ , while baroclinic production maximum is of the order of $3\times 10^{-3}$ . The pressure velocity–divergence correlation, which would act as a dissipation term during the RT regime, is also negligible. In the FD regime, it produces turbulence at a very small rate, the maximum being $2\times 10^{-4}$ obtained at $t\approx 8$ . Classically, the turbulent kinetic energy decay may be fitted with a power law of the form $t^{-n}$  (Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014). However, it was not possible to obtain a fit on a reasonable time interval, probably due the perturbation brought by the acoustic waves.

The mixing layer disequilibrium is defined as the ratio of the production to the dissipation (Chung & Pullin Reference Chung and Pullin2010)

A turbulent mixing layer is said to be at the equilibrium when the turbulence is produced and dissipated at the same rate, i.e. when ${\mathcal{D}}_{equil.}=1$ . This quantity is plotted in figure 21(b) for two values of the coefficient $\unicode[STIX]{x1D6FD}=0.8$ and 0.2. It turns out that a RT compressible mixing layer is never in equilibrium. Instead strong departures from equilibrium are observed during the RT regime. At the end of this regime the smallest value is ${\mathcal{D}}_{equil.}=1.22$ at $t\approx 4.35$ , which is also the vorticity maximum. At large times, in the FD regime, the source term changes sign and so does the equilibrium parameter.

Figure 22. Vertical mean profiles of the r.m.s.-density source terms (equation (6.4)). (a) Density-gradient source term $-2\overline{\unicode[STIX]{x1D70C}^{\prime }u_{z}^{\prime }}\unicode[STIX]{x2202}_{z}\bar{\unicode[STIX]{x1D70C}}$ . (b) Velocity-gradient source term $-2\overline{\unicode[STIX]{x1D70C}^{\prime 2}}\unicode[STIX]{x2202}_{z}\tilde{u} _{z}$ .

The turbulent kinetic energy equation (6.1) uses the mass flux $\overline{\unicode[STIX]{x1D70C}^{\prime }u_{z}^{\prime \prime }}=-\overline{\unicode[STIX]{x1D70C}}\overline{u_{z}^{\prime \prime }}$ in the baroclinic production. It is thus of interest to derive the evolution equation of this quantity and to analyse the production terms. It turns out that this equation contains the r.m.s. density. We thus also derive the evolution equation for this quantity. The equation for the r.m.s. density is

The two source terms of the right-hand side are displayed in figure 22 at six selected times. The term proportional to the density gradient follows the RT regime and vanishes at time $t\approx 5$ . At time $t=3.25$ , its contribution is small and for larger times, this term becomes a sink for the r.m.s. density. This density-gradient source term is one order of magnitude larger than the velocity-gradient source term. It is a source in the light fluid but a sink in the heavy fluid part. The equation for the mass flux is written (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002, § 5.5)

Equations (6.4) and (6.5) for the density fluctuations and the mass flux are coupled through their source terms. Density fluctuations and mass flux are mainly created by the density gradient and the pressure gradient, respectively. In a second step, the mass flux and the pressure gradient produce turbulent kinetic energy. The evolution of the three source term mean values of equation (6.5) are represented in figure 23. The term proportional to the pressure gradient is dominant and always positive. Its maximum is reached at $t\approx 3.25$ . Production due to the velocity gradient is always negligible. Production due to the density gradient is small, positive in the RT regime and negative in the FD regime. Its maximum is reached very early at $t\approx 2.80$ .

Figure 23. Mass flux source terms of equation (6.5). (b) Mean values of the three production terms proportional to the density, vertical velocity and pressure gradients, respectively. The pressure-gradient term is dominant and reaches its maximum at $t=3.25$ . (b) Vertical profiles of the mass flux source term proportional to the pressure at six different times.

The energy equation is also written as an equation for the Favre-averaged temperature, $\widetilde{T}$ . It is written

The six terms of the right-hand side have been evaluated and the maximum values observed in the interval $3.25\leqslant t\leqslant 18.55$ are the following. The $p\unicode[STIX]{x2202}_{j}u_{j}$ term is approximately $0.10\times 10^{-1}$ . The dissipation function is smaller, i.e. $0.07\times 10^{-1}$ , the term due to variations of the specific heat is $0.20\times 10^{-1}$ , the enthalpic part of the heat flux contribution is $1.50\times 10^{-1}$ , the diffusive contribution is $0.10\times 10^{-1}$ and the turbulent energy flux, $\unicode[STIX]{x2202}_{z}\overline{\unicode[STIX]{x1D70C}u_{z}^{\prime \prime }T^{\prime \prime }}$ is $0.50\times 10^{-1}$ . The enthalpic part of the heat flux contribution and the turbulent energy flux are therefore the two main contributions. The enthalpic contribution has very large values at the initial times and very quickly decays, while the turbulent thermal flux is zero at the initial time and grows with temperature and velocity fluctuations. To be more precise, at times $t=3.25$ , 4.15 and 8.40, these two contributions have the same magnitude, $(0.015,0.040)$ , $(0.080,0.055)$ and $(0.035,0.035)$ , respectively. These relative values can be explained by the modest Reynolds number reached in this simulation. On the other hand, at short times the full right-hand side is negative in the light fluid side and positive in the heavy fluid. At larger times, this is the opposite, the right-hand side is positive in the light fluid and negative in the heavy fluid. Finally the difference between Favre and Reynolds averaging has been evaluated on the temperature. The difference $\overline{T^{\prime \prime }}=\overline{T}-\widetilde{T}$ is found to be approximately $-2\times 10^{-3}$ in the middle of the mixing layer and $2\times 10^{-3}$ at the edges. Compressibility effects on the temperature field are therefore quite small.

Figure 24. Vertical profiles of r.m.s. thermodynamic fluctuations, density $\unicode[STIX]{x1D70C}_{rms}/\overline{\unicode[STIX]{x1D70C}}$ (blue), pressure $p_{rms}/\overline{p}$ (green), temperature $T_{rms}/\overline{T}$ (red) and concentration $c_{rms}/2$ (yellow), at times $t=3.25$ ; 4.15; 6.85; 7.45; 8.15; and 18.55.

7 Kovàsznay modes: vorticity, entropy and acoustics

The r.m.s. density, pressure, temperature and concentration are first detailed. Vertical profiles of the relative r.m.s.-thermodynamic quantities are displayed in figure 24 at six different times. For the concentration fluctuations the absolute value has been used, and divided by 2 for the sake of clarity. Density and concentration fluctuations take the most important values in the RT regime and quickly decrease in the FD regime. At the contrary, temperature fluctuations are larger in the FD regime. In the RT regime concentration fluctuations follow closely the density fluctuations, while in the FD regime, temperature fluctuations follow closely the concentration fluctuations. Pressure fluctuations remain small along the whole phenomena. Temperature fluctuations are higher in the heavy fluid than in the light fluid. The asymmetry between r.m.s. temperature in the light and heavy fluid side has already been seen in the mean temperature profile (figure 3). Recall that the heavy and light fluids are slightly cooled and heated, respectively (see figure 3 b). The maximum is reached much later than the density fluctuations. Averaged value evolution of density, pressure, temperature and concentration fluctuations, versus time, is given in figure 25. They confirm the previous conclusions drawn from figure 24. Density fluctuations reach their maximum $\langle \unicode[STIX]{x1D70C}_{rms}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}|_{max}=8.6\,\%$ at $t\approx 3.25$ , while temperature fluctuations reach their maximum $\langle T_{rms}/\overline{T}\rangle _{\unicode[STIX]{x1D6FD}}|_{max}=7.8\,\%$ much later, at $t\approx 8.85$ . Pressure fluctuations remain at a low level and reach their maximum $\langle p_{rms}/\overline{p}\rangle _{\unicode[STIX]{x1D6FD}}|_{max}=0.46\,\%$ at $t\approx 4.05$ . Unsurprisingly, concentration fluctuations reach their maximum $\langle \unicode[STIX]{x1D70C}_{rms}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}|_{max}=21\,\%$ at approximately the same time as the density fluctuation maximum, i.e. $t\approx 3.45$ .

Figure 25. (a) Evolution of r.m.s. density, $\langle \unicode[STIX]{x1D70C}_{rms}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}$ , pressure, $\langle p_{rms}/\overline{p}\rangle _{\unicode[STIX]{x1D6FD}}$ , temperature, $\langle T_{rms}/\overline{T}\rangle _{\unicode[STIX]{x1D6FD}}$ and concentration fluctuations, $\langle c_{rms}\rangle _{\unicode[STIX]{x1D6FD}}/2$ versus time. (b) Entropic and acoustic modes of the density and temperature versus time, i.e. $\langle \overline{\unicode[STIX]{x1D70C}_{en}^{\prime 2}}^{1/2}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}$ (blue), $\langle \overline{\unicode[STIX]{x1D70C}_{ac}^{\prime 2}}^{1/2}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}$ (green), $\langle \overline{T_{en}^{\prime 2}}^{1/2}/\overline{T}\rangle _{\unicode[STIX]{x1D6FD}}$ (red), $\langle \overline{T_{ac}^{\prime 2}}^{1/2}/\overline{T}\rangle _{\unicode[STIX]{x1D6FD}}$ (yellow), respectively.

Dynamic compressibility effects may be well understood through the Kovàsznay-mode decomposition. Indeed it is well known that a small-amplitude motion in a compressible fluid may be decomposed into three modes, i.e. vorticity, acoustic and entropic modes. The vorticity behaviour has been detailed in § 4.5. Chassaing et al. (Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002, § 9.6.2) recall that this decomposition is not unique and we follow here their guidelines, in which the acoustic mode is given by the pressure fluctuations, as

and the entropic mode is thus given by

This decomposition for the pressure and density fluctuations has also been used for RT flows by Mellado et al. (Reference Mellado, Sarkar and Zhou2005). The r.m.s. densities of the acoustic and entropic mode are computed from (7.1) and written

and the r.m.s. temperatures of the acoustic and entropic parts are

Figure 25 shows the evolution of the entropic, $\langle \overline{\unicode[STIX]{x1D70C}_{en}^{\prime 2}}^{1/2}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}$ and $\langle \overline{T_{en}^{\prime 2}}^{1/2}/\overline{T}\rangle _{\unicode[STIX]{x1D6FD}}$ , and acoustic part, $\langle \overline{\unicode[STIX]{x1D70C}_{ac}^{\prime 2}}^{1/2}/\overline{\unicode[STIX]{x1D70C}}\rangle _{\unicode[STIX]{x1D6FD}}$ and $\langle \overline{T_{ac}^{\prime 2}}^{1/2}/\overline{T}\rangle _{\unicode[STIX]{x1D6FD}}$ of the density and temperature modes, respectively. These averages have been calculated with $\unicode[STIX]{x1D6FD}=0.8$ . As such, it gives a complementary point of view of figure 26 where mean spatial profiles are given. It shows that the entropic parts are one or two orders of magnitude larger than the acoustic parts. As a result, they are almost equal to their r.m.s. values. The density entropic mode grows in the RT regime and starts decaying right after, while the temperature entropic mode grows much more tardily, as we have already seen above. The density and temperature acoustic modes evolve closely and they maintain their magnitudes during the acoustic regime in the beginning of the FD regime. Figure 26 displays the density, $\overline{\unicode[STIX]{x1D70C}_{ac}^{\prime 2}}$ , and temperature, $\overline{T_{ac}^{\prime 2}}$ , acoustic modes computed from (7.3) and (7.4). We now study the correlation between thermodynamic variables, density, pressure and temperature.

Figure 26. Acoustic-mode vertical profiles of the density $\overline{\unicode[STIX]{x1D70C}_{ac}^{\prime 2}}^{1/2}/\overline{\unicode[STIX]{x1D70C}}$ (a) and temperature $\overline{T_{ac}^{\prime 2}}^{1/2}/\overline{T}$ (b).

Figure 27. (a) Density–pressure correlation, $R_{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x1D70C},p)$ , defined by (7.5) as a function of the vertical coordinate $z$ , at five different times. (b) Density–temperature correlation, $S(\unicode[STIX]{x1D70C},T)$ , defined by (7.7) as a function of the vertical coordinate $z$ , at three different times.

The density–pressure correlation coefficient $R_{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x1D70C},p)$ is defined as (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002, § 9.6.2)

When density fluctuations are essentially of acoustic type, i.e. $\unicode[STIX]{x1D70C}^{\prime }\propto p^{\prime }$ , this coefficient is close to unity. However, since some thermodynamic fluctuations vanish outside the turbulent mixing layer, correlation (7.5), and also correlation coefficients defined by (7.6) and (7.7), are plotted only at collocation points where they are unambiguously defined. This explains the shape of figure 27. Profiles of the density–pressure correlation $R_{\unicode[STIX]{x1D70C}}$ are displayed in figure 27, at five different times. At short times, up to $t\approx 2.7$ , in the early transitional regime, the correlation is positive in the light fluid side ( $z<0$ ) and negative in the heavy fluid ( $z>0$ ). From $t\approx 2.7$ this correlation begins changing sign, i.e. it is negative in the light fluid and positive in the heavy fluid. At larger times, from $t\approx 6$ , the behaviour is non-monotonic in space. Outside the mixing layer boundaries, the correlation $R_{\unicode[STIX]{x1D70C}}$ is not defined, although one may infer that the weak density fluctuations are associated with the acoustic mode, i.e. pressure fluctuations.

A temperature–pressure correlation coefficient may also be defined as

The behaviour of this temperature–pressure correlation is very close to the density–pressure correlation and calls for the same comments.

Finally a density–temperature correlation coefficient is defined as

Figure 27 shows that density and temperature fluctuations are strongly correlated inside the turbulent mixing layer. However, on the boundaries of the mixing layer, at large times when temperature fluctuations are well developed, density and temperature fluctuations are anti-correlated.

8 Statistical study

A statistical study is now conducted. We have successively computed skewness, which measures the PDF asymmetry, and flatness, which gives indications about the Gaussianity of a random variable. Probability density functions are also calculated. They are built through the kernel method with a Gaussian kernel. The interpolation from the spectral grid to a uniform grid is achieved by computing the Lagrange interpolation polynomial locally (Schneider & Gauthier Reference Schneider and Gauthier2016c , appendix B).

Figure 28. Velocity statistics. Skewness, $\overline{u_{z}^{\prime 3}}/\overline{u_{z}^{\prime 2}}^{3/2}$ (a) and flatness, $\overline{u_{z}^{\prime 4}}/\overline{u_{z}^{\prime 2}}^{2}$ (b), of the vertical velocity. The flatness of the Gaussian process is represented by the horizontal dashed line.

Figure 29. Velocity PDFs. (a) Horizontal velocity, $u_{x}^{\prime }$ and vertical velocity, $u_{z}^{\prime }$ . (b) Vertical gradient of the horizontal velocity and vertical velocity components. On each plot, the dashed lines show Gaussian distribution with the same mean and variance.

Figure 30. (a) Pressure PDFs for two different intervals defined by the mean concentration $\overline{c}$ . (b) PDF of the horizontal and vertical potential vorticity components.

The skewness of one horizontal velocity component, $\overline{u_{x}^{\prime 3}}/\overline{u_{x}^{\prime 2}}^{3/2}$ , which refers to large scales, has no particular structure (not included). Skewnesses of the horizontal and vertical derivative of the horizontal velocity, $\overline{(\unicode[STIX]{x2202}_{x,z}u_{x}^{\prime })^{3}}/\overline{(\unicode[STIX]{x2202}_{x,z}u_{x}^{\prime 2})}^{3/2}$ , which refers to intermediate scales, show two peaks of same sign at the layer boundaries. It means that this quantity is weakly skewed in the turbulent region. The skewness and flatness of the vertical velocity component, $\overline{u_{z}^{\prime 3}}/\overline{u_{z}^{\prime 2}}^{3/2}$ and $\overline{u_{z}^{\prime 4}}/\overline{u_{z}^{\prime 2}}^{2}$ , are represented in figure 28. Two peaks of different sign are observed for the skewness, matched by a zero crossing. This zero crossing is slightly shifted to the light fluid. These peaks have the sign of the dominant velocities, i.e. positive on the upper side of the mixing layer and negative on the lower side. Two peaks are also observed for the flatness and departure from Gaussianity is clearly seen. The flatness of the vertical derivative of the vertical velocity $\unicode[STIX]{x2202}_{z}u_{z}$ exhibit a larger departure from Gaussianity (not included). The horizontal velocity flatness (not included) is similar to the vertical velocity flatness but the former is clearly closer to a Gaussian process than the vertical velocity. We may thus anticipate the velocity component PDFs. Figure 29(a) displays the horizontal and vertical velocity component PDFs. Horizontal velocity PDF is very close to the Gaussian process and the vertical velocity component shows small departures from Gaussianity. The distribution is indeed slightly asymmetric and flatter than the Gaussian for low velocity values. The vertical derivative of the vertical velocity PDFs are shown in figure 29(b). The tails are clearly closer to an exponential than a Gaussian, except the derivative of the vertical velocity for positive values. The PDF of the horizontal velocity gradient, $\unicode[STIX]{x2202}_{z}u_{x}^{\prime }$ , is symmetric and the tails are exponential, i.e. far from Gaussianity, as opposed to the PDF of the vertical velocity gradient, $\unicode[STIX]{x2202}_{z}u_{z}^{\prime }$ , which is not symmetric. Positive events are close to Gaussianity, but negative values of the gradient are on an exponential tail. Exponential tails are also observed in figure 30, where the PDFs of the horizontal and vertical potential vorticity components are displayed. These departures from Gaussianity are associated with intermittency. Pressure PDFs are shown in figure 30(a) for two different proportions of the mixing layer, defined by $0.4\leqslant \overline{c}\leqslant 0.6$ and $0.2\leqslant \overline{c}\leqslant 0.8$ , respectively. These two PDFs are very close to each other, this is a general result which holds for most of the variables. Positive events are Gaussian, while on the negative side, there is an exponential tail. This behaviour is consistent with earlier observations on various turbulent flows (Lesieur Reference Lesieur1997, chap. VI, § 8). Indeed, it occurs in decaying isotropic turbulence (Kalelkar Reference Kalelkar2006), statistically steady turbulence (Cao, Chen & Doolen Reference Cao, Chen and Doolen1999) and in RT Boussinesq turbulence (Schneider & Gauthier Reference Schneider and Gauthier2016c ). As is recalled in this former reference, this asymmetry is related to the fact that vortex centres and pressure minima correspond. However the tail is only exponential as opposed to the Boussinesq case at a Taylor–Reynolds number $Re_{\unicode[STIX]{x1D706}_{zz}}\approx 142$ where the tail is superexponential (Pumir Reference Pumir1994a ; Schneider & Gauthier Reference Schneider and Gauthier2016c ). Figure 30(b) shows the PDFs of the potential vorticity components, $\unicode[STIX]{x1D714}_{x,z}$ . Tails are of exponential type, which is characteristic of the small-scale intermittency (Siggia Reference Siggia1981; She, Jackson & Orszag Reference She, Jackson and Orszag1990; Vincent & Meneguzzi Reference Vincent and Meneguzzi1991, Reference Vincent and Meneguzzi1994; Pumir Reference Pumir1994b ). Concentration skewness and flatness versus the horizontal coordinate $z$ are displayed in figure 31. Skewness takes large values at the mixing layer edges, positive on the light fluid side, where $\overline{c}=0$ , and negative on the heavy fluid side, where $\overline{c}=1$ . The profiles are however slightly skewed, values on the light fluid side are larger than the values of the heavy fluid. This is clearly a variable-density effect (Livescu & Ristorcelli Reference Livescu and Ristorcelli2008). Inside the mixing layer, the skewness is small and the plateau is increasing with time. The flatnesses are close, at least from $t=4.15$ , but not equal to 3, the value of a Gaussian process. This behaviour should lead to a quasi-symmetric PDF. The flatnesses of the concentration fluctuation $x$ - and $z$ -gradients are displayed in figure 32, in semi-log scale. Very high values are reached on the mixing layer edges, while inside the mixing layer, flatness profiles are rather flat, but far away from a Gaussian process. As a result strong departures from Gaussianity may be inferred. Figure 33(a) displays the PDFs of the concentration calculated in the interval $z_{1}\leqslant z\leqslant z_{2}$ , in such a way that $\overline{c}(z_{1})=0.4$ and $\overline{c}(z_{2})=0.6$ and $\overline{c}(z_{1})=0.2$ and $\overline{c}(z_{2})=0.8$ . These results are close to each other as we have previously noted. They are slightly asymmetric and close to a Gaussian process for small-amplitude events. Moreover, the largest interval ( $\overline{c}(z_{1})=0.2$   and   $\overline{c}(z_{2})=0.8$ )  is closer to Gaussianity than the small interval. This is

Figure 31. Mixing statistics. Skewness, $\overline{c^{\prime 3}}/\overline{c^{\prime 2}}^{3/2}$ (a) and flatness, $\overline{c^{\prime 4}}/\overline{c^{\prime 2}}^{2}$ (b), of the concentration fluctuation at six different times, $t=$ 3.25, 4.15, 4.35, 6.85, 7.45 and 8.30. The horizontal black dashed line stands for the Gaussian process flatness.

Figure 32. Mixing statistics. Flatness of the concentration fluctuation gradient versus the vertical coordinate $z$ at six different times, $t=3.25$ , 4.15, 4.35, 6.85, 7.45 and 8.30. (a) Horizontal gradient $\overline{(\unicode[STIX]{x2202}_{x}c^{\prime })^{4}}/\overline{(\unicode[STIX]{x2202}_{x}c^{\prime })^{2}}^{2}$ . (b) Vertical gradient $\overline{(\unicode[STIX]{x2202}_{z}c^{\prime })^{4}}/\overline{(\unicode[STIX]{x2202}_{z}c^{\prime })^{2}}^{2}$ . The horizontal black dashed line stands for the Gaussian process flatness.

Figure 33. Mixing statistics. (a) PDFs of the concentration fluctuations $c^{\prime }$ at time $t=7.45$ , for two different intervals. (b) PDFs of the concentration fluctuation $x$ - and $z$ -gradients, at the same time.

Figure 34. Density statistics. (a) PDFs of the density fluctuations $c^{\prime }$ at time $t=7.45$ , for two different intervals defined by the mean concentration $\overline{c}$ . (b) Density fluctuation $x$ - and $z$ -gradient PDFs, at the same time.

Figure 35. Temperature statistics. Skewness $\overline{T^{\prime 3}}/\overline{T^{\prime 2}}^{3/2}$ (a) and flatness $\overline{T^{\prime 4}}/\overline{T^{\prime 2}}^{2}$ (b) of the temperature fluctuations versus the vertical coordinate $z$ , at six different times, $t=3.25$ , 4.15, 4.35, 6.85, 7.45 and 8.30. The horizontal black dashed line stands for the Gaussian process flatness.

consistent with the skewnesses and flatnesses displayed in figure 31. The PDFs of the horizontal and vertical gradients of the concentration fluctuations are displayed in figure 33 and shown a different behaviour. The PDF of the horizontal gradient $\unicode[STIX]{x2202}_{x}c^{\prime }$ is symmetric but the tails are rather exponential than Gaussian. It means that concentration fluctuations show some organization at intermediate scales. The PDF of the horizontal gradient $\unicode[STIX]{x2202}_{x}c^{\prime }$ is strongly asymmetric and the tails are rather exponential than Gaussian. These tails are piecewise exponential, with a smaller slope for positive values than for the negative values. It means that very large positive values of the vertical gradient concentration are more frequent than the negative ones. To summarize, the concentration statistics in this fully compressible case is not very different from the Boussinesq case (Schneider & Gauthier Reference Schneider and Gauthier2016c ). However in this compressible flow, one has access to the density whose PDFs are displayed in figure 34(a). They have been obtained from the data points that are such that $0.40\leqslant \overline{c}\leqslant 0.60$ and $0.20\leqslant \overline{c}\leqslant 0.80$ at time $t\approx 7.45$ . As opposed to concentration PDFs, they are rather skewed and positive events are clearly favoured with respect to negative ones. The PDFs of the horizontal and vertical gradients, $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D70C}^{\prime }$ and $\unicode[STIX]{x2202}_{z}\unicode[STIX]{x1D70C}^{\prime }$ , of the density fluctuations are displayed in this same figure 34(b). They show a similar behaviour to the concentration PDFs, i.e. the $\unicode[STIX]{x2202}_{x}c^{\prime }$ -PDF is rather symmetric with exponential tails, while the $\unicode[STIX]{x2202}_{z}c^{\prime }$ -PDF is strongly asymmetric with exponential tails. Large-amplitude positive events are more frequent than the negative. Negative density and temperature gradients are thus more probable. The largest gradients are positive, but rare. Temperature statistics is now detailed and as it has already been stressed, temperature behaviour is important for the applications in high energy density physics, particularly in inertial confinement fusion and in astrophysics. Indeed the temperature field is modified in a RT mixing layer, even with an isothermal initial state. We will show that the temperature statistics follows closely the mixing statistics. Temperature skewness, $\overline{T^{\prime 3}}/\overline{T^{\prime 2}}^{3/2}$ , and flatness, $\overline{T^{\prime 4}}/\overline{T^{\prime 2}}^{2}$ , versus the horizontal coordinate $z$ are displayed in figure 35. Skewnesses take large values at the turbulent layer edges, and as the concentration skewness, it is positive on the light fluid side and negative on the heavy fluid side. However, as opposed to the concentration skewness, it is strongly skewed, positive values are much larger than negative ones, in absolute values. The plateau of almost zero skewness is increasing with time. Flatness is given in figure 35(b). Temperature flatness is similar to concentration flatness and they are close to a Gaussian process only at the very interior of the mixing layer, during the RT regime. Consequently, the temperature PDF is also skewed and asymmetric. The flatnesses of the temperature fluctuation $x$ - and $z$ -gradients versus the vertical coordinate $z$ , $\overline{(\unicode[STIX]{x2202}_{i}T^{\prime })^{4}}/\overline{(\unicode[STIX]{x2202}_{i}T^{\prime })^{2}}^{2}$ , $i=x,z$ , at six different times are represented in figure 36. Asymmetry is clearly observed in both $x$ - and $z$ -directions, flatnesses are stronger in the light fluid side. Departures from Gaussianity are also visible and more marked in the vertical direction. These skewness and flatness characteristics are related to the shape of the PDFs, which are displayed in figure 37. In brief, temperature PDFs are close to Gaussianity, although rare large-amplitude events deviate from this Gaussianity. Temperature fluctuation $x$ - and $z$ -gradient PDFs are clearly non-Gaussian, they are skewed, with exponential tails. Large positive-amplitude events are more frequent than negative events of the same amplitude. One observes a strong coherence between the three PDFs of the vertical gradients of concentration, density and temperature (see figures 33–35). Moreover, the density and temperature $z$ -gradient maxima are slightly shifted toward negative values. While these $z$ -gradient PDFs bear some resemblance with PDFs obtained in Pumir (Reference Pumir1994b ) for a passive scalar, they are are much more skewed. However, the present configuration of a two-component mixing with variable temperature, linked with density and concentration through an EOS, is quite different from the canonical situation of a passive scalar in an incompressible flow. It has been shown by Ashurst et al. (Reference Ashurst, Kerstein, Kerr and Gibson1987) and for various turbulent flows that the scalar gradient is preferentially aligned with the strain rate eigenvector corresponding to the smallest – the third – eigenvalue, while vorticity is preferentially aligned with the second eigenvector. As a result, vorticity and scalar gradient are mostly orthogonal so that the cosine values of the angle between vorticity and scalar gradient, i.e. $\cos \unicode[STIX]{x1D703}=\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c/|\unicode[STIX]{x1D714}||\unicode[STIX]{x1D735}c|$ , are mostly zero, where $c$ stands for a scalar field. In other words, the PDF of the cosine

Figure 36. Temperature statistics. Flatness of the temperature fluctuation $x$ - and $z$ gradients versus the vertical coordinate $z$ , at six different times, $t=$ 3.25, 4.15, 4.35, 6.85, 7.45 and 8.30. (a) Horizontal gradient $\overline{(\unicode[STIX]{x2202}_{x}T^{\prime })^{4}}/\overline{(\unicode[STIX]{x2202}_{x}T^{\prime })^{2}}^{2}$ . (b) Vertical gradient $\overline{(\unicode[STIX]{x2202}_{z}T^{\prime })^{4}}/\overline{(\unicode[STIX]{x2202}_{z}T^{\prime })^{2}}^{2}$ . The horizontal black dashed line stands for the Gaussian process flatness.

Figure 37. Temperature statistics. (a) PDFs of the temperature fluctuations $T^{\prime }$ at time $t=7.45$ , for two different intervals defined by the mean concentration, $\overline{c}$ . (b) PDFs of the temperature fluctuation $x$ - and $z$ -gradient, at the same time.

of this angle is peaked near $\unicode[STIX]{x1D703}=0$ . This is true for the mixing of a passive scalar in the presence of a mean gradient (Pumir Reference Pumir1994b ) or in a RT turbulent mixing layer, within the Sandoval model (Cabot & Zhou Reference Cabot and Zhou2013). The PDFs of the cosine of the angle between vorticity and temperature, density and concentration gradient are displayed in figure 38(a) for the $Sr6\text{-}Re6\times 10^{4}$ simulation. These three PDFs are almost superimposed, which shows that these three scalars are correlated to the vorticity field and at the same high level. The PDF of the cosine of the angle between vorticity and pressure gradient is also displayed (figure 38 b), which shows that these two quantities are also correlated, but the maximum is smaller than the one obtained for the concentration–vorticity correlation. Moreover, these values are significantly smaller than those obtained from a homogeneous isotropic turbulence at moderate Reynolds number (Kalelkar Reference Kalelkar2006, figure 9). Figure 38(a) also displays the PDF of the cosine of the angle between concentration and density gradient. It is strongly peaked near 1, i.e. for $\unicode[STIX]{x1D703}\approx 0$ . These two gradients are thus essentially aligned. Finally, the PDF of the cosine of the angle between pressure and density gradient is also represented in the same figure. It shows clearly the absence of correlation between these two quantities, except for $\unicode[STIX]{x1D703}=0$ and $\unicode[STIX]{x03C0}$ . These two gradients are thus almost randomly distributed, the alignment and anti-alignment are slightly more frequent.

Figure 38. Statistics. PDFs of the cosine of the angle between various quantities at time $t=7.45$ , the turbulent kinetic energy maximum. (a) Superimposition of the PDFs of the cosine of the angle between vorticity and temperature (blue), density (green) and concentration (red) gradient, respectively. (b) PDFs of the cosine of the angle between vorticity and pressure gradient (blue) and the gradients of density and concentration (green) and density and pressure (red).

9 Flow visualization

RT turbulent flow visualization has been provided by several authors (efluids 2017). In Schneider & Gauthier (Reference Schneider and Gauthier2016b ), isosurfaces and slices of concentration, vorticity, $Q$ -criterion, turbulent kinetic energy, Taylor–Reynolds number, dissipation, pressure and temperature have been displayed for the Boussinesq, anelastic and fully compressible models. To illustrate the flow detailed in this article, temperature and helicity isosurfaces are displayed in figure 39, in addition to the palinstrophy isosurface displayed in figure 12. Temperature appears in blue in the top layer, since the heavy fluid is cooled and the light fluid is heated. Mushroom-like patterns are also clearly visible on the temperature field. The helicity isosurfaces coloured by the concentration are displayed in figure 39(b), where small scales participating in the mixing are also visible.

Figure 39. Visualization. (a) Temperature isosurfaces (top view) $T=0.95$ (blue) and $T=1.05$ (red). Temperature at rest is $T=1$ . Time is $t=4.6$ . (b) Helicity isosurfaces for the two values, ${\mathcal{H}}=\pm 0.20$ , coloured by the concentration, at time $t=4.35$ , top view.

10 Summary and conclusions

We have presented a detailed analysis of a large-scale DNS performed with the full NSEs for two Newtonian fluids, in which the initial equilibrium state is strongly stratified. These results have been obtained with a pseudo-spectral (Chebyshev–Fourier–Fourier) numerical method with an auto-adaptive Chebyshev multidomain method. The spatial resolution is varied during the calculation from $(9\times 64)\times 576^{2}\approx 191M$ to $(9\times 100)\times 1000^{2}=900M$ collocation points. In such a stratified RT configuration, the flow is unsteady. After the linear regime and the transition to turbulence, the RT regime is qualitatively close to the classical RT growth. However, quantitatively, this regime is affected by the stratification, since the mixing layer starts to smooth the density jump so that, at some instant, physical variables start to decay, but at different times and at different rates. In particular, turbulence erases the density jump and leads to the mixing layer growth termination. The final step is a freely decaying turbulence regime in a stable stratification. Various compressibility effects are involved in the present flow, static and dynamic. The stratification is a variable-density effect. This static compressibility has a little influence on the linear regime but a strong influence on the nonlinear behaviour since it leads to the termination of the turbulent mixing layer growth. Compressibility effects due to the mixing or non-zero velocity divergence effects are dynamic. There are also acoustic effects or EOS effects, due to the finite speed of sound. Indeed a strong acoustic production is observed, but the Mach number and the velocity divergence are rather small.

Two time-variable Atwood numbers have been used. A large-scale global Atwood number is built from the mean density profile, while a small-scale local Atwood number is provided by the r.m.s.-density fluctuations (Cook et al. Reference Cook, Cabot and Miller2004). The large-scale Atwood number decays quickly in the RT regime and vanishes at the time where the small-scale Atwood number reaches its maximum. As a result, some local buoyancy are still present in the freely decaying regime. Beyond some point in this regime, the decaying turbulence homogenizes the mixing and the concentration profile flattens. In the freely decaying regime, the mean value of the spike thickness saturates while bubble thickness still grows at a very small rate. None of these thicknesses follow the $t^{2}$ -scaling, since the density length scale, $L_{\unicode[STIX]{x1D70C}H,L}$ , is of the same order as the mixing thickness. One also observes that the mixing region located in the heavy fluid side is significantly cooled, while the mixture located in the light fluid is heated. In addition to the definition of two effective Atwood numbers, a local condition for the stability of a compressible fluid is given by the gradient of entropy. It shows that the instability region is located in the upper side of the mixing layer.

The molecular mixing fraction reaches the value $\langle \unicode[STIX]{x1D703}\rangle _{\unicode[STIX]{x1D6FD}}=0.80$ at the end of the RT regime, close to the value obtained from the Boussinesq approximation (Schneider & Gauthier Reference Schneider and Gauthier2016c ). This is within the interval obtained in experiments with incompressible fluids, i.e. 0.75–0.80 (Linden et al. Reference Linden, Redondo and Youngs1994; Dalziel et al. Reference Dalziel, Linden and Youngs1999; Ramaprabhu & Andrews Reference Ramaprabhu and Andrews2004). The final value, $\langle \unicode[STIX]{x1D703}\rangle _{\unicode[STIX]{x1D6FD}}=0.99$ , is close to homogeneity. Since the mixing layer stop to grow quite early, no more pure fluids enter in the turbulent layer. However turbulence and vortical motions continuously mix the two fluids and increases the homogeneity.

The vertical baroclinic vorticity production is non-zero, as opposed to the incompressible case, but is much smaller than the horizontal components, since these components involve vertical density and pressure gradients. The influence of the acoustic waves are also clearly seen on the vertical vorticity component during the FD regime. The helicity evolution and the three source term behaviours have been investigated and dissipation appears to be dominant. The baroclinic and pressure source terms – typical of compressible flows – are of same order of magnitude. This helicity evolution shares some resemblance with the kinetic energy and the Mach number evolution. Palinstrophy evolution of the three simulations has also been described and compared to each other.

Comparisons between numerical results and classical power laws in the inertial range are often difficult, since inertial ranges are often narrow. This is the case in this simulation since the Taylor–Reynolds number maximum is modest. Moreover, this flow is unsteady and compressibility effects, i.e. variable density and acoustics, are present, even if the latter is weak. One-dimensional spectra of various quantities have nevertheless been displayed. It turns out that the concentration and temperature spectra are very similar to each other, which shows that concentration and temperature behave and evolve closely. The agreement of the power spectra of the concentration and the temperature with the classical $k^{-17/3}$ -scaling is encouraging, although it is observed on a very narrow spectral region.

Anisotropy is an important issue in RT turbulent mixing layers. This characteristic have been investigated with the mean value of the Reynolds stress tensor and a spectral anisotropy indicator. It is confirmed that the anisotropy is very high in the RT regime. At the r.m.s.-density maximum, the mean Reynolds stress tensor reaches 0.55, which means that $88\,\%$ of the turbulent kinetic energy is in the third direction, while at the maximum of the kinetic energy, $64\,\%$ of the turbulent kinetic energy is in the third direction. Velocity field spectral anisotropy shows a persistent anisotropy at all scales, as opposed to the Boussinesq case, where intermediate scales are clearly isotropic, while small scales are anisotropic. Concentration and temperature spectral anisotropy are very close to each other and exhibit intermediate-scale isotropy and small-scale anisotropy, as it has been observed in the Boussinesq turbulence. This isotropy/anisotropy is present from short times to late times.

Favre-averaged equations – turbulent kinetic energy, r.m.s. density, mass flux and internal energy equations – have also been used to analyse the numerical data. The equilibrium, i.e. the ratio between turbulence production and dissipation has been investigated. It turns out that such a turbulent RT compressible mixing layer in a stratified configuration is never in equilibrium. Instead strong departures from equilibrium are observed during the RT and the FD regimes. The source terms of the r.m.s.-density and mass flux $\overline{\unicode[STIX]{x1D70C}^{\prime }u_{i}^{\prime }}$ equations have been studied. Concerning the r.m.s.-density equation, the source term proportional to the vertical density gradient is one order of magnitude larger than the term proportional to the velocity gradient. Concerning the mass flux equation, the term proportional to the pressure gradient is dominant and always positive. Velocity-gradient production is always negligible, while density-gradient production is small, and becomes quickly negative.

In this compressible flow, the three Kovásznay modes, namely, vorticity, entropy and acoustic are present. These two latter modes are actually compressibility effects since the simplest RT model, the Boussinesq approximation does not contain a temperature. The r.m.s. density, pressure, temperature and concentration have been first displayed. The vorticity mode has been analysed with the vorticity, helicity and palinstrophy behaviours. It has been assumed that the acoustic mode is given by the pressure fluctuations, and the entropy mode is thus deducted. It turns out that the density and temperature entropic parts are one or two orders of magnitude larger than the acoustic part. The density entropy mode grows in the RT regime and start to decay right after, while the temperature entropy mode grows much more tardily.

A statistical study has been conducted. Skewnesses, flatnesses and PDFs have been plotted and commented. Some of the conclusions drawn in the Boussinesq case apply in this weakly compressible turbulent flow. Indeed large scales have a Gaussian behaviour, while intermediate scales, associated with the gradients of various quantities, show departures from Gaussianity, i.e. PDF are strongly asymmetric and wings are exponential. It has been confirmed that the PDF of the cosine of the angle between vorticity and concentration gradient is peaked near $\unicode[STIX]{x1D703}=0$ . Moreover, the PDFs of the angle between vorticity and concentration, density and temperature gradients are superimposed, so that these quantities are mostly aligned.

Temperature field appears to be the slave of the mixing. Indeed the thermal layer has the same thickness than the mixing layer. Concentration and temperature spectra are very close to each other, along the time. Anisotropy spectra of density and temperature are also very close to each other. Concentration and temperature PDFs are very similar. Moreover the correlation coefficients $R_{\unicode[STIX]{x1D70C}}(\unicode[STIX]{x1D70C},p)$ and $R_{T}(T,p)$ have very similar behaviour and the density–temperature correlation coefficient $S(\unicode[STIX]{x1D70C},T)$ shows that these two quantities are essentially correlated or anti-correlated. However there is a significant time lag between the density and temperature evolution and, in particular, the occurrence of their maxima in time. Acoustics has a little influence on the concentration and temperature fields.

Finally, let us recall that the Boussinesq model leads to well-defined results, i.e. self-similar behaviours and universal scalings. This approximation may thus be investigated with only one large-scale simulation. By contrast, each compressible case with a stratified initial equilibrium state is specific. A large number of simulations is thus required to accumulate results and to obtain firm conclusions. In particular, different equations of state, various initial temperature profiles and boundary conditions, among others, have to be investigated.