1 Introduction
The irreversible mixing at a molecular level of two fluids of different densities
$\unicode[STIX]{x1D70C}_{2}^{\ast }>\unicode[STIX]{x1D70C}_{1}^{\ast }$
is a fluid dynamical process of great fundamental interest and practical importance, especially when the fluids are turbulent. Such turbulent mixing flows occur in many different circumstances. A particularly important class arises when the buoyancy force associated with the effects of statically unstable variations in fluid density in a gravitational field actually drives both the turbulence and the ensuing mixing itself. Such flows, commonly referred to as ‘Rayleigh–Taylor instability’ (RTI) flows due to the form of the initial linear instability (Rayleigh Reference Rayleigh1900; Taylor Reference Taylor1950), have been very widely studied (see Sharp Reference Sharp1984; Youngs Reference Youngs1984, Reference Youngs1989; Glimm et al.
Reference Glimm, Grove, Li, Oh and Sharp2001; Dimonte et al.
Reference Dimonte, Youngs, Dimits, Weber, Marinak, Wunsch, Garasi, Robinson, Andrews and Ramaprabhu2004; Dimotakis Reference Dimotakis2005; Hyunsun et al.
Reference Hyunsun, Hyeonseong, Yan and Glimm2008; Lee et al.
Reference Lee, Jin, Yu and Glimm2008; Andrews & Dalziel Reference Andrews and Dalziel2010), not least because of their relevance in astrophysics (Cabot & Cook Reference Cabot and Cook2006) and fusion (Petrasso Reference Petrasso1994).
A key characteristic of RTI flows is that the developing turbulence is not driven by some external forcing mechanism, but rather is supplied with kinetic energy by the conversion of ‘available’ potential energy stored in the initial density field. This kinetic energy naturally drives turbulent disorder and a cascade to small scales, with an attendant increase in the dissipation rate of kinetic energy. Such small scales also lead to ‘filamentation’, i.e. enhanced surface area of contact between the two miscible fluids and, crucially, substantially enhanced gradients in the density field, which thus also leads to irreversible mixing, and hence modification in the density distribution. There has been an explosion in interest in investigating the ‘efficiency’ of this mixing, i.e., loosely, the proportion of the converted available potential energy that leads to irreversible mixing, as opposed to viscous dissipation (see the recent review of Tailleux (Reference Tailleux2013)), although the actual definition and calculation of the efficiency are subtle and must be performed with care – see, for example, Davies-Wykes & Dalziel (Reference Davies-Wykes and Dalziel2014) for further discussion.
Nevertheless, there is accumulating evidence that buoyancy-driven turbulence is particularly efficient in driving mixing (Lawrie & Dalziel Reference Lawrie and Dalziel2011; Davies-Wykes & Dalziel Reference Davies-Wykes and Dalziel2014), and certainly more efficient than externally forced turbulent flow. This evidence poses the further question of whether there are distinguishing characteristics of the buoyancy-driven turbulent flow that are different from the flow associated with an external forcing; in particular, whether these characteristics can be identified as being responsible for the enhanced and efficient mixing.
The situation is further complicated by the observation that, even when the two fluids undergoing mixing are themselves incompressible, since molecular mixing generically changes the specific volume of the mixture, the velocity fields of such ‘variable-density’ (VD) flows (following the nomenclature suggested by Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007)) are, in general, not divergence-free. This is certainly the case when the two densities are sufficiently different that the Boussinesq approximation may not be applied. Commonly, the Boussinesq approximation is applied when the Atwood number
$At$
, defined as
$$\begin{eqnarray}At=\frac{\unicode[STIX]{x1D70C}_{2}^{\ast }-\unicode[STIX]{x1D70C}_{1}^{\ast }}{\unicode[STIX]{x1D70C}_{2}^{\ast }+\unicode[STIX]{x1D70C}_{1}^{\ast }},\end{eqnarray}$$
is small; i.e.
$At\ll 1$
. However, as discussed in detail in Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), non-Boussinesq effects may occur when gradients in the density field become large. Following Cook & Dimotakis (Reference Cook and Dimotakis2001) and Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), the composition density
$\unicode[STIX]{x1D70C}^{\ast }(\boldsymbol{x},t)$
of a mixture of two constant fluid densities
$\unicode[STIX]{x1D70C}_{1}^{\ast }$
and
$\unicode[STIX]{x1D70C}_{2}^{\ast }$
(
$\unicode[STIX]{x1D70C}_{2}^{\ast }>\unicode[STIX]{x1D70C}_{1}^{\ast }$
) is expressed in dimensionless form by
$$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D70C}^{\ast }(\boldsymbol{x},t)}=\frac{Y_{1}(\boldsymbol{x},t)}{\unicode[STIX]{x1D70C}_{1}^{\ast }}+\frac{Y_{2}(\boldsymbol{x},t)}{\unicode[STIX]{x1D70C}_{2}^{\ast }},\end{eqnarray}$$
where the
$Y_{i}(\boldsymbol{x},t)$
(
$i=1,2$
) are the mass fractions of the two fluids, and
$Y_{1}+Y_{2}=1$
. Equation (1.2) shows that the composition density
$\unicode[STIX]{x1D70C}^{\ast }$
is bounded by
$\unicode[STIX]{x1D70C}_{1}^{\ast }\leqslant \unicode[STIX]{x1D70C}^{\ast }(\boldsymbol{x},t)\leqslant \unicode[STIX]{x1D70C}_{2}^{\ast }$
. Assuming that there is Fickian diffusion, the mass transport equations for the two species are
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}Y_{i})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }Y_{i}\boldsymbol{u})=Pe_{0}^{-1}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}Y_{i}),\end{eqnarray}$$
where
$Pe_{0}$
is the Péclet number; the dimensionless Reynolds, Schmidt and Péclet numbers are defined in table 1. Since the specific volume
$1/\unicode[STIX]{x1D70C}^{\ast }$
changes due to mixing, a non-zero divergence is induced in the velocity field, which affects the conventional continuity equation for mass conservation,
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}^{\ast }+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }\boldsymbol{u})=0,\end{eqnarray}$$
which is derived from the sum over the two species in (1.3). In appendix A, it is shown that the full equations for
$\unicode[STIX]{x1D70C}^{\ast }$
derived from (1.2) and (1.3) take the form (see (A 6) and (A 7))
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D70C}^{\ast }=Pe_{0}^{-1}\unicode[STIX]{x1D70C}^{\ast }\unicode[STIX]{x1D6E5}(\ln \unicode[STIX]{x1D70C}^{\ast }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=-Pe_{0}^{-1}\unicode[STIX]{x1D6E5}(\ln \unicode[STIX]{x1D70C}^{\ast }). & \displaystyle\end{eqnarray}$$
As discussed in Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), the Boussinesq approximation leads to the requirement that the velocity field is divergence-free,
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0$
, in which case the mass conservation equation becomes
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}^{\ast }+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }=Pe_{0}^{-1}\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D70C}^{\ast }.\end{eqnarray}$$
This relies on the requirement that the second (nonlinear) term on the right-hand side of (1.5) can be ignored compared with the first term, i.e. that
$$\begin{eqnarray}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }|^{2}\ll \unicode[STIX]{x1D70C}^{\ast }|\unicode[STIX]{x1D6E5}\unicode[STIX]{x1D70C}^{\ast }|.\end{eqnarray}$$
As noted by Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), this condition may be violated if substantial gradients develop in the density field. It is not a priori clear, even when the Atwood number is very small, that the non-divergence-free nature of the velocity field qualitatively changes the properties of the turbulent flow in ways that are significant to the mixing, and specifically whether regions may develop in the flow where the condition (1.8) is violated. This issue can be explored by careful numerical simulation, as reviewed by Livescu (Reference Livescu2013), with a key observation (see Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007) for more details) being that the pressure distribution is substantially modified by non-Boussinesq effects. Furthermore, the central role played by intermittency and anisotropy, as discussed in Livescu & Ristorcelli (Reference Livescu and Ristorcelli2008), suggests that it would be instructive to focus carefully on the time-dependent evolution of nonlinearity within such buoyancy-driven variable-density flows using the full equations (1.5) and (1.6) without the approximation (1.8).
Simulation parameters.
Indeed, using a normalization density
$\unicode[STIX]{x1D70C}_{0}^{\ast }=(\unicode[STIX]{x1D70C}_{1}^{\ast }+\unicode[STIX]{x1D70C}_{2}^{\ast })/2$
, an interesting observation is that with
$$\begin{eqnarray}\unicode[STIX]{x1D703}(\boldsymbol{x},t)=\ln \unicode[STIX]{x1D70C},\quad \unicode[STIX]{x1D70C}=\frac{\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x1D70C}_{0}^{\ast }},\end{eqnarray}$$
(1.5) becomes a diffusion-like equation coupled to a
$\boldsymbol{u}$
-field whose divergence is linked to two derivatives in
$\unicode[STIX]{x1D703}$
,
$$\begin{eqnarray}(\unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D703}=Pe_{0}^{-1}\unicode[STIX]{x0394}\unicode[STIX]{x1D703},\quad \text{with}~\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=-Pe_{0}^{-1}\unicode[STIX]{x0394}\unicode[STIX]{x1D703}.\end{eqnarray}$$
While
$\unicode[STIX]{x1D70C}^{\ast }$
itself is bounded above and below, the two gradients
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$
and
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
are of much greater interest. Mathematically, an understanding of the properties of these gradients may reveal something about the regularity or otherwise of solutions of the governing equations. Conversely, physically, the magnitudes of these gradients are central to the rate of irreversible mixing within the flow.
Following Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), there is another way of looking at the growth of these gradients. Consider the equation for
$\unicode[STIX]{x1D703}$
and introduce a new velocity field
$\boldsymbol{v}=\boldsymbol{u}+Pe_{0}^{-1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
. The Hopf–Cole-like transformation
$\unicode[STIX]{x1D703}=\ln \unicode[STIX]{x1D70C}$
in (1.9) then leads to an exact cancellation of the nonlinear terms in (1.10) to give
$$\begin{eqnarray}(\unicode[STIX]{x2202}_{t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D70C}=Pe_{0}^{-1}\unicode[STIX]{x0394}\unicode[STIX]{x1D70C},\quad \text{with}~\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0.\end{eqnarray}$$
This is the linear advection diffusion equation driven by a divergence-free velocity field. It should be noted that
$\unicode[STIX]{x1D74E}=\text{curl}\,\boldsymbol{u}=\text{curl}\,\boldsymbol{v}$
. However, the fact that
$\boldsymbol{v}$
is actually an (explicit) function of
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
makes (1.11) less simple than it first appears. Nevertheless, this equation provides a hint as to how we might look at the dynamics in a descriptive way. Consider a one-dimensional horizontal section through a rightward moving wave of
$\unicode[STIX]{x1D70C}$
at a snapshot in time. In the frame of the advecting velocity
$\boldsymbol{u}$
the relevant component of
$\boldsymbol{v}$
is greater on the back face of any part of the wave (where
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}>0$
) than on the front face (where
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}<0$
). Thus, in the advecting frame, (1.11) implies that not only is there the usual advection and diffusion but also a natural tendency for the back of a wave to catch up with the front, thus leading to a natural and inevitable steepening of
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$
.
While our explanation of this steepening process is heuristic, it is consistent with experimental observations, and raises several questions. To explore these, we have undertaken a reanalysis of a data set of Livescu, arising from the simulation of a buoyancy-driven flow very similar to that reported in Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), which is freely available at the Johns Hopkins turbulence database (JHTDB). Using this reanalysis, there are three central questions which we wish to explore as the primary aims of this paper.
-
(i) First, can the growth of gradients in the density field be bounded or controlled in any meaningful way? Such bounds could yield valuable insights into the structure and regularity of the density field and the uniform validity of the Boussinesq approximation for flows with
$At\ll 1$
, which may explain the ‘efficiency’ of mixing associated with buoyancy-driven turbulence? -
(ii) Alternatively, if the growth of density gradients fails to saturate and instead accelerates by drawing on available kinetic energy, is it fast enough for a singularity to form in a finite time? The rapid growth observed by Luo & Hou (Reference Luo and Hou2014a ,Reference Luo and Hou b ) in the vorticity field of a symmetric incompressible Euler flow in a cylinder (with
$\boldsymbol{u}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=0$
on the boundary) provides strong evidence of a singularity in that system. It is possible that buoyancy-driven variable-density turbulence is another example, despite the difference in boundary conditions and the presence of diffusion. However, it ought to be pointed out that neither the data nor the methods used here are sufficient to answer this question definitively. -
(iii) Does the nonlinear ‘depletion’ (a concept we define below) in three-dimensional Navier–Stokes turbulence survive when it is coupled to the density field? In tandem with question (i), is there any nonlinear depletion in this density field?
The method of analysis used to address these questions is based on that used to analyse multiple data sets generated from the incompressible three-dimensional Navier–Stokes equations (Donzis et al.
Reference Donzis, Gibbon, Kerr, Gupta, Pandit and Vincenzi2013; Gibbon et al.
Reference Gibbon, Donzis, Kerr, Gupta, Pandit and Vincenzi2014; Gibbon Reference Gibbon2015) on a cubical domain
$[0,L]^{3}$
. Here, we review this method for completeness. The method involves taking higher
$L^{p}$
-norms of the vorticity
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{m}=\left(L^{-3}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D74E}|^{2m}\,\text{d}V\right)^{1/2m},\quad \text{for}~m\geqslant 1,\end{eqnarray}$$
which are, in effect, higher moments of the enstrophy field, each having the dimensions of a frequency. Symmetry considerations in the three-dimensional Navier–Stokes equations suggest that the following scaling is appropriate (
$\unicode[STIX]{x1D71B}_{0}=\unicode[STIX]{x1D708}L^{-2}$
):
$$\begin{eqnarray}D_{m}=(\unicode[STIX]{x1D71B}_{0}^{-1}\unicode[STIX]{x1D6FA}_{m})^{\unicode[STIX]{x1D6FC}_{m}},\quad \text{where}~\unicode[STIX]{x1D6FC}_{m}=\frac{2m}{4m-3}.\end{eqnarray}$$
In most theoretical analyses, it is impossible to avoid gradients of
$\unicode[STIX]{x1D74E}$
in expressions for the higher moments
$\unicode[STIX]{x1D6FA}_{m}$
, thus causing great difficulties with closure. It has been shown in Gibbon (Reference Gibbon2015) that the sequence of
$D_{m}$
can be connected with
$D_{1}$
in the following way (
$m\geqslant 2$
):
$$\begin{eqnarray}D_{m}=C_{m}D_{1}^{A_{m,\unicode[STIX]{x1D706}}}\quad \text{and}\quad A_{m,\unicode[STIX]{x1D706}}=\frac{\unicode[STIX]{x1D706}_{m}(t)(m-1)+1}{4m-3},\end{eqnarray}$$
where the set of exponents
$\{\unicode[STIX]{x1D706}_{m}(t)\}$
, subject to
$1\leqslant \unicode[STIX]{x1D706}_{m}\leqslant 4$
, are time-dependent, and the
$C_{m}$
are a set of positive constants. One way of explaining the use of the relation in (1.14) is as follows. Instead of using higher
$L^{2m}$
-norms in the sequence
$\{D_{1};D_{m}(t)\}$
, one examines the sequence
$\{D_{1};\unicode[STIX]{x1D706}_{m}(t);C_{m}\}$
. Thus, the enstrophy
$D_{1}$
is taken as the main variable and the exponents
$\{\unicode[STIX]{x1D706}_{m}(t)\}$
are then monitored numerically. Following how these exponents vary in time thus captures how the scaled moments
$D_{m}$
vary in time. The lower bound
$\unicode[STIX]{x1D706}_{m}\geqslant 1$
comes from the fact that
$\unicode[STIX]{x1D6FA}_{1}\leqslant \unicode[STIX]{x1D6FA}_{m}$
, which can be shown to be equivalent to
$D_{1}^{\unicode[STIX]{x1D6FC}_{m}/2}\leqslant D_{m}$
. The case
$\unicode[STIX]{x1D706}_{m}=1$
lies at the lower bound. At
$\unicode[STIX]{x1D706}_{m}=4$
,
$A_{m}=1$
, and so
$D_{m}=C_{m}D_{1}$
in (1.14). In this situation, it is apparent that the Navier–Stokes nonlinearity is fully saturated. Computations in Donzis et al. (Reference Donzis, Gibbon, Kerr, Gupta, Pandit and Vincenzi2013) and Gibbon et al. (Reference Gibbon, Donzis, Kerr, Gupta, Pandit and Vincenzi2014) have shown that the
$\unicode[STIX]{x1D706}$
-parameter used there lies in the range
$1.15\leqslant \unicode[STIX]{x1D706}\leqslant 1.5$
. It is this that we refer to as ‘nonlinear depletion’, because the higher-order appropriately scaled moments
$D_{m}$
do not ‘saturate’ nonlinearly.
As shown in § 3.1, in the buoyancy-driven flow considered here, it is possible to define not only equivalent higher
$L^{p}$
-norms of the vorticity
$\unicode[STIX]{x1D6FA}_{m}$
and hence
$D_{m}$
to those defined in (1.12) and (1.13) (which we relabel as
$\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D706}}$
and
$D_{m,\unicode[STIX]{x1D714}}$
to make the dependence on vorticity explicit), but also equivalent norms of
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
, which we refer to as
$\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}}$
and
$D_{m,\unicode[STIX]{x1D703}}$
. We can then investigate the equivalent nonlinear depletion in the density gradients through considering the rate of growth of
$D_{1,\unicode[STIX]{x1D703}}$
and then monitoring the evolution of the exponents
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
in the nonlinear growth terms, while also taking into account
$D_{1,\unicode[STIX]{x1D714}}$
. While it is not currently possible to answer definitively the question posed about blow-up of the density gradients at intermediate times, before diffusive mixing might be presumed to smooth out the density distribution, this paper is meant to set the scene as a way of analysing and interpreting future computations.
The rest of the paper is organized as follows. In § 2, we describe in detail the properties of the simulation data set that we reanalyse, and we then present the results of this reanalysis in § 3. Finally, we draw our conclusions in § 4.
2 Description of the database
As noted in the introduction, we use the JHTDB (Livescu et al. Reference Livescu, Canada, Kanov, Burns and Pulido2014), a publicly available direct numerical simulation (DNS) database, to explore nonlinear depletion in buoyancy-driven turbulence. We have used the homogeneous buoyancy-driven turbulence data set for our study. For more information, please see http://turbulence.pha.jhu.edu/datasets.aspx.
The equations used for this problem are the miscible two-fluid incompressible Navier–Stokes equations given by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}^{\ast }+(\unicode[STIX]{x1D70C}^{\ast }u_{j})_{,j}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}^{\ast }u_{i})+(\unicode[STIX]{x1D70C}^{\ast }u_{i}u_{j})_{,j}=-p_{,i}+\unicode[STIX]{x1D70F}_{ij,j}+\frac{1}{Fr^{2}}\unicode[STIX]{x1D70C}^{\ast }g_{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle u_{j,j}=-\frac{1}{Re_{0}Sc}(\ln \unicode[STIX]{x1D70C}^{\ast })_{,jj}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70F}_{ij}=\unicode[STIX]{x1D70C}^{\ast }\,Re_{0}^{-1}\left(u_{i,j}+u_{j,i}-{\textstyle \frac{2}{3}}\unicode[STIX]{x1D6FF}_{ij}u_{k,k}\right), & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}^{\ast }$
is the non-dimensional density of the mixture.
The individual densities of the two components,
$\unicode[STIX]{x1D70C}_{1}^{\ast }$
and
$\unicode[STIX]{x1D70C}_{2}^{\ast }$
, are constant, but due to changes in mass fractions of each species, the density of the mixture can change (1.2). For this reason, the divergence of velocity is dependent on the density, as seen in (2.3). The value of the Atwood number,
$At$
, which characterizes the density difference, is 0.05 and represents a small departure from the Boussinesq approximation. Some of the other important simulation parameters are displayed in table 1, where
$L_{0}$
is the non-dimensionalization length,
$U_{0}$
is the reference velocity scale,
$\unicode[STIX]{x1D707}_{0}$
is the dynamic viscosity and
$D$
is the mass diffusivity.
In the simulation run by Livescu and colleagues, the fluids are initialized as random blobs with periodic boundaries in each direction and an initial diffusion layer at the interface. At sufficiently late time, the statistically homogeneous turbulent flow generated by such conditions resembles the interior of the mixing layer (away from the edge effects) of the Rayleigh–Taylor instability at the turbulent stage (see Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007) for further discussion).
The inhomogeneities in the transport terms are important only at the edge and, thus, it is safe to assume that the homogeneous simulation data under consideration describe the core of a fully developed mixing layer. Eventually, the turbulent behaviour dies out as the fluids become sufficiently mixed at the molecular level.
The variable-density version of the petascale CFDNS code (Livescu et al.
Reference Livescu, Mohd-Yusof, Petersen and Grove2009) was used with an ideal-gas equation of state to carry out the DNS on
$1024^{3}$
grid points (for more information on a similar numerical study, refer to Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007)). The solution algorithm in the code is based on the pseudospectral method. The pressure is computed from the Poisson equation using the projection method. Time integration is achieved using the third-order predictor–corrector Adams–Bashforth–Moulton method. The resultant high-resolution data have been stored as a sequence of 1015 files, ranging from
$t=0$
to
$t=40.56$
, each representing time-step data for
$32^{3}$
spatial points.
The velocity gradients in the database are calculated as a post-processing step using a fourth-order central finite differencing approximation. In general, for calculation of the gradients or the state variables at a particular spatial location between the stored grid points, either a fourth-order spatial interpolation (for gradients) or a sixth-order Lagrangian interpolation (for state variables) is used. To get the temporal values other than those that are stored in the database, a piecewise cubic Hermite interpolation is employed. More information on the specific data set functions used for interpolation and calculation of the gradients is given in appendix B.
3 Results
3.1 Definitions
It is clear from (1.2) that the composition density
$\unicode[STIX]{x1D70C}^{\ast }$
is bounded by
$\unicode[STIX]{x1D70C}_{1}^{\ast }\leqslant \unicode[STIX]{x1D70C}^{\ast }\leqslant \unicode[STIX]{x1D70C}_{2}^{\ast }$
. Moreover, in appendix C, it is also shown that every
$\Vert \unicode[STIX]{x1D70C}^{\ast }\Vert _{L^{2m}}$
is bounded above by its initial data provided the advecting
$\boldsymbol{u}$
-field is regular. However, our interest lies more in
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }$
, both mathematically through its analogy with the vorticity field
$\unicode[STIX]{x1D714}$
and physically through its central role in irreversible mixing, although it is difficult to work with this quantity alone. To circumvent this problem, the variable
$\unicode[STIX]{x1D703}$
introduced in (1.9) and satisfying (1.10) is an easier variable with which to work. The idea is to consider both
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
and
$\unicode[STIX]{x1D74E}=\text{curl}\,\boldsymbol{u}$
in the higher (unscaled) norms
$L^{2m}({\mathcal{V}})$
(analogously, but slightly differently from those defined in (1.12) and (1.13)), defined for (
$1\leqslant m<\infty$
) by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}}=\left(\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2m}\,\text{d}V\right)^{1/2m}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D714}}=\left(\int _{{\mathcal{V}}}|\unicode[STIX]{x1D74E}|^{2m}\,\text{d}V\right)^{1/2m}, & \displaystyle\end{eqnarray}$$
where
$L_{0}$
is the non-dimensionalization length in the database, which rescales the volume of integration here to be
${\mathcal{V}}=[0,L/L_{0}]^{3}$
. The natural sequence of Hölder inequalities
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}}\leqslant \left(\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2(m+1)}\,\text{d}V\right)^{1/(2(m+1))}\left(\int _{{\mathcal{V}}}1\,\text{d}V\right)^{1/(2m(m+1))}=\left(\frac{L}{L_{0}}\right)^{3/(2m(m+1))}\unicode[STIX]{x1D6FA}_{m+1,\unicode[STIX]{x1D703}}\end{eqnarray}$$
has a multiplicative factor that is only unity when
$L=L_{0}$
. Remembering the definition of
$\unicode[STIX]{x1D6FC}_{m}$
from (1.13),
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{m}=\frac{2m}{4m-3},\end{eqnarray}$$
the exponent on
$L/L_{0}$
in (3.3) is related to
$\unicode[STIX]{x1D6FC}_{m}$
and
$\unicode[STIX]{x1D6FC}_{m+1}$
by
$$\begin{eqnarray}\frac{3}{2m(m+1)}=\frac{1}{\unicode[STIX]{x1D6FC}_{m+1}}-\frac{1}{\unicode[STIX]{x1D6FC}_{m}}.\end{eqnarray}$$
In turn, this leads us to define a natural dimensionless length
$$\begin{eqnarray}\ell _{m}=(L/L_{0})^{1/\unicode[STIX]{x1D6FC}_{m}},\end{eqnarray}$$
which turns (3.3) into
$\ell _{m}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}}\leqslant \ell _{m+1}\unicode[STIX]{x1D6FA}_{m+1,\unicode[STIX]{x1D703}}$
. Use of these length scales allows us to investigate the relationship between the various moments of the density gradients in the most natural fashion, to allow us to determine whether they exhibit nonlinear depletion (or not). The aim is to assume that there exists a solution of (1.10) in tandem with the vorticity field
$\unicode[STIX]{x1D74E}$
. Motivated by the depletion properties studied in Donzis et al. (Reference Donzis, Gibbon, Kerr, Gupta, Pandit and Vincenzi2013) and Gibbon et al. (Reference Gibbon, Donzis, Kerr, Gupta, Pandit and Vincenzi2014) for the Navier–Stokes equations, the following definitions are made:
$$\begin{eqnarray}D_{m,\unicode[STIX]{x1D703}}=(\ell _{m}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}})^{\unicode[STIX]{x1D6FC}_{m}},\quad D_{m,\unicode[STIX]{x1D714}}=(\ell _{m}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D714}})^{\unicode[STIX]{x1D6FC}_{m}}.\end{eqnarray}$$
In the JHT database the dimensionless domain size is
$2\unicode[STIX]{x03C0}$
, thus indicating that
$L/L_{0}=2\unicode[STIX]{x03C0}$
. The
$\unicode[STIX]{x1D6FC}_{m}$
-scaling in (3.7) has its origins in scaling properties of the three-dimensional Navier–Stokes equations (see Gibbon et al.
Reference Gibbon, Donzis, Kerr, Gupta, Pandit and Vincenzi2014). Then, as in (1.13) above, we consider
$$\begin{eqnarray}D_{m,\unicode[STIX]{x1D703}}(t)=D_{1,\unicode[STIX]{x1D703}}^{A_{m,\unicode[STIX]{x1D703}}(t)},\end{eqnarray}$$
where the multiplicative set of constants
$C_{m}$
have here been taken as unity. Following this, the JHT database shows that the relation between
$D_{m,\unicode[STIX]{x1D703}}$
and
$D_{1,\unicode[STIX]{x1D703}}$
takes the form of (3.8). The data are consistent with
$A_{m,\unicode[STIX]{x1D703}}(t)$
being expressed as
$$\begin{eqnarray}A_{m,\unicode[STIX]{x1D703}}(t)=\frac{\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)(m-1)+1}{4m-3},\quad 1\leqslant \unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}\leqslant 4.\end{eqnarray}$$
The
$D_{m,\unicode[STIX]{x1D703}}(t)$
are thus defined in terms of the set
$\{D_{1,\unicode[STIX]{x1D703}}(t),\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)\}$
. There is also an equivalent formula for
$D_{m,\unicode[STIX]{x1D714}}$
in terms of the set
$\{D_{1,\unicode[STIX]{x1D714}}(t),\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}(t)\}$
. There is a lower bound
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}\geqslant 1$
for the same reasons given in the introduction for the equivalent quantities
$D_{m}$
and
$D_{1}$
, as defined in (1.13). Numerically, the
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
can be calculated from the JHT database by considering
$\ln D_{m,\unicode[STIX]{x1D703}}/\ln D_{1,\unicode[STIX]{x1D703}}$
. It should be noted that the ordering observed in (3.3) does not necessarily hold for the
$D_{m,\unicode[STIX]{x1D703}}$
or the
$D_{m,\unicode[STIX]{x1D714}}$
because
$\unicode[STIX]{x1D6FC}_{m}$
decreases with
$m$
.
3.2 The evolution of
$D_{1,\unicode[STIX]{x1D703}}$
Formally, consider the time evolution of
$\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}\,\text{d}V$
using (1.10) and the relation for
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}$
in (1.10),
$$\begin{eqnarray}\displaystyle \frac{1}{2}\frac{\text{d}}{\text{d}t}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}\,\text{d}V & = & \displaystyle \int _{{\mathcal{V}}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }(Pe_{0}^{-1}\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D735}\boldsymbol{u})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\,\text{d}V+\frac{1}{2}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})\,\text{d}V\nonumber\\ \displaystyle & {\leqslant} & \displaystyle -\,Pe_{0}^{-1}\int _{{\mathcal{V}}}|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|^{2}\,\text{d}V+\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x1D735}\boldsymbol{u}|\,\text{d}V+\frac{1}{2}Pe_{0}^{-1}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|\,\text{d}V.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Using the fact that
$D_{1,\unicode[STIX]{x1D703}}=2\unicode[STIX]{x03C0}\,\unicode[STIX]{x1D6FA}_{1,\unicode[STIX]{x1D703}}^{2}$
together with the
$D_{m,\unicode[STIX]{x1D703}}-D_{1,\unicode[STIX]{x1D703}}$
relation in (3.8) and (3.9), appendix D shows how this can be manipulated into the differential inequality
$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}D_{1,\unicode[STIX]{x1D703}}\leqslant -\frac{D_{1,\unicode[STIX]{x1D703}}^{2}}{4\unicode[STIX]{x03C0}Pe_{0}\Vert \unicode[STIX]{x1D703}\Vert _{2}^{2}}+4\unicode[STIX]{x03C0}c_{m}Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)}+2\unicode[STIX]{x03C0}Pe_{0}D_{1,\unicode[STIX]{x1D714}},\end{eqnarray}$$
where the
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
appear in the exponent of the middle term on the right-hand side. Plots of
$\ell _{m}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}}(t)$
,
$D_{m,\unicode[STIX]{x1D703}}(t)$
and
$A_{m,\unicode[STIX]{x1D703}}$
are shown in figure 1, with plots of the corresponding
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
in figure 2(a); the relevant constants have not been estimated. It should be noted that the set
$\{\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)\}$
fan out with time, with no tendency to coincide. Nonlinear depletion occurs when
$A_{m,\unicode[STIX]{x1D703}}<1$
, which figure 1 shows to be the case. Because
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)\geqslant 1$
, the exponent
$1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)\geqslant 2$
. Given the values of
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
in figure 2(a), it is clear that, despite the boundedness of
$\Vert \unicode[STIX]{x1D703}\Vert _{2}^{2}$
, the negative
$D_{1,\unicode[STIX]{x1D703}}^{2}$
-term will always be smaller than the middle nonlinear term
$Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}}$
for all non-small values of
$D_{1,\unicode[STIX]{x1D703}}$
.
There is also the issue of the size of the additive
$2\unicode[STIX]{x03C0}Pe_{0}D_{1,\unicode[STIX]{x1D714}}$
-term in (3.11). To turn this term into a term depending on
$D_{1,\unicode[STIX]{x1D703}}$
alone, a relation is required between
$D_{1,\unicode[STIX]{x1D703}}$
and
$D_{1,\unicode[STIX]{x1D714}}$
, with the latter representing the fluid vorticity. Analytically, we have been unable to establish a relation between them, but the JHT database provides us with a relation with a further exponent which we call
$\unicode[STIX]{x1D6FD}(t)$
,
$$\begin{eqnarray}D_{1,\unicode[STIX]{x1D714}}=D_{1,\unicode[STIX]{x1D703}}^{\unicode[STIX]{x1D6FD}(t)},\end{eqnarray}$$
where the growth in the exponent
$\unicode[STIX]{x1D6FD}(t)$
is shown in figure 2(b). Figure 3 shows that the
$Pe_{0}D_{1,\unicode[STIX]{x1D703}}^{\unicode[STIX]{x1D6FD}(t)}$
-term (plotted with green squares) in (D 6) is dominant over the
$Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)}$
-term (plotted with red triangles), even when
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
is chosen to be the maximum across
$m$
at each particular time step. The plots of
$1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}$
and
$\unicode[STIX]{x1D6FD}(t)$
show that the values of these two quantities are both greater than two and thus cannot be controlled by the
$-D_{1,\unicode[STIX]{x1D703}}^{2}$
-term in (3.11). The conclusion is that these estimates do not point to any mechanism of saturation, and so the possibility of exponential growth in or even blow-up of
$D_{1,\unicode[STIX]{x1D703}}$
in a finite time cannot be discounted. It is important to stress that this is not inconsistent with the physical intuition that at late times, diffusion should smooth out all density gradients. This analysis only points to the possibility that at intermediate times it is possible that the density gradients may become extremely large. Furthermore, the ‘fanning out’ of the
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
for various
$m$
suggests a highly complex time-dependent fine structure in the density field, although figure 3 shows that this variation is not essential to our central conclusion concerning the time dependence of
$D_{1,\unicode[STIX]{x1D703}}$
.
Time variation of the
$Pe_{0}D_{1,\unicode[STIX]{x1D703}}^{\unicode[STIX]{x1D6FD}(t)}$
-term (plotted with green squares) and the
$Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)}$
-term (plotted with red triangles), where
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)$
is chosen to be the maximum value over
$m$
at each time step, with
$Pe_{0}=12\,500$
.
Finally, figure 4 shows the equivalent set of plots of the time variation of
$\ell _{m}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D714}}(t)$
,
$D_{m,\unicode[STIX]{x1D714}}(t)$
and
$A_{m,\unicode[STIX]{x1D714}}(t)$
, defined as
$$\begin{eqnarray}A_{m,\unicode[STIX]{x1D714}}(t)=\ln D_{m,\unicode[STIX]{x1D714}}/\ln D_{1,\unicode[STIX]{x1D714}}.\end{eqnarray}$$
In figure 5, we also show the time variation of the corresponding
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}(t)$
, calculated using the analogous relationship
$$\begin{eqnarray}A_{m,\unicode[STIX]{x1D714}}(t)=\frac{\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}(t)(m-1)+1}{(4m-3)}.\end{eqnarray}$$
It is apparent that the turbulent fluid part of the problem, which drives and dominates the system, has corresponding
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}(t)$
that are flat in time and sit in the range
$1<\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}<2$
. This is consistent with the behaviour found in three-dimensional Navier–Stokes flow described in Donzis et al. (Reference Donzis, Gibbon, Kerr, Gupta, Pandit and Vincenzi2013), Gibbon et al. (Reference Gibbon, Donzis, Kerr, Gupta, Pandit and Vincenzi2014) and Gibbon (Reference Gibbon2015). It should be noted that this contrasts strongly with the behaviour of the
$\unicode[STIX]{x1D703}$
-variable, where the
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}$
fan out and grow in time, as shown in figure 2, and suggests, in some as yet unexplored sense, that the density gradient field is even more complex than the vorticity field in this buoyancy-driven flow, in particular that the nonlinear depletion is appreciably weaker in the density gradient field.
Time variation of
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}(t)$
, calculated using the relation (3.14).
4 Conclusions
We have studied the growth of
$D_{1,\unicode[STIX]{x1D703}}=2\unicode[STIX]{x03C0}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}\,\text{d}V$
(and
$D_{1,\unicode[STIX]{x1D714}}=2\unicode[STIX]{x03C0}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D74E}|^{2}\,\text{d}V$
) based on the evolution of
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
with
$\unicode[STIX]{x1D703}=\ln (\unicode[STIX]{x1D70C}^{\ast }/\unicode[STIX]{x1D70C}_{0}^{\ast })$
as expressed in (1.10). The differential inequality for
$D_{1,\unicode[STIX]{x1D703}}$
shown in (3.11) contains the set of exponents
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}$
(
$m=2,\ldots ,8$
) together with
$\unicode[STIX]{x1D6FD}(t)$
, which is the exponent that connects
$D_{1,\unicode[STIX]{x1D703}}$
with
$D_{1,\unicode[STIX]{x1D714}}$
. These have been calculated directly from the JHTDB data which have been employed in this paper. The fan-like shape of the set
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}$
(
$m=2,\ldots ,8$
) and its growth in time are shown in figure 2(a) while the growth of
$\unicode[STIX]{x1D6FD}$
is shown in figure 2(b).
The extraction of statistics for the higher moments
$D_{m}$
for
$m\geqslant 2$
raises a similar question to that addressed in Donzis et al. (Reference Donzis, Gibbon, Kerr, Gupta, Pandit and Vincenzi2013) for the three-dimensional Navier–Stokes equations which was based on four data sets. In three of them (a
$2048\times 1024^{2}$
and two
$512^{3}$
calculations), statistics were able to be reliably extracted for up to
$m=9$
but not beyond. In the
$4096^{3}$
calculation, reliable extraction only up to
$m=6$
was achievable. Likewise, only reliable extraction up to
$m=9$
was achievable from the JHTDB. In the three-dimensional Navier–Stokes calculations, a strong degree of convergence was observed for
$m\geqslant 4$
, suggesting, but not proving, that this set might ultimately converge to a finite
$D_{\infty }$
. Given the fan-like separation of the
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}$
in figure 2, it is not clear one way or the other whether convergence ultimately occurs to a finite
$D_{\infty }$
.
In tandem with the numerical evidence in figure 2(a) consistent with strong growth in
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }$
, varying degrees of nonlinear depletion are also observed in the sense that
$A_{m,\unicode[STIX]{x1D703}}<1$
and
$A_{m,\unicode[STIX]{x1D714}}<1$
(as in figures 4
c and 5). Depletion in
$A_{m,\unicode[STIX]{x1D703}}$
reduces, as the growth of
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}$
to the value
$3.5$
in the final stages attests. Indeed, it should be noted that
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}=4$
would give a linear relation and be equivalent to a full estimate of the nonlinearity. Depletion in
$D_{1,\unicode[STIX]{x1D714}}$
is quite severe, as shown in figures 4(c) and 5, which is consistent with the same effect observed in Navier–Stokes flows. Despite this, the cross-effect of the turbulent fluid flow driving the growth of
$D_{1,\unicode[STIX]{x1D703}}$
through the exponent
$\unicode[STIX]{x1D6FD}(t)$
swamps the term
$D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}}$
in (3.11).
The fact that
$1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}\geqslant 2$
, with
$\unicode[STIX]{x1D6FD}(t)$
growing rapidly up to a value of 6, leaves open the possibility that
$D_{1,\unicode[STIX]{x1D703}}$
could blow up in a finite time. The rapid growth observed by Luo & Hou (Reference Luo and Hou2014a
,Reference Luo and Hou
b
) in their computations of the symmetric Euler equations with a boundary, with a clear indication of the formation of a finite time singularity, suggests that the problem in this paper should be analysed in this light, despite the differences in boundary conditions, the presence of diffusion and the limited grid size. While the blow-up question cannot be directly answered with the available data, the growth is sufficiently strong that the mixing could be driven down to near molecular scales where the validity of the model fails. Interestingly,
$\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D714}}$
simultaneously saturates, as in figure 5, showing that the turbulent fluid component is not amplified by the density growth. This hints that buoyancy-driven turbulence may well be more intense in some real physical and mathematical sense than constant-density turbulence, which may explain the observed extremely efficient mixing possible in such flows.
Acknowledgements
We acknowledge, with thanks, the staff of IPAM UCLA where this collaboration began in the Autumn of 2014 on the programme ‘Mathematics of Turbulence’. We would also like to thank C. Doering and D. Livescu for useful discussions. The research activity of C.P.C. is supported by EPSRC Programme grant EP/K034529/1 (‘Mathematical Underpinnings of Stratified Turbulence’). All of the numerical data used are from the JHTDB (Livescu et al. Reference Livescu, Canada, Kanov, Burns and Pulido2014), a publicly available DNS database. For more information, please see http://turbulence.pha.jhu.edu/. We also thank the referees for suggesting substantial improvements.
Appendix A. The equations for the composite density
Following Cook & Dimotakis (Reference Cook and Dimotakis2001) and Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), the composition density
$\unicode[STIX]{x1D70C}^{\ast }(\boldsymbol{x},t)$
of a mixture of two constant fluid densities
$\unicode[STIX]{x1D70C}_{1}^{\ast }$
and
$\unicode[STIX]{x1D70C}_{2}^{\ast }$
(
$\unicode[STIX]{x1D70C}_{2}^{\ast }>\unicode[STIX]{x1D70C}_{1}^{\ast }$
) is expressed by (1.2), where
$Y(\boldsymbol{x},t)=Y_{2}$
is the mass fraction of the heavier fluid. It is important to stress that the two fluids are assumed to be incompressible, yet we do not make the Boussinesq approximation and so the difference between the two densities is allowed to take arbitrary values. Under the transport of a (dimensionless) velocity field
$\boldsymbol{u}(\boldsymbol{x},t)$
,
$\unicode[STIX]{x1D70C}^{\ast }$
obeys the equation of conservation of mass,
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}^{\ast }+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }\boldsymbol{u})=0,\end{eqnarray}$$
and the species transport equation,
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}^{\ast }Y)+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }Y\boldsymbol{u}+\boldsymbol{j}_{F})=0,\end{eqnarray}$$
where the divergence of the flux
$\boldsymbol{j}_{F}$
represents Fickian diffusion, i.e.
$$\begin{eqnarray}\boldsymbol{j}_{F}=-Pe_{0}^{-1}\unicode[STIX]{x1D70C}^{\ast }\unicode[STIX]{x1D735}Y.\end{eqnarray}$$
The Péclet number has been defined in table 1. Given that the solution of (1.2) shows that
$\unicode[STIX]{x1D70C}^{\ast }Y$
is linear in
$\unicode[STIX]{x1D70C}^{\ast }$
such that
$$\begin{eqnarray}\unicode[STIX]{x1D70C}^{\ast }Y=a\unicode[STIX]{x1D70C}^{\ast }+b,\quad a=\frac{\unicode[STIX]{x1D70C}_{2}^{\ast }}{\unicode[STIX]{x1D70C}_{2}^{\ast }-\unicode[STIX]{x1D70C}_{1}^{\ast }},\quad b=-\frac{\unicode[STIX]{x1D70C}_{1}^{\ast }\unicode[STIX]{x1D70C}_{2}^{\ast }}{\unicode[STIX]{x1D70C}_{2}^{\ast }-\unicode[STIX]{x1D70C}_{1}^{\ast }},\end{eqnarray}$$
equation (A 2) simplifies to
$$\begin{eqnarray}b\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=Pe_{0}^{-1}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }\unicode[STIX]{x1D735}Y).\end{eqnarray}$$
Noting from (A 4) that
$\unicode[STIX]{x1D70C}^{\ast }\unicode[STIX]{x1D735}Y=-b\unicode[STIX]{x1D735}(\ln \unicode[STIX]{x1D70C}^{\ast })$
, the coefficient
$b$
cancels to make (A 5) and (A 1) into
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=-Pe_{0}^{-1}\unicode[STIX]{x1D6E5}(\ln \unicode[STIX]{x1D70C}^{\ast }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\unicode[STIX]{x1D70C}^{\ast }=Pe_{0}^{-1}\unicode[STIX]{x1D70C}^{\ast }\unicode[STIX]{x1D6E5}(\ln \unicode[STIX]{x1D70C}^{\ast }), & \displaystyle\end{eqnarray}$$
Appendix B. Details of the database
We have used the matlab analysis tools provided at the webpage http://turbulence.pha.jhu.edu/matlabanalysistools.aspx to fetch the data from the database. The specific dataset functions used for interpolating and calculating gradients are given in table 2. For details regarding the computational domain, boundary conditions, numerical method, etc., please refer to http://turbulence.pha.jhu.edu/docs/README-HBDT.pdf.
Function names in the database corresponding to the interpolation of various quantities used for data analysis in the paper.
Appendix C. Every
$\Vert \unicode[STIX]{x1D70C}^{\ast }\Vert _{L^{2m}({\mathcal{V}})}$
decays in time
To prove the boundedness of each
$\Vert \unicode[STIX]{x1D70C}^{\ast }\Vert _{L^{2m}({\mathcal{V}})}$
under a sufficiently regular advecting field
$\boldsymbol{u}$
, we write
$$\begin{eqnarray}\frac{1}{2m}\frac{\text{d}}{\text{d}t}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D70C}^{\ast }|^{2m}\,\text{d}V=-\int _{{\mathcal{V}}}\unicode[STIX]{x1D70C}^{\ast (2m-1)}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }\boldsymbol{u})\,\text{d}V\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D70C}^{\ast (2m-1)}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast }\boldsymbol{u})=\left(1-\frac{1}{2m}\right)\unicode[STIX]{x1D70C}^{\ast 2m}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}+\frac{1}{2m}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}^{\ast 2m}\boldsymbol{u}).\end{eqnarray}$$
Equation (C 1) then becomes
$$\begin{eqnarray}\frac{1}{2m}\frac{\text{d}}{\text{d}t}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D70C}^{\ast }|^{2m}\,\text{d}V=-\left(1-\frac{1}{2m}\right)\int _{{\mathcal{V}}}\unicode[STIX]{x1D70C}^{\ast 2m}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}\,\text{d}V,\end{eqnarray}$$
where the volume integral of the second term in (C 2) is zero through the divergence theorem. Using (2.3), (C 3) becomes
$$\begin{eqnarray}\displaystyle \frac{1}{2m}\frac{\text{d}}{\text{d}t}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D70C}^{\ast }|^{2m}\,\text{d}V & = & \displaystyle Pe_{0}^{-1}\left(1-\frac{1}{2m}\right)\int _{{\mathcal{V}}}\unicode[STIX]{x1D70C}^{\ast 2m}\unicode[STIX]{x1D6E5}(\ln \unicode[STIX]{x1D70C}^{\ast })\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle -Pe_{0}^{-1}(2m-1)\int _{{\mathcal{V}}}\unicode[STIX]{x1D70C}^{\ast 2(m-1)}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast }|^{2}\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle -Pe_{0}^{-1}\frac{(2m-1)}{m^{2}}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{\ast m}|^{2}\,\text{d}V\leqslant 0.\end{eqnarray}$$
Thus, every
$L^{2m}$
-norm decays in time.◻
Appendix D. The differential inequality for
$D_{1,\unicode[STIX]{x1D703}}$
Let us repeat (3.10) and proceed from there:
$$\begin{eqnarray}\frac{1}{2}\frac{\text{d}}{\text{d}t}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}\,\text{d}V\leqslant -Pe_{0}^{-1}\int _{{\mathcal{V}}}|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|^{2}\,\text{d}V+\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x1D735}\boldsymbol{u}|\,\text{d}V+\frac{1}{2}Pe_{0}^{-1}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|\,\text{d}V.\end{eqnarray}$$
Before going further, it should be noted that a standard result connects
$\unicode[STIX]{x1D735}\boldsymbol{u}$
,
$\unicode[STIX]{x1D74E}$
and
$\text{div}\,\boldsymbol{u}$
,
$$\begin{eqnarray}\displaystyle \int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\,\text{d}V & = & \displaystyle \int _{{\mathcal{V}}}|\unicode[STIX]{x1D74E}|^{2}\,\text{d}V+\int _{{\mathcal{V}}}|\text{div}\,\boldsymbol{u}|^{2}\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \int _{{\mathcal{V}}}|\unicode[STIX]{x1D74E}|^{2}\,\text{d}V+Pe_{0}^{-2}\int _{{\mathcal{V}}}|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|^{2}\,\text{d}V.\end{eqnarray}$$
For
$m\geqslant 2$
, and noting that
$(m-2)/(2(m-1))+m/(2(m-1))=1$
, consider the term
$$\begin{eqnarray}\displaystyle \int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x1D735}\boldsymbol{u}|\,\text{d}V & {\leqslant} & \displaystyle c_{1,m}\unicode[STIX]{x1D6FA}_{1,\unicode[STIX]{x1D703}}^{(m-2)/(m-1)}\unicode[STIX]{x1D6FA}_{m,\unicode[STIX]{x1D703}}^{m/(m-1)}[\unicode[STIX]{x1D6FA}_{1,\unicode[STIX]{x1D714}}+Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}]\nonumber\\ \displaystyle & {\leqslant} & \displaystyle c_{2,m}D_{1,\unicode[STIX]{x1D703}}^{(m-2)/(2(m-1))}D_{m,\unicode[STIX]{x1D703}}^{m/(\unicode[STIX]{x1D6FC}_{m}(m-1))}[(2\unicode[STIX]{x03C0})^{-1/2}D_{1,\unicode[STIX]{x1D714}}^{1/2}+Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}],\qquad\end{eqnarray}$$
where the factors of
$\ell _{m}$
have been absorbed into the dimensionless constants
$c_{1,m}$
and
$c_{2,m}$
. Now, we turn to the idea introduced in § 3.1 which connects
$D_{m,\unicode[STIX]{x1D703}}$
with
$D_{1,\unicode[STIX]{x1D703}}$
by using the formulae (3.8) and (3.9) for
$A_{m,\unicode[STIX]{x1D703}}$
and the set
$\{\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}(t)\}$
. Insertion of (3.8) into the right-hand side of (D 3) gives (factors of
$2\unicode[STIX]{x03C0}$
have also been absorbed into the constants)
$$\begin{eqnarray}\displaystyle \int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x1D735}\boldsymbol{u}|\,\text{d}V & {\leqslant} & \displaystyle c_{2,m}D_{1,\unicode[STIX]{x1D703}}^{(1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}})/2}[(2\unicode[STIX]{x03C0})^{-1/2}D_{1,\unicode[STIX]{x1D714}}^{1/2}+Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}]\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \frac{1}{2}Pe_{0}D_{1,\unicode[STIX]{x1D714}}+c_{3,m}Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}}+\frac{1}{2}Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}^{2},\end{eqnarray}$$
where the use of a Hölder inequality has split up the terms of the right-hand side of (D 3). The same idea is used on the last term in (D 1) with
$|\unicode[STIX]{x1D735}\boldsymbol{u}|$
replaced by
$|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|$
:
$$\begin{eqnarray}\displaystyle Pe_{0}^{-1}\int _{{\mathcal{V}}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}|^{2}|\unicode[STIX]{x0394}\unicode[STIX]{x1D703}|\,\text{d}V & {\leqslant} & \displaystyle (Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}^{2})^{1/2}(c_{3,m}Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}})^{1/2}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle \frac{1}{2}Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}^{2}+c_{4,m}Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}}.\end{eqnarray}$$
Altogether, (D 1) becomes
$$\begin{eqnarray}\frac{1}{4\unicode[STIX]{x03C0}}{\dot{D}}_{1,\unicode[STIX]{x1D703}}\leqslant -\frac{1}{4}Pe_{0}^{-1}\Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}^{2}+c_{5,m}Pe_{0}^{-1}D_{1,\unicode[STIX]{x1D703}}^{1+\unicode[STIX]{x1D706}_{m,\unicode[STIX]{x1D703}}}+\frac{1}{2}Pe_{0}D_{1,\unicode[STIX]{x1D714}}.\end{eqnarray}$$
A simple integration by parts shows that
$$\begin{eqnarray}\Vert \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\Vert _{2}^{2}\leqslant \Vert \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\Vert _{2}\Vert \unicode[STIX]{x1D703}\Vert _{2},\end{eqnarray}$$
which leads to (3.11) in which the constant
$c_{5,m}$
has been relabelled as
$c_{m}$
.
