1. Introduction
The Rayleigh–Taylor instability occurs for variable-density fluids in which a net acceleration applies from the lighter to the heavier fluid, or when heavy fluid is placed above light fluid in a gravitational field (Rayleigh Reference Rayleigh1882; Taylor Reference Taylor1950; Chandrasekhar Reference Chandrasekhar1961; Sharp Reference Sharp1984). The corresponding destabilization of the interface between the two fluids thus produces a mixing zone. These configurations are found in many applications, for instance in inertial confinement fusion (Lindl Reference Lindl1995), in astrophysics for supernovae (Cook & Cabot Reference Cook and Cabot2006), in geophysical flows with upwelling sites in the ocean (Cui & Street Reference Cui and Street2004) or with turbulence in the ionosphere (Molchanov Reference Molchanov2004). After the linear regime of initiation of the instability, nonlinearities gain importance, and the mixing zone becomes fully turbulent so that its width
$L$
eventually reaches an asymptotic self-similar state (Youngs Reference Youngs1984)
$$\begin{eqnarray}L=2{\it\alpha}\mathscr{A}gt^{2}\end{eqnarray}$$
where
$\mathscr{A}$
is the Atwood number,
$g>0$
the gravitational acceleration,
$t$
the time and
${\it\alpha}$
the parameter characterizing the mixing zone growth rate. Many experimental and numerical studies of Rayleigh–Taylor mixing zones have been performed, but no definite value of the parameter
${\it\alpha}$
was found (Dimonte et al.
Reference Dimonte, Youngs, Dimits, Weber, Mariak, Wunch, Garasi, Robinson, Andrews, Ramaprabhu, Calder, Fryxell, Biello, Dursi, Macneice, Olson, Ricker, Rosner, Timmes, Tufo, Young and Zingale2004; Youngs Reference Youngs2013). This shows that initial conditions have a strong influence on the late dynamics of the flow and make the prediction of the mixing zone evolution difficult.
Recent theoretical works propose linking the growth rate parameter to different turbulent statistics of the mixing layer: in Poujade & Peybernes (Reference Poujade and Peybernes2010),
${\it\alpha}$
is expressed in terms of the large-scale distribution of the turbulent structures, through the infrared slope of the turbulent kinetic energy spectrum. This approach to Rayleigh–Taylor turbulence can be interpreted as an extension of classical theories introduced to evaluate the decay of homogeneous isotropic turbulence (HIT) (Batchelor Reference Batchelor1949; Lesieur & Ossia Reference Lesieur and Ossia2000). Alternatively, since the infrared slope of the spectrum is not always easily measured, Gréa (Reference Gréa2013) proposes expressing the growth rate parameter as a function of the mixing intensity and its anisotropy. Unfortunately, it is difficult to check these predictions on the growth rate
${\it\alpha}$
in simulations or in experiments, since self-similar regimes are not always well established. These conclusions suggested a new approach focusing on the characterization of turbulent quantities in the mixing zone, with two salient features.
-
(a) First, in place of a mixing zone, we model unstably stratified homogeneous turbulence (USHT) (Thoroddsen, Van Atta & Yampolsky Reference Thoroddsen, Van Atta and Yampolsky1998; Griffond, Gréa & Soulard Reference Griffond, Gréa and Soulard2014; Soulard, Griffond & Gréa Reference Soulard, Griffond and Gréa2014), thus permitting a better handling of turbulence in a simplified framework that discards inhomogeneous effects. For that, one needs to assume that the ratio between the integral length scale of turbulence
$\ell$
and the characteristic scale
$L$
of the mixing region is small. This hypothesis is shown by Vladimirova & Chertkov (Reference Vladimirova and Chertkov2009) to be valid in well-developed Rayleigh–Taylor mixing zones (
$\ell /L\approx 0.2$
). -
(b) Second, the Rayleigh–Taylor mixing zone is therefore viewed as a stratified turbulent flow with a buoyancy frequency
$N=\sqrt{2\mathscr{A}g{\it\Gamma}}$
(
${\it\Gamma}$
being the mean concentration gradient of heavy fluid
${\sim}1/L$
). It characterizes the strength of the buoyancy force applied on the flow along the direction of acceleration. The resulting anisotropy has therefore to be included in models for turbulence statistics, so that classical models for isotropic homogeneous turbulence need to be modified. In the following, we will assume that the flow is axisymmetric with mirror symmetry about the plane orthogonal to the vertical axis, meaning that its statistics do not depend on either horizontal direction in physical space or in spectral space. The axis of symmetry is along the gravity vector.
Hence, USHT can be understood as an idealized configuration to study buoyancy-driven turbulence, and in particular Rayleigh–Taylor mixing. Also, it can be interpreted as an extension of previous studies on buoyancy-induced turbulence, for instance by Batchelor, Canuto & Chasnov (Reference Batchelor, Canuto and Chasnov1992), who first introduced this approach, and by Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007) and Chung & Pullin (Reference Chung and Pullin2010). In these works, the adopted point of view is slightly different: a homogeneous flow is considered, as in USHT, but they model statistically stationary or decaying variable-density turbulence and do not investigate the development of the instability.
Let us now describe briefly some features of USHT. Figure 1 shows visualizations of the fluctuating scalar field from USHT simulations corresponding to run 5 (detailed later in § 3.1). The work of unstable buoyancy forces can be identified from anisotropic structures that are elongated along the vertical direction. From the plots at different times, one observes that these structures quickly develop from the initial isotropic state. The dynamics of USHT thus corresponds to an increase of turbulent kinetic energy and Reynolds number due to the input of energy by mean stratification, and to a growth of integral length scale
$\ell$
characterizing the large-scale eddies. However, beyond these phenomenological aspects, the question of how to predict the dynamics of USHT is still open, and the related unsolved question is: what is the growth rate of turbulent quantities in unstably stratified flows? We therefore investigate whether a self-similar regime is reached, and how long after the beginning of the development of the turbulent flow.
Figure 1. Visualizations of the fluctuating buoyancy field extracted from DNS in USHT at different times (
$t^{\ast }=0,2.5,5$
) and for vertical and horizontal planes corresponding to run 5 in table 1.
In the case of constant stratification frequency
$N$
, the expected self-similar regime is characterized by an exponential kinetic energy evolution,
${\sim}\text{e}^{{\it\beta}Nt}$
. One can estimate
${\it\beta}$
using arguments proposed by Soulard et al. (Reference Soulard, Griffond and Gréa2014) for time-dependent
$N$
, namely that the dynamics of large scales is self-similar, determined by the buoyancy force, and that nonlinear effects can be discarded. In that case,
${\it\beta}$
for USHT – equivalent to
${\it\alpha}$
in Rayleigh–Taylor – can be related to the slope
$s$
of infrared spectra (see figure 2) as (see derivation in appendix A)
$$\begin{eqnarray}{\it\beta}=\frac{4}{s+3}.\end{eqnarray}$$
In order to describe further the self-similar characteristics of buoyancy-induced turbulence, let us sketch what the energy spectrum should be like. At smaller scales and at high Reynolds number, several authors, for instance Chertkov (Reference Chertkov2003) and Soulard (Reference Soulard2012), argue that a
$k^{-5/3}$
inertial range appears in the energy spectrum
$E(k,t)$
, consistent with Kolmogorov theory. On the contrary, rather than a scaling invariant based on dissipation and leading to the
$k^{-5/3}$
scaling, other invariants are proposed by Zhou (Reference Zhou2001) and Abarzhi (Reference Abarzhi2010) yielding different spectral scalings. Moreover, the vertical component of the anisotropic spectra should scale as
$k^{-3}$
according to the perturbative approach based on linear response theory (Soulard & Griffond Reference Soulard and Griffond2012). This suggests that the energy spectrum should straighten when approaching the integral scale. However, it is difficult to predict the precise extent of this zone, bounded by the integral and the Ozmidov scales. For large scales, it can be shown from quasi-normal approximations that, at least for
$s>4$
, backscatter plays a decisive role in the possible variation of the infrared slope
$s$
in time. Therefore, the controversial question of permanence of big eddies in HIT (Batchelor Reference Batchelor1953) has its counterpart in buoyancy-induced turbulence. Taken together, these elements yield a self-similar spectrum as in figure 2, with typical wavenumbers
$k_{\ell }$
,
$k_{O}$
and
$k_{{\it\eta}}$
associated with respectively the integral, the Ozmidov and the Kolmogorov scales. Such spectra for USHT have never been clearly observed as they require high Reynolds number values. In addition, most of the above arguments are based on dimensional analysis that does not account for the significant flow anisotropy.
Figure 2. Schematic representation of a self-similar kinetic energy spectrum
$E(k,t)$
, in USHT with
$N=\text{const}$
. The different regimes are shown that would correspond to separate spectral ranges if enough scale separation is achieved. The infrared range corresponds to the largest structures with a scaling
$k^{s}$
.
In this paper, our goal is to investigate USHT phenomenology and predict the ultimate value asymptotically reached by the growth parameter, and never really measured, as already mentioned, due to experimental and numerical limitations. For instance, in the experiment on USHT by Thoroddsen et al. (Reference Thoroddsen, Van Atta and Yampolsky1998), overturning due to Rayleigh–Bénard-like instability prevented the observation of the long-term evolution of the unstably stratified flow. We shall still use direct numerical simulations (DNS) as a reference for the model developed below, but the resolved spectral range in DNS is limited by computational resources, and so is the Reynolds number. This is due to the trade-off between a sufficient extent of the infrared spectrum, and sufficiently well-resolved small scales. Since the available spectral range evolves linearly with the number of points in each direction in pseudo-spectral DNS, no clear regime with scale separation is reached even in the largest existing DNS (for instance by Cook & Cabot Reference Cook and Cabot2006).
The imposed truncation of the infrared part of the kinetic energy spectrum, or, in other words, the numerical confinement of the largest turbulent scales, prevents a reliable study of the dependence of
${\it\beta}$
on the large-scale properties of the flow, as well as a proper identification of scaling laws in the self-similar spectra. We propose in this work to address these two limitations by introducing a new model for USHT that is able to reach larger Reynolds number values than in DNS, and to represent a wider spectral range. Since we focus on the statistics of USHT, we develop a new statistical two-point model for the dynamical equations, based on the eddy damped quasi-normal Markovian closure (EDQNM). This model was first derived by Orszag (Reference Orszag1969) for isotropic homogeneous turbulence (see also Pouquet et al.
Reference Pouquet, Lesieur, André and Basdevant1975; André & Lesieur Reference André and Lesieur1977). New derivations extended the model to decaying axisymmetric turbulence (Herring Reference Herring1974) and to anisotropic turbulence subjected to the Coriolis force (Cambon, Mansour & Godeferd Reference Cambon, Mansour and Godeferd1997), to the Lorentz force (Favier et al.
Reference Favier, Godeferd, Cambon, Delache and Bos2011), to stably stratified flows with uniform mean density gradient (Godeferd & Cambon Reference Godeferd and Cambon1994), and even to variable-viscosity flows (Gréa, Griffond & Burlot Reference Gréa, Griffond and Burlot2014). In the latter cases, the EDQNM model was compared to experimental and numerical results with excellent agreement.
The plan of the article is as follows. The dynamical equations and some elements of the closure are presented in § 2, as well as its numerical implementation, a short presentation of our DNS algorithm, and how solutions of linearized equations are computed. We compare in § 3 the EDQNM predictions to DNS results with the same parameters and initial conditions. From this first set of simulations, we conclude that the model can better represent the nonlinear dynamics by a simple modification of its core parameter, the eddy damping, to account for the intensity of stratification. Thus, the amended version of the model is proposed in § 4. Finally, predictions of the EDQNM model for a turbulent regime at high Reynolds number are presented in § 5.
2. Governing equations for USHT
In this section, we introduce the dynamical equations for USHT, we present the main features of the two-point EDQNM closure for axisymmetric turbulence, and describe our different predictive approaches: the EDQNM closure, DNS and linearized solutions.
2.1. Basic equations in physical space
The system considered is the Navier–Stokes equations for an incompressible flow in the Boussinesq approximation. The flow is assumed to be statistically homogeneous, and the background concentration gradient for the heavy fluid, of amplitude
${\it\Gamma}=\partial \langle C\rangle /\partial x_{3}$
, is applied along the
$x_{3}$
(vertical) direction. Here,
$\langle \,\rangle$
denotes ensemble averaging. The gravity acceleration and the mean stratification gradient are aligned but opposite for an unstable configuration. Statistical axisymmetry is assumed about this
$x_{3}$
axis, and the flow distribution is isotropic in the
$(x_{1},x_{2})$
plane. In addition, the
$x_{3}$
dependence of the mean flow (appearing through
${\it\Gamma}$
) can be neglected since we assume that the largest scales of turbulence are small with respect to the size of the mixing zone. Due to the Boussinesq approximation and symmetries, the mean velocity is uniform and can be taken as zero,
$\langle u_{i}\rangle =0$
. The velocity field
$\boldsymbol{u}(\boldsymbol{x},t)$
and the fluctuating buoyancy rescaled as a velocity
${\it\vartheta}(\boldsymbol{x},t)$
(with
$\langle {\it\vartheta}\rangle =0$
) evolve dynamically following
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}}{\partial x_{i}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{i}}{\partial t}+u_{j}\frac{\partial u_{i}}{\partial x_{j}}=-\frac{1}{{\it\rho}_{ref}}\frac{\partial p}{\partial x_{i}}+{\it\nu}\frac{\partial ^{2}u_{i}}{\partial x_{j}\partial x_{j}}+N{\it\vartheta}{\it\delta}_{i3}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\vartheta}}{\partial t}+u_{j}\frac{\partial {\it\vartheta}}{\partial x_{j}}=\mathscr{D}\frac{\partial ^{2}{\it\vartheta}}{\partial x_{j}\partial x_{j}}+Nu_{3}, & \displaystyle\end{eqnarray}$$
${\it\rho}_{light}$
of light and
${\it\rho}_{heavy}$
of heavy fluid as
${\it\rho}_{ref}=({\it\rho}_{light}+{\it\rho}_{heavy})/2$
. The fluctuating pressure is denoted
$p$
(with
$\langle p\rangle =0$
), the kinematic viscosity
${\it\nu}$
, the molecular diffusivity
$\mathscr{D}$
. In addition,
$N=\sqrt{2\mathscr{A}g{\it\Gamma}}$
, characterizing the instability, is an inverse timescale analogous to the Brunt–Väisälä frequency for stably stratified turbulence, with
$\mathscr{A}=({\it\rho}_{heavy}-{\it\rho}_{light})/({\it\rho}_{heavy}+{\it\rho}_{light})$
the Atwood number that represents the relative density contrast between the two mixed fluids. The relation between the buoyancy velocity and fluctuations of concentration
$c$
is therefore
${\it\vartheta}=2\mathscr{A}gc/N$
.
2.2. Equations in Fourier space and axisymmetric representation
Anticipating that we are going to propose a model for spectra of two-point correlations, and from the homogeneity assumption, we start by transforming the dynamical equations to Fourier space before defining the spectra. This method proved to be simpler than Fourier-transforming equations for two-point correlation functions (Craya Reference Craya1958).
We denote
$\hat{u} _{i}\!\left(\boldsymbol{k},t\right)$
and
$\hat{{\it\vartheta}}\!\left(\boldsymbol{k},t\right)$
the Fourier coefficients respectively of
$u_{i}\!\left(\boldsymbol{x},t\right)$
and
${\it\vartheta}\!\left(\boldsymbol{x},t\right)$
, associated to the wavevector
$\boldsymbol{k}$
. Fourier-transforming (2.1)–(2.3) provides the following equations for the spectral coefficients:
$$\begin{eqnarray}\displaystyle & k_{i}\hat{u} _{i}(\boldsymbol{k})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\partial _{t}+{\it\nu}k^{2}\right)\hat{u} _{i}(\boldsymbol{k})=-\text{i}k_{j}P_{il}\!\left(\boldsymbol{k}\right)\int _{\mathbb{R}^{3}}\hat{u} _{j}\!\left(\boldsymbol{p}\right)\hat{u} _{l}\!\left(\boldsymbol{k}-\boldsymbol{p}\right)\!\text{d}\boldsymbol{p}+NP_{i3}\!\left(\boldsymbol{k}\right)\hat{{\it\vartheta}}(\boldsymbol{k}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \!\left(\partial _{t}+\mathscr{D}k^{2}\right)\hat{{\it\vartheta}}(\boldsymbol{k})=-\text{i}k_{j}\int _{\mathbb{R}^{3}}\hat{u} _{j}(\boldsymbol{p})\hat{{\it\vartheta}}(\boldsymbol{k}-\boldsymbol{p})\text{d}\boldsymbol{p}+N\hat{u} _{3}(\boldsymbol{k}), & \displaystyle\end{eqnarray}$$
$\text{i}^{2}=-1$
. The projector
$P_{ij}\!\left(\boldsymbol{k}\right)={\it\delta}_{ij}-k_{i}k_{j}/k^{2}$
with
$k^{2}=k_{i}k_{i}$
appears after taking the divergence of (2.2). It leads to a Poisson equation for pressure which can be simplified thanks to the incompressibility hypothesis stated in (2.1) with a differential operator that becomes a geometrical one in the spectral equation (2.4). In the following, we shall also use twice the Kraichnan projector
$P_{ilm}\!\left(\boldsymbol{k}\right)=k_{l}P_{im}\!\left(\boldsymbol{k}\right)+k_{m}P_{il}\!\left(\boldsymbol{k}\right)$
.
The last terms in (2.5) and (2.6) represent a linear operator that couples velocity and scalar fields. The eigenvalues of this linear operator are: one neutral eigenvalue
${\it\sigma}=0$
associated with a toroidal eigenmode; two non-zero eigenvalues
${\it\sigma}=\pm Nk_{3}/(k_{1}^{2}+k_{2}^{2})^{1/2}=\pm N\sin {\it\theta}$
that represent the time-evolving modes. When projected onto the frame of reference defined by these eigenmodes, the linearized inviscid equations (2.5) and (2.6) for the Fourier components become diagonal, and so do the linearized equations for higher-order moments of the Fourier components. We shall use this in § 2.2.3.
The
$n$
th-order statistics of the velocity (or scalar) field can be computed in a general way through correlations such as
$\langle u_{i_{1}}(\boldsymbol{x}_{1})u_{i_{2}}(\boldsymbol{x}_{2})\cdots u_{i_{n}}(\boldsymbol{x}_{n})\rangle$
. Velocity/scalar cross-correlations can also be defined in the same way. The spectra of these statistics may be computed from the Fourier components themselves, such as, for the two-point second-order moments
$$\begin{eqnarray}\displaystyle & \langle \hat{u} _{i}(\boldsymbol{k})\hat{u} _{j}(\boldsymbol{k}^{\prime })\rangle =\mathscr{R}_{ij}(\boldsymbol{k}){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle \hat{{\it\vartheta}}(\boldsymbol{k})\hat{u} _{i}(\boldsymbol{k}^{\prime })\rangle =\mathscr{F}_{i}(\boldsymbol{k}){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle \hat{{\it\vartheta}}(\boldsymbol{k})\hat{{\it\vartheta}}(\boldsymbol{k}^{\prime })\rangle =\mathscr{B}(\boldsymbol{k}){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime }), & \displaystyle\end{eqnarray}$$
$\mathscr{R}_{ij}$
is the Fourier transform of
$R_{ij}(\boldsymbol{r})=\langle u_{i}(\boldsymbol{x})u_{j}(\boldsymbol{x}+\boldsymbol{r})\rangle$
with respect to the separation vector
$\boldsymbol{r}$
(homogeneity implies that
$R_{ij}$
does not depend on
$\boldsymbol{x}$
), and
${\it\delta}$
is the three-dimensional Dirac function. Likewise, one obtains the two-point buoyancy correlation spectrum
$\mathscr{B}$
as the Fourier transform of
$\langle {\it\vartheta}(\boldsymbol{x}){\it\vartheta}(\boldsymbol{x}+\boldsymbol{r})\rangle$
and the two-point velocity/scalar cross-correlation spectrum
$\mathscr{F}_{i}$
corresponding to
$\langle u_{i}(\boldsymbol{x}){\it\vartheta}(\boldsymbol{x}+\boldsymbol{r})\rangle$
. Note that the dependence of these spectra on
$\boldsymbol{k}$
reflects the two-point nature of the correlation. Similarly, third-order three-point correlation spectra can be defined as 
$$\begin{eqnarray}\displaystyle & \langle \hat{u} _{i}(\boldsymbol{k})\hat{u} _{j}(\boldsymbol{k}^{\prime })\hat{u} _{n}(\boldsymbol{k}^{\prime \prime })\rangle =S_{ijn}(\boldsymbol{k},\boldsymbol{k}^{\prime }){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime }+\boldsymbol{k}^{\prime \prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle \hat{{\it\vartheta}}(\boldsymbol{k})\hat{u} _{i}(\boldsymbol{k}^{\prime })\hat{u} _{j}(\boldsymbol{k}^{\prime \prime })\rangle =S_{ij}(\boldsymbol{k},\boldsymbol{k}^{\prime }){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime }+\boldsymbol{k}^{\prime \prime }), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \langle \hat{{\it\vartheta}}(\boldsymbol{k})\hat{{\it\vartheta}}(\boldsymbol{k}^{\prime })\hat{u} _{i}(\boldsymbol{k}^{\prime \prime })\rangle =S_{i}(\boldsymbol{k},\boldsymbol{k}^{\prime }){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime }+\boldsymbol{k}^{\prime \prime }). & \displaystyle\end{eqnarray}$$
$S_{ijn}$
,
$S_{ij}$
and
$S_{i}$
depend on both
$\boldsymbol{k}$
and
$\boldsymbol{k}^{\prime }$
as an indication that the corresponding correlations are taken at three physical points. In spectral space, this is reflected by the triad (triangle) of wavevectors
$\boldsymbol{k}+\boldsymbol{k}^{\prime }+\boldsymbol{k}^{\prime \prime }=\mathbf{0}$
.
2.2.1. Equations for second-order moments
Here we shall consider flows with unit Schmidt number
$\mathit{Sc}={\it\nu}/\mathscr{D}$
so that
$\mathscr{D}={\it\nu}$
. This assumption simplifies the diagonalization of the viscous–diffusive operator in the following equations, although it is not necessary, e.g. if one wishes to retain differential diffusive effects.
The equations for the two-point velocity correlation spectra are obtained by adding the dynamical equations for
$\hat{u} _{i}(\boldsymbol{k})$
and for
$\hat{u} _{j}(\boldsymbol{k}^{\prime })$
respectively weighted by
$\hat{u} _{j}(\boldsymbol{k}^{\prime })$
and
$\hat{u} _{i}(\boldsymbol{k})$
, and taking the ensemble average (similarly for scalar and cross-correlation spectra). This yields:
$$\begin{eqnarray}\left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right)\mathscr{R}_{ij}(\boldsymbol{k},t)=NP_{i3}\!\left(\boldsymbol{k}\right)\mathscr{F}_{j}(\boldsymbol{k})+NP_{j3}\!\left(\boldsymbol{k}\right)\mathscr{F}_{i}(-\boldsymbol{k})+T_{ij}^{\mathscr{R}}(\boldsymbol{k},t)\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle & \displaystyle T_{ij}^{\mathscr{R}}(\boldsymbol{k},t)=-\frac{\text{i}}{2}P_{inm}(\boldsymbol{k})\int _{\mathbb{R}^{3}}S_{jnm}(-\boldsymbol{k},\boldsymbol{p})\text{d}\boldsymbol{p}+\frac{\text{i}}{2}P_{jmn}(\boldsymbol{k})\int _{\mathbb{R}^{3}}S_{inm}(\boldsymbol{k},\boldsymbol{p})\text{d}\boldsymbol{p}; & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right)\mathscr{F}_{i}(\boldsymbol{k})=N\mathscr{R}_{3i}(\boldsymbol{k})+NP_{j3}\!\left(\boldsymbol{k}\right)\mathscr{B}(\boldsymbol{k})+T_{i}^{\mathscr{F}}(\boldsymbol{k},t) & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle T_{i}^{\mathscr{F}}(\boldsymbol{k},t)=\frac{\text{i}}{2}P_{inm}(\boldsymbol{k})\int _{\mathbb{R}^{3}}S_{nm}(\boldsymbol{k},\boldsymbol{p})\text{d}\boldsymbol{p}-\text{i}k_{n}\int _{\mathbb{R}^{3}}S_{in}(\boldsymbol{p},-\boldsymbol{k})\text{d}\boldsymbol{p}; & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right)\mathscr{B}(\boldsymbol{k},t)=N\!\left(\mathscr{F}_{3}(\boldsymbol{k})+\mathscr{F}_{3}(-\boldsymbol{k})\right)+T^{\mathscr{B}}(\boldsymbol{k},t) & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}T^{\mathscr{B}}(\boldsymbol{k},t)=-\text{i}k_{n}\int _{\mathbb{R}^{3}}S_{n}(-\boldsymbol{k},\boldsymbol{p})\text{d}\boldsymbol{p}+\text{i}k_{n}\int _{\mathbb{R}^{3}}S_{n}(\boldsymbol{k},\boldsymbol{p})\text{d}\boldsymbol{p}.\end{eqnarray}$$
Nonlinear transfer terms
$T_{ij}^{\mathscr{R}}$
,
$T_{i}^{\mathscr{F}}$
and
$T^{\mathscr{B}}$
in (2.13)–(2.17) represent triadic exchanges computed from third-order correlations. The additional right-hand-side terms are linear couplings between velocity-, scalar- and cross-spectra. Equations (2.13)–(2.17) are thus not closed and do not provide a means of computing two-point correlation spectra unless one is able to estimate third-order correlations from second-order ones. We shall do this in § 2.2.3, after introducing a convenient reference frame accounting for the statistical axisymmetry of the flow in § 2.2.2.
2.2.2. Polar-spherical decomposition for axisymmetric spectra
We shall now transpose the spectral description to a polar-spherical frame of reference, thus taking advantage of the solenoidality of the velocity field and axisymmetric property of the spectra about the axis of gravity which is along the unit vector
$\boldsymbol{x}_{3}$
, and hence providing a self-consistent description with respect to the symmetries of the flow. A second advantage is a reduction of the number of spectral functions in the polar-spherical description with respect to the Cartesian frame of reference:
$\hat{\boldsymbol{u}}$
has a priori three components in the latter frame, but, from (2.4), only two components are required to describe it in a reference frame of the polar-spherical coordinates in which one of the unit vectors is along
$\boldsymbol{k}$
. This frame of reference has proved to be physically meaningful, in addition to being mathematically useful, for various axisymmetric contexts such as stably stratified flows (Godeferd & Cambon Reference Godeferd and Cambon1994), rotating and stratified flows (Liechtenstein, Godeferd & Cambon Reference Liechtenstein, Godeferd and Cambon2005), and flows of conducting fluids with an imposed external magnetic field (Favier et al.
Reference Favier, Godeferd, Cambon, Delache and Bos2011). It is also called the Craya–Herring frame (Craya Reference Craya1958; Herring Reference Herring1974) (see also Sagaut & Cambon Reference Sagaut and Cambon2008 for its close relationship to the poloidal/toroidal decomposition in physical space). This frame is defined by the following unit vectors:
$$\begin{eqnarray}\boldsymbol{e}^{(1)}\!\!\left(\boldsymbol{k}\right)=\frac{\boldsymbol{k}\times \boldsymbol{x}_{3}}{\Vert \boldsymbol{k}\times \boldsymbol{x}_{3}\Vert },\quad \boldsymbol{e}^{(2)}\!\!\left(\boldsymbol{k}\right)=\frac{\boldsymbol{k}\times \boldsymbol{e}^{(1)}\!\!\left(\boldsymbol{k}\right)}{\Vert \boldsymbol{k}\times \boldsymbol{e}^{(1)}\!\!\left(\boldsymbol{k}\right)\!\Vert },\quad \boldsymbol{e}^{(3)}\!\!\left(\boldsymbol{k}\right)=\frac{\boldsymbol{k}}{\Vert \boldsymbol{k}\Vert },\end{eqnarray}$$
onto which the velocity field is projected, and hence has two components
$$\begin{eqnarray}\hat{u} _{tor}(\boldsymbol{k})=\hat{u} _{i}(\boldsymbol{k})e_{i}^{(1)}\!\!\left(\boldsymbol{k}\right)\!,\quad \hat{u} _{pol}(\boldsymbol{k})=\hat{u} _{i}(\boldsymbol{k})e_{i}^{(2)}\!\!\left(\boldsymbol{k}\right)\!,\end{eqnarray}$$
the third-one, along
$\boldsymbol{e}^{(3)}$
, being zero from incompressibility. Therefore, the algebra involves a single global velocity–buoyancy vector (
$\hat{u} _{tor}$
,
$\hat{u} _{pol}$
,
$\hat{{\it\vartheta}}$
) (see Godeferd & Cambon Reference Godeferd and Cambon1994).
We can then rewrite the second-order moments of velocity, buoyancy and the velocity–buoyancy covariance projected onto the polar-spherical unit vectors as
$$\begin{eqnarray}\displaystyle & \mathscr{R}_{ij}(\boldsymbol{k})={\it\Phi}_{1}(\boldsymbol{k})e_{i}^{(1)}\!\!\left(\boldsymbol{k}\right)e_{j}^{(1)}\!\!\left(\boldsymbol{k}\right)+{\it\Phi}_{2}(\boldsymbol{k})e_{i}^{(2)}\!\!\left(\boldsymbol{k}\right)e_{j}^{(2)}\!\!\left(\boldsymbol{k}\right)\!, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \mathscr{F}_{i}(\boldsymbol{k})={\it\Psi}_{r}(\boldsymbol{k})e_{i}^{(2)}\!\!\left(\boldsymbol{k}\right)\!, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \mathscr{B}(\boldsymbol{k})={\it\Phi}_{3}(\boldsymbol{k}), & \displaystyle\end{eqnarray}$$
${\it\Phi}_{1}(\boldsymbol{k}){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime })=\langle \hat{u} _{tor}(\boldsymbol{k})\hat{u} _{tor}(\boldsymbol{k}^{\prime })\rangle$
, poloidal energy spectrum
${\it\Phi}_{2}(\boldsymbol{k}){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime })=\langle \hat{u} _{pol}(\boldsymbol{k})\hat{u} _{pol}(\boldsymbol{k}^{\prime })\rangle$
, and the velocity/buoyancy co-spectrum
${\it\Psi}(\boldsymbol{k}){\it\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime })=\langle \hat{u} _{pol}(\boldsymbol{k})\hat{{\it\vartheta}}(\boldsymbol{k}^{\prime })\rangle$
with real part
${\it\Psi}_{r}(\boldsymbol{k})$
. The modal spectra
${\it\Phi}_{1}$
and
${\it\Phi}_{2}$
add up to twice the kinetic energy spectrum
$E$
. Also, as we assume mirror symmetry in the flow, helicity is zero and no additional scalar functions are needed to represent
$\mathscr{R}_{ij}$
,
$\mathscr{F}_{i}$
,
$\mathscr{B}$
.
One can then recast the three statistical equations (2.13), (2.14) and (2.15) into dynamical equations for the
${\it\Phi}$
and
${\it\Psi}$
spectra, and express the axisymmetric dependence of the spectra in terms of
$k=|\boldsymbol{k}|$
and the polar angle
${\it\theta}$
between
$\boldsymbol{k}$
and
$\boldsymbol{x}_{3}$
. This yields
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right){\it\Phi}_{1}(k,{\it\theta})=T^{{\it\Phi}_{1}}(k,{\it\theta}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right){\it\Phi}_{2}(k,{\it\theta})=T^{{\it\Phi}_{2}}(k,{\it\theta})+2N\sin ({\it\theta}){\it\Psi}_{r}(k,{\it\theta}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right){\it\Phi}_{3}(k,{\it\theta})=T^{{\it\Phi}_{3}}(k,{\it\theta})+2N\sin ({\it\theta}){\it\Psi}_{r}(k,{\it\theta}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial }{\partial t}+2{\it\nu}k^{2}\right){\it\Psi}_{r}(k,{\it\theta})=T^{{\it\Psi}_{r}}(k,{\it\theta})+N\sin ({\it\theta})\left({\it\Phi}_{2}(k,{\it\theta})+{\it\Phi}_{3}(k,{\it\theta})\right), & \displaystyle\end{eqnarray}$$
$\sin {\it\theta}$
arises from inner products between unit vectors of the Craya–Herring frame. Equations (2.24)–(2.27) are the exact evolution equations for second-order moments, in which the transfer terms on the right-hand sides contain all the triadic interactions (triple correlations of velocity and buoyancy components). The linear terms on the right-hand sides of (2.25) and (2.26) correspond to the direct coupling between velocity and buoyancy mediated by a buoyancy flux which determines the dynamics of
${\it\Phi}_{2}+{\it\Phi}_{3}$
in the absence of nonlinearity.
From the energy-density modal spectra whose dynamics is described by (2.24)–(2.27), other statistics can be obtained. For instance, spectra integrated over spherical shells are computed by integration over the angle
${\it\theta}$
, with the azimuthal dependence analytically integrated into a
$2{\rm\pi}k^{2}$
weight from axial symmetry: the turbulent kinetic energy spectrum is therefore
$E(k)={\rm\pi}k^{2}\int _{0}^{{\rm\pi}}({\it\Phi}_{1}(k,{\it\theta})+{\it\Phi}_{2}(k,{\it\theta}))\sin {\it\theta}\text{d}{\it\theta}$
; likewise the vertical flux spectrum
$F(k)=2{\rm\pi}k^{2}\int _{0}^{{\rm\pi}}{\it\Psi}_{r}(k,{\it\theta})\sin ^{2}{\it\theta}\text{d}{\it\theta}$
, and buoyancy spectrum
$B(k)=2{\rm\pi}k^{2}\int _{0}^{{\rm\pi}}{\it\Phi}_{3}(k,{\it\theta})\sin {\it\theta}\text{d}{\it\theta}$
. And one-point statistics are obtained by implementing integrals over
$k$
. For instance, the vertical mass flux is
$\langle u_{3}{\it\vartheta}\rangle =\int _{0}^{\infty }F(k)\text{d}k$
.
2.2.3. Two-point statistics equations closure
We briefly present here the technique used to model the statistical equations (2.24)–(2.27), using a closure at the level of third-order correlation spectra. The method, described by Orszag (Reference Orszag1973) for isotropic turbulence, provides a closed dynamical equation for a spherically averaged kinetic energy spectrum
$E(k)$
. It is more complicated in axisymmetric turbulence, which requires the inclusion of a spectral dependence on the polar angle
${\it\theta}$
in addition to the dependence on the wavenumber
$k$
. The derivation is only described here in the diagonalizing frame of the linear operator for brevity, but the closed equations in the axisymmetric frame are obtained by projecting their expressions on the Craya–Herring frame. (The reader is referred to e.g. Sagaut & Cambon Reference Sagaut and Cambon2008 for details.)
One needs to estimate third-order moments
$S$
defined in (2.10)–(2.12), upon which the transfer terms in (2.24)–(2.27) are based. In the reference frame of the eigenmodes of the linear stratification operator, mentioned previously, the dynamical equations for these third-order moments are
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\frac{\partial }{\partial t}+{\it\nu}\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)-N\left(s\sin {\it\theta}+s^{\prime }\sin {\it\theta}^{\prime }+s^{\prime \prime }\sin {\it\theta}^{\prime \prime }\right)\right)S_{ss^{\prime }s^{\prime \prime }}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)\nonumber\\ \displaystyle & & \displaystyle \qquad =Q_{ss^{\prime }s^{\prime \prime }}^{\langle uuuu\rangle }(\boldsymbol{k},\boldsymbol{k}^{\prime },t).\end{eqnarray}$$
The different signs
$s$
,
$s^{\prime }$
,
$s^{\prime \prime }\in \{-1,0,1\}$
characterize the eigenvalues
${\it\sigma}=0,\pm N\sin {\it\theta}$
, and
$\boldsymbol{k}^{\prime \prime }=-(\boldsymbol{k}+\boldsymbol{k}^{\prime })$
for a triadic interaction. The right-hand-side term involves fourth-order correlations. In the quasi-normal approximation, one expresses the latter in terms of second- and third-order correlations, through cumulants that represent the departure from an exactly Gaussian distribution (hence the ‘quasi-normal’ QN), written symbolically as
$$\begin{eqnarray}Q^{(IV)}=\underbrace{\langle uuuu\rangle }_{Q^{\langle uuuu\rangle }}-\underbrace{\langle uu\rangle \langle uu\rangle }_{Q^{(QN)}}\end{eqnarray}$$
(the exact tensorial development is detailed in Orszag Reference Orszag1973). Therefore (2.28) can be rewritten as
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\frac{\partial }{\partial t}+{\it\nu}\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)-N\!\left(s\sin {\it\theta}+s^{\prime }\sin {\it\theta}^{\prime }+s^{\prime \prime }\sin {\it\theta}^{\prime \prime }\right)\right)S_{ss^{\prime }s^{\prime \prime }}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)\nonumber\\ \displaystyle & & \displaystyle \qquad =Q_{ss^{\prime }s^{\prime \prime }}^{(QN)}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)+Q_{ss^{\prime }s^{\prime \prime }}^{(IV)}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)\end{eqnarray}$$
separating on the right-hand side the closed QN contribution from the contribution from fourth-order cumulants. The latter are expressed in terms of third-order moments, as did Orszag (Reference Orszag1969) in isotropic turbulence, by introducing an eddy-damping (ED) in the model to represent the scrambling effect of nonlinearities, to which we attribute a ‘functional’ damping role. This yields
$$\begin{eqnarray}Q_{ss^{\prime }s^{\prime \prime }}^{(IV)}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)=-\underbrace{\left({\it\eta}(k,t)+{\it\eta}(k^{\prime },t)+{\it\eta}(k^{\prime \prime },t)\right)}_{{\it\mu}_{kk^{\prime }k^{\prime \prime }}(t)}S_{ss^{\prime }s^{\prime \prime }}(\boldsymbol{k},\boldsymbol{k}^{\prime },t),\end{eqnarray}$$
introducing the ED rate
${\it\mu}_{kk^{\prime }k^{\prime \prime }}$
based on the definition of a characteristic turbulence damping timescale
${\it\eta}(k,t)=a_{0}(\int _{0}^{k}p^{2}E(p)\text{d}p)^{1/2}$
(Pouquet et al.
Reference Pouquet, Lesieur, André and Basdevant1975) computed from the kinetic energy spectrum
$E(k)$
. Classically, the constant introduced is fixed at
$a_{0}=0.36$
to recover the value of the Kolmogorov constant (André & Lesieur Reference André and Lesieur1977).
Although we have retained the linear stratification effects in the third-order moments equation (2.30), we further assume in the following model that these are negligible, with two supporting arguments: (a) linear stratification terms are already retained in second-order moments equations; (b) we want to keep the model mathematically as simple as possible. In other spectral models such as that by Canuto, Dubovikov & Dienstfrey (Reference Canuto, Dubovikov and Dienstfrey1997) for convective turbulence, the characteristic timescale used to model triple correlation damping did not include the explicit effect of stratification either, nor did they include the angular dependence of the two-point velocity spectra, which we shall retain in the present closure. The analytical solution for the triple correlations equation (2.30) is thus obtained assuming the model is Markovian (M), as
$$\begin{eqnarray}S_{ss^{\prime }s^{\prime \prime }}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)={\it\Theta}_{kk^{\prime }k^{\prime \prime }}(t)Q_{ss^{\prime }s^{\prime \prime }}^{(QN)}(\boldsymbol{k},\boldsymbol{k}^{\prime },t),\end{eqnarray}$$
with
$$\begin{eqnarray}{\it\Theta}_{kk^{\prime }k^{\prime \prime }}(t)=\frac{1-\text{exp}\left[-({\it\nu}\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}(t))t\right]}{{\it\nu}(k^{2}+k^{\prime 2}+k^{\prime \prime 2})+{\it\mu}_{kk^{\prime }k^{\prime \prime }}(t)}.\end{eqnarray}$$
The closure (2.32) expresses third-order moments in terms of second-order ones. It can be written in any frame of reference, and provides the transfer terms required to solve (2.13)–(2.17), which become closed integro-differential equations for two-point correlation spectra. In practice, one uses the polar-spherical Craya–Herring frame of reference.
2.3. Numerical implementation of the EDQNM model equations
The closed equations (2.24)–(2.27) are numerically solved by discretizing the
$\boldsymbol{k}$
-space using a polar-spherical mesh of discretized
$k$
and
${\it\theta}$
values. An important advantage of the model is that the discretization of
$k$
between
$k_{min}$
and
$k_{max}$
uses a logarithmic distribution with
$N_{k}$
modes, thus allowing a better representation of large scales than in the linear distribution in pseudo-spectral DNS. The angle
${\it\theta}$
in the model is uniformly discretized between
$[0,{\rm\pi}/2]$
(assuming vertical mirror symmetry in the absence of net vertical helicity) with
$N_{{\it\theta}}$
angles. One then needs to solve a system of
$N_{k}\times N_{{\it\theta}}$
coupled equations for the discretized
${\it\Phi}_{i}(k_{m},{\it\theta}_{n})$
,
$i=1,2,3$
, and
${\it\Psi}_{r}(k_{m},{\it\theta}_{n})$
,
$m=1\dots \,N_{k}$
,
$n=1\dots \,N_{{\it\theta}}$
.
As in Cambon & Jacquin (Reference Cambon and Jacquin1989) and Godeferd & Cambon (Reference Godeferd and Cambon1994), time-marching is done with a first-order Euler scheme over the time interval
$[t,t+{\rm\Delta}t]$
. The linear terms, viscous and stratification, are treated implicitly. For instance, for the toroidal energy mode
${\it\Phi}_{1}$
, the implicit treatment of viscosity yields
$$\begin{eqnarray}{\it\Phi}_{1}(t+{\rm\Delta}t)={\it\Phi}_{1}(t)\text{e}^{-2{\it\nu}k^{2}{\rm\Delta}t}+\frac{1-\text{e}^{-2{\it\nu}k^{2}{\rm\Delta}t}}{2{\it\nu}k^{2}{\rm\Delta}t}T^{{\it\Phi}_{1}}(t),\end{eqnarray}$$
where
${\rm\Delta}t$
is the timestep. The implicit treatment of other terms is similar although more complex due to the non-diagonal structure of the linear system, the corresponding eigenvectors being expressed by
$({\it\Phi}_{1},({\it\Phi}_{2}+{\it\Phi}_{3})/2\pm {\it\Psi}_{r},({\it\Phi}_{2}-{\it\Phi}_{3})/2)$
. Note that the different transfer terms are evaluated at
$t$
. At this point, an integration of third-order correlations over wavevectors is required to evaluate the transfers, from their expressions in (2.13)–(2.17). In practice, the integral is recast using a change of variables, replacing the three-dimensional integration over
$\boldsymbol{p}$
by a three-dimensional integration over the wavenumbers
$p$
,
$q$
and the angle
${\it\lambda}$
, as
$$\begin{eqnarray}\int _{\mathbb{R}^{3}}S(\boldsymbol{k},\boldsymbol{p})\text{d}^{3}\boldsymbol{p}=\iint _{{\it\Delta}}\frac{pq}{k}\left(\int _{0}^{2{\rm\pi}}S(\boldsymbol{k},{\it\lambda},p,q)\text{d}{\it\lambda}\right)\text{d}p\text{d}q,\end{eqnarray}$$
in which
${\it\lambda}$
denotes the angle rotating the plane of the triad (
$\boldsymbol{k},\boldsymbol{p},\boldsymbol{q}=-\boldsymbol{k}-\boldsymbol{p}$
) about
$\boldsymbol{k}$
. In so doing, the integration domain
${\it\Delta}$
for
$p,q$
as well as the Jacobian
$pq/k$
are the same as in the isotropic case. The integral can then be evaluated by the numerical method proposed by Leith (Reference Leith1971).
As in DNS, the stability of the EDQNM model numerical implementation imposes a threshold
${\rm\Delta}t$
, as well as a sufficient resolution in
$k$
(at least four wavenumbers per decade). The discretization in
${\it\theta}$
seems to impose a less stringent stability criterion, provided enough angles are considered for a reasonable representation of the angular dependence of the spectra. For the present computations,
$N_{k}=128$
and
$N_{{\it\theta}}=21$
. Using these figures, the computations are well resolved, checked by a grid convergence study.
2.4. Direct numerical simulations algorithm
In addition to computing solutions of USHT with the above EDQNM model, we perform DNS by solving the Navier–Stokes–Boussinesq equations (2.1)–(2.3) on a cubic box of size
$2{\rm\pi}$
, with triply periodic boundary conditions. We use a classical pseudo-spectral Fourier collocation method (see Orszag & Patterson Reference Orszag and Patterson1972; Rogallo Reference Rogallo1981; Soulard et al.
Reference Soulard, Griffond and Gréa2014) with phase-shift dealiasing. This allows us to keep all Fourier modes with wavenumber
$k$
satisfying
$k<\sqrt{2}N_{p}/3$
where
$N_{p}$
is the number of points in one direction to discretize the geometry. The code is parallelized using slab-shaped domain decomposition, and, for instance, the current simulations with
$1024^{3}$
points are run on 512 cores. The timestepping uses a third-order low-storage strong-stability-preserving Runge–Kutta scheme, with implicit viscous terms.
2.5. Solutions of the linearized equations
In the following, we shall compare systematically the predictions of EDQNM and DNS in the nonlinear regime to the solutions of the linearized equations (2.1)–(2.3). This permits us to assess the level of nonlinearity in the flow. It also provides a reference for comparing the departure between EDQNM and DNS to the nonlinear/linear departure. The linearized approach for computing two-point velocity spectral evolution is also sometimes called rapid distortion theory (RDT).
RDT solutions can be obtained in three ways: (a) by a semi-analytical solution of the spectral equations (2.24)–(2.27) in which the transfer terms are zeroed out; (b) by using the numerical implementation of the EDQNM model over the spherical discretization mesh, again discarding the integral transfer terms; (c) by using the DNS pseudo-spectral code in which nonlinear terms are removed, so that one obtains an RDT solution over a Cartesian discretization mesh. In the latter two ways, numerical inaccuracies can still be present due to discretization, unlike the solution (a). We have nonetheless chosen solution (b), after checking for a test set of parameters that sampling effects are negligible.
3. EDQNM results and comparison to DNS
We compare DNS and EDQNM simulations in order to evaluate the closure predictions in a range of parameters accessible to DNS. EDQNM can be used with a wider range of parameters not accessible to DNS, in terms of Reynolds number and of spectral range. We thus start with a parametric comparison of DNS and EDQNM (in § 3.2) that also permits us to assess the influence of the initial Froude number and mixing intensity on the evolution of USHT. The differences between DNS and EDQNM results suggest that nonlinear turbulent damping should depend on the intensity of stratification in the model, so that we improve the EDQNM closure with a simple yet efficient modification (following § 4).
Figure 3. Spherically averaged kinetic energy spectrum
$E(k)$
and buoyancy energy spectrum
$B(k)$
at
$t_{0}$
for
$1024^{3}$
DNS. See table 1 for the corresponding parameters.
3.1. Methodology and parametric cases
We build initial conditions for comparing EDQNM and DNS results as follows. Initial velocity and buoyancy fields
$\boldsymbol{u}$
and
${\it\vartheta}$
used to initialize USHT simulations at
$t=t_{0}$
come from a DNS of HIT with passive scalar
$c$
. Fully developed isotropic turbulence is obtained by running a simulation starting at
$t=0$
from energy and scalar random fields without dominant modal structure which could influence the development of the flow (e.g. the effect of initial spatial or frequency distribution on the evolution of Rayleigh–Taylor instability is discussed in Zhou, Robey & Buckingham Reference Zhou, Robey and Buckingham2003). We choose here Gaussian distributed fields with a von Karman spectral distribution
${\sim}k^{4}\exp [-(k/k_{peak})^{2}]$
. The peak wavenumber
$k_{peak}$
is chosen as large as possible in order to delay confinement effects due to the growth of the integral length scale. In the simulations,
$k_{peak}=40$
and the kinetic energy is
$\mathscr{K}(t=0)=0.75$
computed as
$\mathscr{K}=\int _{0}^{+\infty }E(k)\text{d}k$
. At
$t=t_{0}$
, USHT is obtained by applying gravity (
$\mathscr{A}g\neq 0$
) and the mean concentration gradient (
${\it\Gamma}\neq 0$
) parameterized by
$N$
. At this time, the scalar
$c$
becomes active, and triggers buoyancy
${\it\vartheta}$
, with
$\langle u_{i}{\it\vartheta}\rangle (t=t_{0})=0$
from isotropy. The time
$t_{0}$
is chosen such that turbulence is in a self-similar regime driven by large scales, for our
$k^{4}$
infrared spectrum at initial conditions. In terms of the initial eddy turnover frequency
$f_{ed}=(\mathscr{K}(0))^{1/2}k_{peak}$
, we have
$t_{0}\times f_{ed}\approx 20$
. The initialization of EDQNM is done at
$t=t_{0}$
using spectra extracted from DNS fields at this time. Unlike DNS, the model builds instantly fully developed triple correlations. The same kinematic viscosity
${\it\nu}$
and diffusivity
$\mathscr{D}$
is used both in DNS and EDQNM simulations, here set at
$2.5\times 10^{-3}$
.
For USHT simulations with unit Schmidt number, the non-dimensional parameters are the Reynolds and Froude numbers
$\mathit{Re}=\mathscr{K}^{2}/{\it\varepsilon}{\it\nu}$
and
$\mathit{Fr}={\it\varepsilon}/\mathscr{K}N$
, based on mean turbulent kinetic energy and its dissipation
${\it\varepsilon}$
. Initial conditions at
$t_{0}$
are also characterized by the mixing intensity
${\it\Lambda}$
, defined as the relative amount of kinetic energy
$\mathscr{K}$
and variance of buoyancy velocity
$\langle {\it\vartheta}{\it\vartheta}\rangle =\int _{0}^{+\infty }B(k)\text{d}k$
:
${\it\Lambda}=\langle {\it\vartheta}{\it\vartheta}\rangle /\mathscr{K}$
. We explore different values of
$\mathit{Fr}$
and
${\it\Lambda}$
by varying
$N$
and rescaling passive scalar
$c$
at
$t_{0}$
. This could lead to an issue with the resolution of the small scales of the scalar field if the amplitude of
$c$
becomes large. This is not the case here, since we use moderate
${\it\Lambda}$
, and we choose over-resolved initial conditions, anticipating the exponential energy growth.
The different simulations serving as the basis for direct comparison with EDQNM have a
$1024^{3}$
resolution. The different parameters are presented in table 1. The Reynolds numbers at
$t_{0}$
are small (
$\mathit{Re}\approx 3$
) due to the energy decay during the HIT phase. This is not an issue as it rapidly increases during USHT phase. The Froude numbers range between 0.4 and 1.6, allowing the study of different regimes in which linear buoyancy forces either dominate the short-time dynamics or are comparable to nonlinear effects. Exploring different mixing intensities
${\it\Lambda}$
permits the triggering of different transient phases, thus pushing further the comparison between DNS and EDQNM. In figure 3, we show the initial kinetic energy spectrum
$E$
, and the buoyancy energy spectra
$B$
at
$t_{0}$
, both for DNS and EDQNM. Note that, from
$k_{peak}=40$
at
$t=0$
, the maximum of the spectrum evolves towards
$k_{peak}\simeq 12$
at the end of the HIT phase at
$t_{0}$
, with a
$t^{-2/7}$
decay. Moreover, the spectra are very well resolved for simulations with
$1024^{3}$
points, only to be able to follow the USHT phase as long as possible with a fast increasing Reynolds number. At the small scales, the product between the Kolmogorov scale and the maximum wavenumber is
$k_{max}{\it\eta}=9.9$
at
$t_{0}$
for runs 1 to 9, and, for instance,
$k_{max}{\it\eta}=2.6$
at
$t^{\ast }=N(t-t_{0})=7$
in run 5. This shows that small scales are well resolved throughout the simulations, based on common DNS standards (
$k_{max}{\it\eta}\sim 1$
or 2 in Kaneda et al.
Reference Kaneda, Ishihara, Yokokawa, Itakura and Uno2003,
$k_{max}{\it\eta}>1.5$
in Pope Reference Pope2000).
Table 1. Parameters of the USHT simulations. The Reynolds number is
$\mathit{Re}=\mathscr{K}^{2}/{\it\varepsilon}{\it\nu}$
and the Froude number is
$\mathit{Fr}={\it\varepsilon}/N\mathscr{K}$
, where
$\mathscr{K}$
is the kinetic energy and
${\it\varepsilon}$
its dissipation. The initial mixing intensity is
${\it\Lambda}=\langle {\it\vartheta}{\it\vartheta}\rangle /\mathscr{K}$
.
Figure 4. Time evolution of: (a) Reynolds number
$\mathit{Re}$
, (b) Froude number
$\mathit{Fr}$
, (c) mixing intensity
${\it\Lambda}=\langle {\it\vartheta}{\it\vartheta}\rangle /\mathscr{K}$
; DNS runs 1 to 9.
3.2. Results
In this section, we discuss the results obtained from USHT simulations and we present DNS/EDQNM comparisons. It is important to recall that one DNS run provides a single realization of the flow field, whereas EDQNM results represent ensemble-averaged statistics. Statistics are obtained by space averages in DNS, but this means that statistical convergence may be poor, particularly at large scales. Unfortunately, this issue exists in the vast majority of DNS approaches.
Figure 4 shows the time evolution of non-dimensional numbers in runs 1 to 9 from DNS. The Reynolds number shown in figure 4(a) strongly increases in all the simulations, as the turbulence intensity increases due to the instability (as also observed on scalar field visualizations in figure 1). We clearly observe that initial conditions influence the time evolution of the Reynolds number: large Froude number cases with small
$\langle {\it\vartheta}{\it\vartheta}\rangle$
limit buoyancy production, leading to a short-time decrease of
$\mathit{Re}$
, as in HIT. At larger times, in all cases,
$\mathit{Re}$
strongly increases until
$t^{\ast }\approx 4$
; subsequently, the growth rate is reduced, indicating a change of regime, which we attribute to the appearance of the self-similar state.
In figure 4(b,c) for the time evolution of Froude number and mixing intensity, we observe a short-time transient in which the different curves evolve from well-separated initial values. After
$t^{\ast }>4$
, the Froude number curves (figure 4
b) join at
$\mathit{Fr}\approx 0.6$
, after reaching a minimum at
$t^{\ast }\simeq 3$
. This means that during the initial transient, buoyancy effects are very active in driving the flow, before nonlinear turbulent interactions enter more strongly into play. At
$t^{\ast }\simeq 6$
, the evolution seems to steady to 0.6. A similar behaviour is observed for the mixing intensity
${\it\Lambda}$
in figure 4(c), where a universal partition value
${\it\Lambda}\approx 1.8$
is reached from
$t^{\ast }\simeq 3$
after the initial transient, with a very slow evolution afterwards. As shown by mixing intensity curves, the energy from mean stratification injected into the system is equally distributed between fluctuations of velocity and buoyancy. It can be easily shown that
${\it\Lambda}=2$
corresponds exactly to the ratio obtained from linearized equations if we neglect nonlinear and diffusive terms. Therefore, the fact that
${\it\Lambda}$
is not exactly 2 is the imprint of nonlinear effects.
In the following, we shall use the normalized kinetic energy
$\mathscr{K}^{\ast }(t)=\mathscr{K}(t)/(\mathscr{K}+\langle {\it\vartheta}{\it\vartheta}\rangle )(t_{0})$
; likewise for
$\langle {\it\vartheta}{\it\vartheta}\rangle ^{\ast }(t)$
and
$\langle {\it\vartheta}u_{3}\rangle ^{\ast }(t)$
. Figure 5 shows the time evolution of kinetic energy
$\mathscr{K}^{\ast }(t)$
, variance of buoyancy
$\langle {\it\vartheta}{\it\vartheta}\rangle ^{\ast }(t)$
and vertical buoyancy flux
$\langle {\it\vartheta}u_{3}\rangle ^{\ast }(t)$
for run 5, from DNS and EDQNM. The results from the linear model (RDT) starting from the same initial conditions are also represented in order to show the importance of nonlinear effects. If viscosity and diffusion are discarded, RDT has a simple long-time asymptotic solution – for instance the kinetic energy
$\mathscr{K}(t)\sim \text{e}^{2Nt}$
– which is close to what is observed in figure 5(a) for the RDT numerical solution. The short-time behaviour of RDT matches that of DNS and EDQNM up to
$t^{\ast }\approx 2$
, after which it fails to capture the inflection of the curves due to nonlinear corrections. The global trend in DNS and EDQNM is the same, with a strong growth of the amplitude for all three quantities plotted in figure 5. The EDQNM model captures the damping of energy growth and flux due to nonlinear interactions. The solutions for DNS and EDQNM agree as long as nonlinear effects are not too strong. At large
$t$
, the model clearly underestimates the values of kinetic and buoyancy energy as well as buoyancy flux, compared to DNS. This tendency is observed in all cases of table 1 (not shown here for the sake of brevity).
Figure 5. Time evolution of: (a) kinetic energy
$\mathscr{K}^{\ast }$
, (b) buoyancy variance
$\langle {\it\vartheta}{\it\vartheta}\rangle ^{\ast }$
, (c) buoyancy vertical flux
$\langle {\it\vartheta}u_{3}\rangle ^{\ast }$
for run 5 from DNS, EDQNM and the linear RDT approximation.
Figure 6 allows a comparison of the EDQNM and DNS energy spectra and an evaluation of their scale-by-scale differences. The departure at small scales between EDQNM and DNS spectra is expected, as was observed in previous studies (Cambon et al.
Reference Cambon, Mansour and Godeferd1997; Godeferd & Staquet Reference Godeferd and Staquet2003). It may be attributed to both the EDQNM model in which small-scale dynamics is built-in (as in all statistical two-point spectral closures, see Chen et al.
Reference Chen, Doolen, Herring, Kraichnan, Orszag and She1993), but also in part to DNS large-wavenumber truncation that produces a small cusp in the tail of the spectrum. Here, as shown in the figure, the departure seems to appear first in the energetic wavenumber range, leading us to believe that EDQNM overestimates the large-scale transfer, but one could also point out that the infrared truncation of DNS spectra may affect its dynamics through box confinement. In order to discard this latter hypothesis, we have carried out a numerical study on confinement effects at large scales using different grid resolutions. It appears clearly that before
$t^{\ast }=6$
, integrated quantities are not sensitive to confinement. Therefore, we conclude that the dynamics modelled in EDQNM is underestimated due to a weakness of the closure.
Figure 6. Kinetic energy spectra comparison between EDQNM and DNS for run 5 at
$t^{\ast }=0$
, 3, 6.
Finally, can one say that a self-similar regime for USHT appears in our DNS results? In figure 5, we have plotted an additional curve corresponding to the expected exponential self-similar growth
${\sim}\text{e}^{{\it\beta}Nt}$
, with
${\it\beta}=4/7$
obtained from (1.2) with
$s=4$
, since the infrared scaling in our initial spectra is
$k^{4}$
. One could consider the last inflection observed in figure 5(a,b) for the energies to indicate the appearance of a self-similar regime, and deduce the actual value of the growth rate from results at
$t^{\ast }=6$
. Unfortunately, the duration of the late phase is very short in our simulations and we know that confinement effects appear quickly, so this clearly jeopardizes any conclusions. Increasing the available spectral range of the simulation would certainly give slightly better results, although at a tremendous computational cost.
4. Correction for explicit
$N$
-dependent triple correlations
Results presented in the previous section suggest that the EDQNM model captures the linear and nonlinear stages in USHT evolution, although with an incorrect nonlinear estimate. In § 4.1, we propose a simple modification of the eddy-damping term. In § 4.2, we compare the improved model and show that it leads to much better agreement with DNS and also that it captures the important anisotropy of USHT.
4.1. Methodology of the correction
Arguments for turbulent decorrelation timescales used in two-point statistical models are derived from phenomenological analyses of the interaction between large and small scales, and the choice of the ED timescale has to include enough correct physics, as discussed by Zhou, Matthaeus & Dmitruk (Reference Zhou, Matthaeus and Dmitruk2004). Two timescales may be considered (see e.g. Zhou Reference Zhou2010). First, related to the local straining effect of turbulent structures on each other, the nonlinear timescale
${\it\tau}_{nl}$
for turbulent interaction permits the recovery of the correct Kolmogorov spectral scaling. In the above model, we used
${\it\eta}(k)^{-1}$
. Second, one can also consider the sweeping of small structures by large ones, but the use of such a timescale based on a mean or root-mean-square (r.m.s.) velocity yields a non-Kolmogorov scaling of the energy spectrum even in isotropic turbulence (
$k^{-3/2}$
instead of
$k^{-5/3}$
). In USHT phenomenology, it is clear that vertical buoyant events may trigger sweeping-like phenomena in the turbulent dynamics and involve an external timescale
${\it\tau}_{ex}$
. We therefore propose to use a hybrid timescale
$1/({\it\tau}_{ex}^{-1}+{\it\tau}_{nl}^{-1})$
in the EDQNM model, as a pragmatic correction that does not require a complete rewriting of the closure, choosing the external timescale
${\it\tau}_{ex}\propto N^{-1}$
.
More precisely, recall that, in deriving the EDQNM model, we neglected stratification terms appearing in (2.28) for the dynamics of triple correlations, that is we did not retain the buoyancy timescale in the evolution of these statistical moments. We propose to modify the closure in order to recreate buoyancy effects in the transfer terms, without resorting to solving the full equation (2.28). The conclusions of § 3 establish that the transfer is overestimated in the basic version of EDQNM. We simply close buoyancy terms in fourth-order cumulants by introducing a stratification-related damping term based on the timescale
$N^{-1}$
. Equation (2.28) then becomes
$$\begin{eqnarray}\!\left(\partial _{t}+{\it\nu}\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}^{N}(t)\right)S_{ss^{\prime }s^{\prime \prime }}\!\left(\boldsymbol{k},\boldsymbol{k}^{\prime },t\right)=Q_{ss^{\prime }s^{\prime \prime }}^{(QN)},\!\left(\boldsymbol{k},\boldsymbol{k}^{\prime },t\right)\end{eqnarray}$$
with a new ED rate
${\it\mu}_{kk^{\prime }k^{\prime \prime }}^{N}(t)={\it\eta}^{N}(k,t)+{\it\eta}^{N}(k^{\prime },t)+{\it\eta}^{N}(k^{\prime \prime },t)$
accounting both for fourth-order cumulants and stratification:
$$\begin{eqnarray}{\it\eta}^{N}(k,t)=a_{0}\left(\int _{0}^{k}p^{2}E(\,p,t)\text{d}p\right)^{1/2}+a_{1}N.\end{eqnarray}$$
The amplitude of the additional term proportional to
$N$
is adjusted using a new constant
$a_{1}$
, to be determined. (We recall that the value of
$a_{0}$
was defined once and for all when EDQNM was developed, and has remained unchanged ever since.) Although a bit coarse, since it does not take into account the explicit dependence of the linear stratification operator on
${\it\theta}$
, this closure ensures the realizability of the model. A similar simplified model is also used for rotating turbulence (Cambon et al.
Reference Cambon, Mansour and Godeferd1997).
In the classical EDQNM approach, the time dependence of the ED rate is neglected during the ‘Markovianization’ process. However, in the modified closure, Markovianization should also consistently account for the switch from isotropic turbulence at
$t<t_{0}$
to stratified turbulence at
$t>t_{0}$
which induces a discontinuity in our ED rate
${\it\mu}_{kk^{\prime }k^{\prime \prime }}^{N}$
. Therefore, the characteristic time
${\it\Theta}_{kk^{\prime }k^{\prime \prime }}(t)$
for
$t>t_{0}$
needs to be changed after taking into account Markovianization to
$$\begin{eqnarray}\displaystyle {\it\Theta}_{kk^{\prime }k^{\prime \prime }}(t) & = & \displaystyle \text{exp}\!\left[-({\it\nu}\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}^{N}(t))\times (t-t_{0})\right]\nonumber\\ \displaystyle & & \displaystyle \times ~\frac{1-\text{exp}\left[-({\it\nu}\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}(t_{0}))\times t_{0}\right]}{{\it\nu}\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}(t_{0})}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1-\text{exp}\left[-({\it\nu}\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}^{N}(t))\times (t-t_{0})\right]}{{\it\nu}\!\!\left(k^{2}+k^{\prime 2}+k^{\prime \prime 2}\right)+{\it\mu}_{kk^{\prime }k^{\prime \prime }}^{N}(t)}.\end{eqnarray}$$
At
$t=t_{0}$
, one recovers the classical expression for the characteristic time of unforced turbulence.
In order to calibrate the new constant
$a_{1}$
, we perform optimization using a cost function
$\mathscr{J}(a_{1})$
based on the time-integrated relative difference between DNS and the new EDQNM model energies:
$$\begin{eqnarray}\displaystyle \mathscr{J}(a_{1}) & = & \displaystyle \int _{t_{0}}^{t_{1}}\left(\left(\frac{\mathscr{K}^{DNS}(s)-\mathscr{K}^{EDQNM}(s,a_{1})}{\mathscr{K}^{DNS}(s)}\right)^{2}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\left(\frac{\langle {\it\vartheta}{\it\vartheta}\rangle ^{DNS}(s)-\langle {\it\vartheta}{\it\vartheta}\rangle ^{EDQNM}(s,a_{1})}{\langle {\it\vartheta}{\it\vartheta}\rangle ^{DNS}(s)}\right)^{2}\right)\text{d}s\end{eqnarray}$$
The upper time bound
$t_{1}$
is such that
$N(t_{1}-t_{0})=5$
to prevent confinement effects from influencing the determination of
$a_{1}$
. We minimize
$\mathscr{J}(a_{1})$
for runs 1 to 9 using a simple Newton algorithm, each iteration necessitating a new EDQNM simulation. The different values of
$a_{1}$
corresponding to the optimization procedure for each case are given in table 2. The convergence is achieved typically within five iterations for an accuracy of
$\Vert {\rm\Delta}\mathscr{J}\Vert /\mathscr{J}<10^{-3}$
and corresponding to a variation of the constant of
$\Vert {\rm\Delta}a_{1}\Vert /a_{1}<10^{-5}$
. The different constants
$a_{1}$
for each case are in the interval
$a_{1}\in [0.15,0.4]$
. It is very reassuring to find values of
$a_{1}$
not too large, especially since they are smaller than the maximum eigenvalue of the linear operator. The variations of
$a_{1}$
are not important from one case to the other but we can notice a sensitivity to the mixing intensity
${\it\Lambda}$
which is not taken into account in the expression for the ED term. We set the final value at the average,
$a_{1}=0.28$
.
Table 2. Values of the constant
$a_{1}$
obtained from minimizing the cost function
$\mathscr{J}$
for each run.
4.2. Updated results
Let us now call the modified model EDQNMc, and compare again to DNS data. Each case described in table 1 is computed with the updated model. We choose to show in figure 7 the evolution of kinetic energy for runs 1 to 9. As expected, the predictions of EDQNMc agree very well with DNS for all cases throughout the time evolution, up to
$t^{\ast }=6$
. At this time the departure from DNS is at most 20 % whereas it was up to 50 % in the unmodified model. This improvement and the relatively small departure is a very good result, considering the fact that the model is able to reproduce a complete evolution over six stratification frequency periods. The transitional state is therefore also well captured with the updated model. Other quantities such as potential energy and buoyancy flux agree as well with DNS (not shown). Note that a uniformly good result was not guaranteed from the start as the constant
$a_{1}=0.28$
is not the optimized constant for each run. This confirms that the dynamics of triple correlations in USHT is significantly modified by the presence of stratification.
Figure 7. Kinetic energy evolution for runs 1 to 9 renormalized by initial sum of potential and kinetic energy from DNS, EDQNM and the updated model EDQNMc.
What is the result of the scale-by-scale comparison between EDQNMc and DNS spectra? We choose to present results for run 5 in figure 8, corresponding to intermediate values of the initial Froude number
$\mathit{Fr}=0.8$
and for mixing intensity
${\it\Lambda}=2$
. We present the kinetic and potential energy spectra in figure 8(a,b) at times
$t^{\ast }=0$
, 3 and 6, and the buoyancy flux spectrum in figure 8(c), at times
$t^{\ast }=0.2$
, 3 and 6, since the flux at
$t^{\ast }=0$
is zero. The agreement between EDQNMc and DNS is impressive in the large scales, especially considering the long duration of the simulation, and the remaining small departure with DNS at the latest time lies within the error in the DNS data due to confinement, although one could also mention the lack of reliable backscatter in the EDQNM closure, which may be present for
$k^{4}$
initial spectra. In that sense, a test with initial
$k^{4}$
spectra (Batchelor type) is more challenging for EDQNM than one with
$k^{2}$
(Saffman type) for instance. Other comparisons with DNS using smaller slopes have been conducted and show an even better DNS/EDQNM agreement. This confirms the important role of backscatter in the results presented.
Figure 8. Evolution of spectra for run 5 at
$\mathit{Fr}=0.8$
and mixing intensity
${\it\Lambda}=2$
. (a) Kinetic energy spectrum; (b) buoyancy energy spectrum; (c) buoyancy flux.
The inertial-range scaling is very well reproduced by the model also. At the small scales, the match is not as good as at the large scales, especially at the end of the simulation. As mentioned previously, EDQNM is not expected to compare very well with DNS in the dissipative range, of less importance in USHT than the highly anisotropic largest scales, as we shall see in figure 12.
We investigate further DNS/EDQNM comparisons by presenting measurements of anisotropy in USHT. We begin with the deviatoric part of the Reynolds stress tensor,
$\unicode[STIX]{x1D623}_{ij}=\langle u_{i}u_{j}\rangle /\langle u_{i}u_{i}\rangle -(1/3)$
, of which we retain only the component
$\unicode[STIX]{x1D623}_{33}$
, since
$\unicode[STIX]{x1D623}_{ii}=0$
from incompressibility, implying
$\unicode[STIX]{x1D623}_{11}=\unicode[STIX]{x1D623}_{22}=-\unicode[STIX]{x1D623}_{33}/2$
in axisymmetric flows. The evolution of
$\unicode[STIX]{x1D623}_{33}$
is shown in figure 9(a). The linear RDT result evolves from the initial zero value to the asymptotic
$2/3$
value, which can be analytically computed. At
$t^{\ast }\approx 2$
, the linear approximation is not valid anymore, and
$\unicode[STIX]{x1D623}_{33}$
begins to decrease from its peak value of approximately 0.47, as if return to isotropy were initiated. However, the later trend shows another inflection in the evolution of
$\unicode[STIX]{x1D623}_{33}$
, which seems to settle to a value close to 0.35. The EDQNMc model prediction is very close to the DNS evolution.
Figure 9. Statistical characterization of anisotropy in run 5. Evolution of (a) component
$\unicode[STIX]{x1D623}_{33}=\langle u_{3}u_{3}\rangle /\langle u_{i}u_{i}\rangle -(1/3)$
of the deviatoric part of the Reynolds-stress tensor; (b)
$\sin ^{2}{\it\gamma}$
.
We also investigate the anisotropy of the scalar, which we quantify using the dimensionality parameter
$$\begin{eqnarray}\sin ^{2}{\it\gamma}=\frac{\displaystyle \int _{0}^{{\rm\pi}}\sin ^{2}{\it\theta}\int _{0}^{\infty }{\it\Phi}_{3}(k,{\it\theta})\text{d}k\text{d}{\it\theta}}{\displaystyle \int _{0}^{{\rm\pi}}\int _{0}^{\infty }{\it\Phi}_{3}(k,{\it\theta})\text{d}k\text{d}{\it\theta}}\end{eqnarray}$$
for two reasons. First, we need to construct a measure of anisotropy based on a scalar field. A way to do this is the above definition, expressing the stretching of scalar structures along the vertical direction (see Gréa Reference Gréa2013), which is clear in figure 1. Second, the Reynolds stress tensor anisotropy evaluated from
$\unicode[STIX]{x1D623}_{33}$
reflects the anisotropy due to the difference between
${\it\Phi}_{1}(\boldsymbol{k})$
and
${\it\Phi}_{2}(\boldsymbol{k})$
, namely the polarization anisotropy, integrated over wavespace, but very much less so the directional dependence of these spectra, i.e. their dependence upon
${\it\theta}$
(see details in e.g. Sagaut & Cambon Reference Sagaut and Cambon2008 and Delache, Cambon & Godeferd Reference Delache, Cambon and Godeferd2014). This dependence is included in the definition of
$\sin ^{2}{\it\gamma}$
. The dimensionality parameter evolves in time as shown in figure 9(b) for run 5. Again, the model reproduces correctly the global trends of DNS. The anisotropy of the scalar structures seems to be linked to that of the Reynolds stress tensor, since the overall behaviour of the curve is similar to that of
$\unicode[STIX]{x1D623}_{33}$
: it increases to a maximum 0.8 at
$t^{\ast }\approx 2$
, then recovers a more moderate value
$\sin ^{2}{\it\gamma}\approx 0.75$
. The agreement between DNS and EDQNM is better for the anisotropy of the Reynolds stress tensor than for the scalar field. This is explained as follows: the anisotropy of the Reynolds stress tensor can be decomposed into contributions from polarization and from directional anisotropies, and, in USHT, it appears that polarization anisotropy plays a very important role, and is linked to the forcing of vertical velocity components by buoyancy. Thus polarization anisotropy is dominant in the Reynolds stress tensor. On the contrary, the dimensionality parameter expresses only directional anisotropy, as seen from its definition (4.5). Since we have chosen not to explicitly include the detailed
${\it\theta}$
-dependence of stratification operators in the eddy damping in the model (§ 4.1), we do not expect to represent it very accurately in the dimensionality parameter. Still, although the DNS/EDQNMc curves for
$\sin ^{2}{\it\gamma}$
do not collapse perfectly in figure 9(b), the similarity is rather good.
5. Unstably stratified homogeneous turbulence at very high Reynolds number
In this section, we present results of a long EDQNM simulation at very high Reynolds number, impossible to reach with DNS. Our goal is to reach the USHT self-similar state and provide a prediction of the energy growth rate, thanks to the anisotropic EDQNM closure. The simulation corresponding to run 5 with
${\it\nu},~\mathscr{D}=2.5\times 10^{-3}$
is continued until
$t^{\ast }=14$
, where the turbulent Reynolds number reaches
$\mathit{Re}\approx 10^{5}$
.
In figure 10(a), the time evolution of Reynolds number is displayed in a semi-logarithmic representation (the curve for kinetic energy is similar). What draws particularly attention is that the self-similar evolution, with an exponential growth
$\mathit{Re}(t)\propto \mathscr{K}(t)\propto \text{e}^{{\it\beta}Nt}$
, clearly appears from
$t^{\ast }\sim 7$
. While this state is hardly established at the end of the DNS, it can be clearly observed through EDQNM simulations. It is important to note that, after the linear transient, there exists another intermediate state between the short-time and the asymptotic long-time regimes, during which the Reynolds number (or energy) growth rate is not as large as the asymptotic value. Therefore, in DNS that are limited in time, it is not enough to observe the flow immediately after the linear regime. The slope corresponding to the growth parameter
${\it\beta}$
agrees with the theoretical prediction
${\it\beta}=4/7$
for
$k^{4}$
infrared spectra.
Figure 10. Time evolution of (a) the Reynolds number
$\mathit{Re}$
in semi-logarithmic representation, (b) the Froude number
$\mathit{Fr}$
and (c) the mixing intensity
${\it\Lambda}$
, all for run 5.
The attenuation of the growth of kinetic energy after the linear transient can be analysed from a different point of view by examining the Froude number evolution in figure 10(b). The evolution of the Froude number during the transient phase is very rapid. After the initial short-time decay, it reaches a maximum
$\mathit{Fr}\simeq 0.7$
at
$t^{\ast }\approx 6$
, and is then reduced to
${\approx}0.48$
when entering the self-similar regime during which it remains almost stationary. This rapid phenomenon indicates a re-adjustment of the dissipation rate in the transient phase after the appearance of stratification effects at
$t_{0}$
. The latter acting rapidly in the large scales, there is a delay between the re-adjustment of kinetic energy and that of dissipation over a timescale proportional to the extent of the inertial subrange, since
$\mathit{Fr}={\it\varepsilon}/(N\mathscr{K})$
. Moreover, the precise evolution of
$\mathit{Fr}$
in this ‘grey zone’, neither linear nor fully self-similar, strongly depends upon the specific initial conditions (energy ratios, non-dimensional parameters, shape of the spectrum, etc.).
The evolution of the mixing intensity
${\it\Lambda}$
is presented in figure 10(c). The slow decrease of the relative importance of buoyancy fluctuations compared to velocity fluctuations may be attributed to the conversion of poloidal kinetic energy generated by unstable stratification to toroidal kinetic energy by nonlinear transfers. In fact, contrary to poloidal, toroidal kinetic energy cannot in return trigger buoyancy fluctuations through linear production.
Finally, we come back to the questions we raised in the Introduction about the presence of different equilibrium ranges in the spectra of USHT, illustrated by the sketch of figure 2. Are we able to identify the corresponding scaling at the very high Reynolds number reached by the self-similar evolution predicted in our last EDQNMc simulation? We plot the kinetic energy spectrum
$E(k)$
, the buoyancy spectrum
$B(k)$
and vertical cross-spectrum
$F(k)$
in figure 11(a–c) at the latest time available in self-similar regime,
$t^{\ast }=14$
, and the evolution from
$t^{\ast }=11$
in figure 11(d–f). The inertial subrange exhibits a clear
$k^{-5/3}$
scaling for
$E(k)$
and
$B(k)$
and a
$k^{-7/3}$
one for
$F(k)$
, as predicted by Lumley (Reference Lumley1964), that spreads out over nearly two decades. The initial
$k^{4}$
infrared subrange seems to persist during the simulation, resulting from combined effects of linear terms and backscatter. It is also interesting to introduce the integral and the Ozmidov scales corresponding respectively to wavenumbers
$k_{\ell }$
and
$k_{O}$
. The ratio of these two quantities is related to the Froude number as
$k_{\ell }/k_{O}=\mathit{Fr}^{3/2}$
(see for instance Lesieur Reference Lesieur2008). In the self-similar regime,
$\mathit{Fr}=0.5$
leads to
$k_{\ell }/k_{O}\approx 1/3$
. Therefore, the two scales are relatively close to each other. This should indicate also that the return to isotropy in the inertial range should be fast, spreading around the small bump which can be seen around the maxima of the spectra. Upon reporting the exact values of the length scales on spectra of figure 11(a–c),
$k_{\ell }$
being associated with the maximum of
$E(k)$
(also
$B(k)$
and
$F(k)$
), it appears that
$k_{O}$
falls at the lower edge of the inertial subrange where stratification and anisotropic effects have to be taken into account, as expected. In figure 11(d–f), we also show the time evolution of energy spectra during the self-similar phase. It can be verified that the maxima for the different spectra follow a
$k^{-3}$
line as expected from the theory of Poujade (Reference Poujade2006) adapted to USHT. We believe that the present spectral characterization predicted by EDQNMc has not yet been observed in the highest-resolved DNS of buoyancy-induced turbulence as Rayleigh–Taylor flows (e.g. Cook & Cabot Reference Cook and Cabot2006).
Figure 11. (a–c) Energy, scalar and cross-spectra as a function of the wavenumber
$k$
for run 5 at
$t^{\ast }=14$
with EDQNMc. The indicated wavenumbers corresponding to the integral and Ozmidov scales,
$k_{\ell }$
and
$k_{O}$
. (d–f) Evolution of energy, scalar and cross-spectra at times
$t^{\ast }=11$
, 12, 13, 14, as a function of the wavenumber
$k$
for the same run 5.
What happens to the anisotropy at different scales in the self-similar state? Angular spectra of energy
${\it\Phi}_{1}$
,
${\it\Phi}_{2}$
,
${\it\Phi}_{3}$
and of the poloidal velocity/buoyancy cross-correlation
${\it\Psi}_{r}$
are shown in figure 12. They permit an estimation of the scale-by-scale anisotropy. A first observation in all spectra is that anisotropy in USHT is mostly confined to the large-scale range. The progressive restoration of isotropy occurs at the beginning of the inertial subrange, after
$k_{O}$
for increasing
$k$
. This is attested first by spectra almost independent of
${\it\theta}$
in this range, and second by the fact that
${\it\Phi}_{1}$
and
${\it\Phi}_{2}$
reach common levels. In their experimental realization of USHT, Thoroddsen et al. (Reference Thoroddsen, Van Atta and Yampolsky1998) seemed also to observe maximal anisotropy in the largest spectral scales by comparing horizontal and vertical spectra of two-point velocity correlations, but some anisotropy was still observable at smaller scales, although not as large.
Figure 12. Angular spectra: (a)
${\it\Phi}_{1}$
, (b)
${\it\Phi}_{2}$
, (c)
${\it\Phi}_{3}$
and (d)
${\it\Psi}_{r}$
are represented for EDQNM simulation run 5 at
$t^{\ast }=14$
as a function of
$k$
for different angles
${\it\theta}=0$
,
${\rm\pi}/8$
,
${\rm\pi}/4$
,
$3{\rm\pi}/8$
,
${\rm\pi}/2$
. Insets: deviatoric part of angular spectra
${\rm\Delta}{\it\Phi}=100\times ({\it\Phi}-\langle {\it\Phi}\rangle _{{\it\theta}})/\langle {\it\Phi}\rangle _{{\it\theta}}$
with
$\langle {\it\Phi}\rangle _{{\it\theta}}=(2/{\rm\pi})\int _{0}^{{\rm\pi}/2}{\it\Phi}\text{d}{\it\theta}$
as a function of
${\it\theta}$
for different wavenumbers
$k=0.1$
, 1, 10.
The second important anisotropic feature is that maximal energy is observed at
${\it\theta}={\rm\pi}/2$
, i.e. for horizontal wavenumbers, in
${\it\Phi}_{2}$
, whereas the anisotropy in the
${\it\Phi}_{1}$
spectrum is not only less monotonic with
${\it\theta}$
but the anisotropy is almost reversed, with lowest energy found in the horizontal at
${\it\theta}={\rm\pi}/2$
. Spectra for
${\it\Phi}_{2}$
,
${\it\Phi}_{3}$
and
${\it\Psi}_{r}$
are maximum at
${\it\theta}={\rm\pi}/2$
, because of the predominance of the linear buoyancy production term (see (2.25)–(2.27)). This matches USHT phenomenology with elongated structures along the vertical direction. However, the dynamics of
${\it\Phi}_{1}$
is mainly driven by nonlinear interactions, as is clear in (2.24) in which no linear production term appears.
Other interesting features of spectral anisotropy are: the angular dependence of the buoyancy density spectrum
${\it\Phi}_{3}$
is less than that of
${\it\Phi}_{2}$
, as also observed in stably stratified homogeneous turbulence (Godeferd & Cambon Reference Godeferd and Cambon1994); although maximal at
$k<k_{O}$
, some anisotropy is still observable throughout the scales in the flux, a hint that buoyancy acts at all scales, even weakly; in the infrared subrange, spectra at
${\it\theta}={\rm\pi}/2$
always scale as
$k^{4}$
when time evolves, while at other orientations
${\it\theta}$
they are subject to a nonlinear redistribution of energy, as noted by Soulard et al. (Reference Soulard, Griffond and Gréa2014).
The insets in figure 12 permit an evaluation of the variation with
${\it\theta}$
of energy at three different wavenumbers:
$k=0.1$
in the infrared range, and
$k=1$
at the beginning and
$k=10$
at the middle of the inertial range. The departure about the mean over
${\it\theta}$
of the energy at this scale is plotted. The curves at
$k=10$
confirm our previous observation that anisotropy is comparably weak at this wavenumber. Anisotropy amplitude is almost the same at
$k=0.1$
and
$k=1$
for
${\it\Phi}_{2}$
,
${\it\Phi}_{3}$
and
${\it\Psi}_{r}$
, and energy is accumulated close to horizontal wave vectors
${\bf\theta}={\rm\pi}/2$
. Only for
${\it\Phi}_{1}$
does anisotropy evolve differently with accumulation at
${\bf\theta}=0$
in the infrared subrange, although with an amplitude less than for the other spectra, as already noted. These insets in figure 12 suggest that simple phenomenological arguments for explaining the anisotropy distribution can be put forward for
${\it\Phi}_{2}$
and
${\it\Phi}_{3}$
, but a more complex analysis is required in order to interpret anisotropy variations in
${\it\Phi}_{1}$
. We note also that obtaining angular-dependent spectra in DNS such as those presented in figure 12 raises the problem of lack of sampling in the small wavenumbers.
6. Conclusion
In this work, we explore USHT from different initial conditions, characterized by a Froude number and mixing intensity. The study by DNS of the late-time dynamics of turbulent statistics and the self-similar regime has limitations because of confinement effects induced by the finite size of the computational box. In particular, boundary conditions constrain the growth of the integral scale. In order to overcome this difficulty and to reach higher Reynolds number regimes, we propose a spectral analysis based on an anisotropic EDQNM model. In a first approach, we consider a closure in which the dynamics for triple correlations does not account at all for buoyancy effects, assuming to first order that the corresponding production terms should modify significantly only the second-order statistics. The resulting model agrees with DNS at short times and gives the correct growth rate of the self-similar regime. However, nonlinear third-order transfers are overestimated during the transient phase from linear to asymptotic regimes, leading to an overly attenuated growth of turbulent energy. This point illustrates the difficulty of predicting USHT as one needs to capture correctly the linear effect of the production term and also to evaluate its importance on the nonlinear interactions between the different scales. We therefore adapt the model in an updated version, in which we introduce a heuristic correction by taking into account stratification effects in the ED term which expresses decorrelation at the level of nonlinear triple correlations. This simple method, while keeping the model tractable, leads to much better agreement with DNS, even if the anisotropy of the flow is not exactly reproduced for the scalar field. The model is nonetheless able to reproduce closely the spectral energy distribution in the energetically important spectral subrange, in the infrared and in the inertial one. During the comparisons between EDQNM and DNS results, we have identified the quantitative effect of confinement of the large-scale structures by the limited size of the computational domain. This is an important issue for DNS, which ought to be undertaken at high resolution and at very long times in order to capture USHT self-similar regimes.
The spectral EDQNM closure does not have this drawback. Spectral EDQNM simulations with turbulent Reynolds number as large as
$\mathit{Re}\approx 10^{5}$
are presented, well beyond the current possibilities of DNS. By comparison, our
$1024^{3}$
DNS simulations with the same initial conditions have to be stopped at
$\mathit{Re}\approx 10^{3}$
, and the
$3072^{3}$
simulations of Rayleigh–Taylor instability by Cook & Cabot (Reference Cook and Cabot2006) reached
$\mathit{Re}\approx 4600$
at a formidable computational price (
${\sim}10^{6}$
CPU hours), compared with EDQNM cost (
${\sim}10^{2}$
CPU hours). The late-time self-similar regime of USHT can be observed with a much better accuracy thanks to the model. The results obtained are of two kinds. First, they support the theoretical formula (1.2) for the growth rate proposed by Soulard et al. (Reference Soulard, Griffond and Gréa2014). The model also predicts a strong anisotropy at large scales and isotropy of scales in the inertial subrange as well as in the dissipative one. This is a rather unexpected result, since one could expect that, at least from local interactions, scales in the lower inertial range would inherit part of the large-scale anisotropy. The transition from large to small scales is therefore characterized by a fast return to isotropy, moving from the integral to the Ozmidov length scales. The separation between those two scales depends on the Froude number and is rather short in our USHT simulations. Nonetheless, we are still able to observe a trend towards the creation of such an intermediate subrange with non-Kolmogorov scaling. An interesting related question is whether such a zone will be observed in different asymptotic states of USHT obtained by changing initial conditions through the slope of the infrared spectra. This is left for future work.
Acknowledgements
A.B. and B.-J.G. would like to thank Dr O. Soulard for discussions about theoretical predictions of the self-similar states of USHT. EDQNM and DNS computations were performed at the computing center of CEA.
Appendix A. Derivation of the self-similar growth rate
${\it\beta}$
in USHT
Soulard et al. (Reference Soulard, Griffond and Gréa2014) propose an analytic formula for the growth rate of turbulent quantities in self-similar USHT with time-varying buoyancy frequency. In this appendix, we show that by similar arguments it is possible to derive an analogous formula, (1.2), for USHT with constant
$N$
.
First, if the flow is in a self-similar regime, the kinetic energy, dissipation and integral length scale evolve as
$$\begin{eqnarray}\mathscr{K}(t)\sim \text{e}^{{\it\beta}Nt},\quad {\it\varepsilon}(t)\sim \text{e}^{{\it\beta}Nt},\quad \ell (t)=\frac{\mathscr{K}^{3/2}(t)}{{\it\varepsilon}}\sim \text{e}^{{\it\beta}Nt/2}.\end{eqnarray}$$
Also, the energy spectrum can be expressed by a function
$f$
as
$$\begin{eqnarray}E(k,t)=N^{2}\ell ^{3}(t)f(k\ell (t)),\end{eqnarray}$$
such that, if we expand
$f$
in power law for large scales and use (A 1), we can write
$$\begin{eqnarray}\text{for}~k\ell \ll 1,\quad E(k,t)\sim N^{2}\ell (t)^{3}(k\ell (t))^{s}\sim \text{e}^{(s+3){\it\beta}Nt/2}.\end{eqnarray}$$
Finally, the central argument to derive the growth formula is the following. As in Soulard et al. (Reference Soulard, Griffond and Gréa2014), it can be shown that the energy spectrum for
$k\ell \ll 1$
is dominated by the contribution of horizontal wavevectors, such that
${\bf\theta}={\rm\pi}/2$
. It can be also shown that the energy corresponding to these wavevectors has a linear evolution when
$1<s<4$
. Their growth rate is then given by the largest eigenvalue of the linear operator of Lin’s equations (last terms on the right-hand side of (2.24)–(2.27)) which is
$2N$
. Thus,
$$\begin{eqnarray}\text{for}~k\ell \ll 1,\quad E(k,t)\sim \text{e}^{2Nt}.\end{eqnarray}$$
Comparing (A 3) and (A 4), we obtain directly the value for
${\it\beta}$
expressed in (1.2).
Note that for
$s>4$
and
${\bf\theta}\neq {\rm\pi}/2$
, the assumption of dominance of linear terms at large scales is wrong and backscatter effects must be taken into account.
