2.1 Equations for the mixing zone (MZ)

We consider an incompressible mixing layer between two miscible fluids. We assume the density of the mixture depends linearly on the mass concentration of heavy fluid $C(\boldsymbol{x},t)\in [0~1]$ with $\boldsymbol{x}$ the position and $t$ the time. This system is governed by the classical Boussinesq/Navier–Stokes equations for $C(\boldsymbol{x},t)$ and the velocity $\boldsymbol{U}(\boldsymbol{x},t)$ :

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}C+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})C={\mathcal{D}}\unicode[STIX]{x0394}C, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{U}+(\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{U}=-\unicode[STIX]{x1D735}P+2{\mathcal{A}}CG(t)\boldsymbol{n}_{3}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\boldsymbol{U}, & \displaystyle\end{eqnarray}$$

(2.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{U}=0, & \displaystyle\end{eqnarray}$$

$P(\boldsymbol{x},t)$

$\unicode[STIX]{x1D707}$

$\boldsymbol{n}_{3}$

${\mathcal{D}}$

The position $\boldsymbol{x}$ is expressed in the usual Cartesian frame with the vertical direction corresponding to the acceleration vector identified by the suffix $_{3}$ . In this problem, quantities are assumed statistically invariant in the horizontal plane (homogeneity) as we do not consider boundary conditions. A quantity, say for example $Q(\boldsymbol{x},t)$ , is averaged along the horizontal as follows:

(2.2) $$\begin{eqnarray}\displaystyle \langle Q\rangle _{H}(x_{3},t)=\lim _{\ell \rightarrow +\infty }\frac{1}{\ell ^{2}}\int _{-\ell /2}^{+\ell /2}\int _{-\ell /2}^{+\ell /2}Q(\boldsymbol{x},t)\,\text{d}x_{1}\,\text{d}x_{2}. & & \displaystyle\end{eqnarray}$$

We introduce the Reynolds decomposition separating the quantity $Q$ between the mean $\langle Q\rangle _{H}$ and its fluctuating part $q$ as

(2.3) $$\begin{eqnarray}\displaystyle Q=\langle Q\rangle _{H}+q, & & \displaystyle\end{eqnarray}$$

with $\langle q\rangle _{H}=0$ .

In this context, the equations for the fluctuating velocity $\boldsymbol{u}(\boldsymbol{x},t)$ and concentration $c(\boldsymbol{x},t)$ deduced from (2.1a–c ) are

(2.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}c+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})c=-\unicode[STIX]{x2202}_{3}\langle C\rangle _{H}u_{3}+{\mathcal{D}}\unicode[STIX]{x0394}c, & \displaystyle\end{eqnarray}$$

(2.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+2{\mathcal{A}}cG(t)\boldsymbol{n}_{3}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

with $p(\boldsymbol{x},t)$ the fluctuating reduced pressure. Note that the mean velocity is zero $\langle \boldsymbol{U}\rangle _{H}=0$ and as also by definition fluctuating quantities $\langle \boldsymbol{u}\rangle _{H}=\langle c\rangle _{H}=\langle p\rangle _{H}=0$ . In addition, the dynamical equation for the mean concentration profile $\langle C\rangle _{H}$ , neglecting molecular diffusion, is provided by

(2.4d ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\langle C\rangle _{H}=-\unicode[STIX]{x2202}_{3}\langle u_{3}c\rangle _{H}. & & \displaystyle\end{eqnarray}$$

The closed system of equations (2.4a–d ) determines the dynamics of a mixing zone driven by an acceleration. Note that for constant positive acceleration $G_{0}>0$ , the stratification is stable if $\unicode[STIX]{x2202}_{3}\langle C\rangle _{H}>0$ and unstable for $\unicode[STIX]{x2202}_{3}\langle C\rangle _{H}<0$ .

We define the mixing zone width classically from the mean concentration profile as in Andrews & Spalding (Reference Andrews and Spalding1990):

(2.5) $$\begin{eqnarray}\displaystyle L=6\int _{-\infty }^{+\infty }\langle C\rangle _{H}(1-\langle C\rangle _{H})\,\text{d}x_{3}. & & \displaystyle\end{eqnarray}$$

The pre-factor in (2.5) is chosen such that it gives the exactly the width for a linear density (concentration) profile inside the mixing zone.

If we assume that the mean concentration remains linear inside the mixing zone and uniform equal to $0$ or $1$ outside, an equation for the mixing zone width can be derived from (2.4d ), (2.5) (see also Gréa Reference Gréa2013) leading to

(2.6) $$\begin{eqnarray}\displaystyle \dot{L}=\frac{12}{L}\int _{-\infty }^{+\infty }\langle u_{3}c\rangle _{H}\,\text{d}x_{3}. & & \displaystyle\end{eqnarray}$$

2.2 The stratified homogeneous turbulence (SHT) model

In order to analyse the mechanisms controlling buoyancy-driven turbulence, we consider a simplified problem in place of system equations (2.4a–d ). To this purpose, varied idealized models have been previously introduced relying on homogeneity assumptions for the turbulent quantities (see for instance Batchelor, Canuto & Chasnov (Reference Batchelor, Canuto and Chasnov1992), Livescu & Ristorcelli (Reference Livescu and Ristorcelli2007), Chung & Pullin (Reference Chung and Pullin2010) or Burlot et al. (Reference Burlot, Gréa, Godeferd, Cambon and Griffond2015)).

The latter one, addressed to stratified homogeneous turbulence (SHT) model, considers a uniform mean concentration gradient ( $\unicode[STIX]{x2202}_{3}\langle C\rangle _{H}=-1/L$ ) leading to

(2.7a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}c+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})c=\frac{1}{L}u_{3}+{\mathcal{D}}\unicode[STIX]{x0394}c, & \displaystyle\end{eqnarray}$$

(2.7b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+2{\mathcal{A}}cG(t)\boldsymbol{n}_{3}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.7c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where the ensemble averages satisfy $\langle \boldsymbol{u}\rangle =\langle c\rangle =0$ . In this flow, averaged quantities do not depend on the position vector $\boldsymbol{x}$ due to homogeneity.

In order to complete the system and to mimic the evolution of the mean density gradient, an equation for the mixing zone width must be added. A possible choice, directly inspired by (2.6), is

(2.7d ) $$\begin{eqnarray}\displaystyle \dot{L}=\frac{12}{L}\ell _{0}\langle u_{3}c\rangle , & & \displaystyle\end{eqnarray}$$

where $\ell _{0}$ is an arbitrary characteristic length introduced for dimensional purposes such that $\ell _{0}\langle u_{3}c\rangle$ represents the concentration flux integrated along the mixing zone. Other choices are possible, such as for instance $\ell _{0}\sim L$ (Griffond, Gréa & Soulard Reference Griffond, Gréa and Soulard2014). The difference between the two versions relies in the interpretation of results with respect to the inhomogeneous problem. For instance, the choice in this work is $\ell _{0}=1$ in (2.7d ), so that the second-order correlations scale as integrated quantities along the mixing zone while with $\ell _{0}\sim L$ they scale as turbulent quantities in the middle of the mixing zone.

It should be stated clearly that the SHT model, although sharing many properties with the complete inhomogeneous problem, cannot be viewed as an approximation of it. Note that the decoupling between turbulent quantities and the mixing width in SHT is advantageous for simulations as will be seen in § 4.1.

2.3 The rapid acceleration (RA) model

In the context of the SHT framework, we can define the various turbulent spectra as Fourier transforms of two-point correlations: see for instance Burlot et al. (Reference Burlot, Gréa, Godeferd, Cambon and Griffond2015). Due to axisymmetry, these spectra depend only on the wave modulus $k$ and the angle $\unicode[STIX]{x1D703}$ between the vertical $\boldsymbol{n}_{3}$ and the wave vector $\boldsymbol{k}$ . We can further introduce $k$ -integrated spectra for the vertical velocity ${\mathcal{E}}_{u_{3}u_{3}}$ , the concentration ${\mathcal{E}}_{cc}$ and the concentration flux ${\mathcal{E}}_{u_{3}c}$ which depend on $\unicode[STIX]{x1D703}$ only and are related to one-point correlations as follows:

(2.8a-c ) $$\begin{eqnarray}\displaystyle \langle u_{3}u_{3}\rangle =\int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{E}}_{u_{3}u_{3}}\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703},\quad \langle cc\rangle =\int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{E}}_{cc}\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703},\quad \langle u_{3}c\rangle =\int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{E}}_{u_{3}c}\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

From (2.7a–c ) and following the same procedure as in Gréa (Reference Gréa2013), it is possible to express simply the dynamics of integrated spectra by neglecting nonlinear and viscous/diffusive terms, retaining only the effects of buoyancy and pressure:

(2.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}{\mathcal{E}}_{u_{3}c}=\frac{1}{L}{\mathcal{E}}_{u_{3}u_{3}}-2{\mathcal{A}}G(t)\sin ^{2}(\unicode[STIX]{x1D703}){\mathcal{E}}_{cc}, & \displaystyle\end{eqnarray}$$

(2.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}{\mathcal{E}}_{u_{3}u_{3}}=-4{\mathcal{A}}G(t)\sin ^{2}(\unicode[STIX]{x1D703}){\mathcal{E}}_{u_{3}c}, & \displaystyle\end{eqnarray}$$

(2.9c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}{\mathcal{E}}_{cc}=\frac{2}{L}{\mathcal{E}}_{u_{3}c}, & \displaystyle\end{eqnarray}$$

(2.9d ) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{L}=\frac{12}{L}\int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{E}}_{u_{3}c}\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

The nonlinearity remaining in the system of equations (2.9a–d ) is due to the evolution of the mixing zone width $L$ . Otherwise the system describes classically the dynamics of gravity waves. In practice, this approach is valid only when the acceleration is strong compared to the inertia of turbulent structures and lasts only for a short period of time before nonlinear terms coupling the turbulent fluctuation modes recover their intensity. This statistical approach for integrated spectra described by (2.9a–d ) can be referred to as the rapid acceleration (RA) model. It has also been introduced in the context of inhomogeneous mixing zones for Rayleigh–Taylor instability in Gréa (Reference Gréa2013).

At this stage, we make an assumption in order to reduce the RA model for the evolving mixing zone into a single differential equation. If the fluctuating concentration $c$ and vertical velocity $u_{3}$ are initially correlated, implying ${\mathcal{E}}_{u_{3}c}^{2}={\mathcal{E}}_{u_{3}u_{3}}{\mathcal{E}}_{cc}$ , it can be easily shown from the RA equations that it remains correlated for all $t$ . In this case, it is possible to simplify the system by introducing ${\mathcal{U}}(t,\unicode[STIX]{x1D703})={\mathcal{E}}_{u_{3}u_{3}}^{1/2}(t,\unicode[STIX]{x1D703})$ and ${\mathcal{B}}(t,\unicode[STIX]{x1D703})={\mathcal{E}}_{cc}^{1/2}(t,\unicode[STIX]{x1D703})$ , leading to

(2.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}{\mathcal{U}}=-2{\mathcal{A}}G(t)\sin ^{2}(\unicode[STIX]{x1D703}){\mathcal{B}}, & \displaystyle\end{eqnarray}$$

(2.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}{\mathcal{B}}=\frac{1}{L}{\mathcal{U}}, & \displaystyle\end{eqnarray}$$

(2.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{L}=\frac{12}{L}\int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{U}}{\mathcal{B}}\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

Alternatively, it may be noted that the system can be expressed in the form of a nonlinear second-order equation for ${\mathcal{B}}$ , as the equation for $L$ can be integrated:

(2.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{tt}{\mathcal{B}}(t,\unicode[STIX]{x1D703})+\frac{\dot{L}(t)}{L(t)}\unicode[STIX]{x2202}_{t}{\mathcal{B}}(t,\unicode[STIX]{x1D703})+\frac{2{\mathcal{A}}G(t)\sin ^{2}(\unicode[STIX]{x1D703})}{L(t)}{\mathcal{B}}(t,\unicode[STIX]{x1D703})=0, & \displaystyle\end{eqnarray}$$

(2.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle L(t)=6\int _{0}^{\unicode[STIX]{x03C0}}{\mathcal{B}}^{2}(t,\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}+L_{0}, & \displaystyle\end{eqnarray}$$

$L_{0}$

$B(t^{\prime },\unicode[STIX]{x1D703})={\mathcal{B}}(t,\unicode[STIX]{x1D703})/(6L_{0})^{1/2}$

$t^{\prime }=N_{0}t$

$N_{0}=(2{\mathcal{A}}G_{0}/L_{0})^{1/2}$

$t$

$t^{\prime }$

(2.12a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{tt}B(t,\unicode[STIX]{x1D703})-2\frac{\dot{\unicode[STIX]{x1D6FA}}(t)}{\unicode[STIX]{x1D6FA}(t)}\unicode[STIX]{x2202}_{t}B(t,\unicode[STIX]{x1D703})+\sin ^{2}(\unicode[STIX]{x1D703})\unicode[STIX]{x1D6FA}^{2}(t)(1+F\cos (\unicode[STIX]{x1D714}/N_{0}t))B(t,\unicode[STIX]{x1D703})=0, & \displaystyle\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}^{2}(t)=\frac{1}{1+\displaystyle \int _{0}^{\unicode[STIX]{x03C0}}B^{2}(t,\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}}. & \displaystyle\end{eqnarray}$$

3.1 Linear analysis

For a specified initial mixing zone width $L_{0}$ , or equivalently an initial frequency $N_{0}$ , the RA model has a trivial solution corresponding to $B(t,\unicode[STIX]{x1D703})=\unicode[STIX]{x2202}_{t}B(t,\unicode[STIX]{x1D703})=0$ . The linear equation for the fluctuations around this state can be simply obtained by replacing $\unicode[STIX]{x1D6FA}=1$ and $\dot{\unicode[STIX]{x1D6FA}}=0$ in (2.12a ). This gives an infinite set of decoupled Mathieu oscillators parametrized by the orientation angle $\unicode[STIX]{x1D703}\in [0~\unicode[STIX]{x03C0}/2]$ , which are further characterized by the continuous frequency $\sin (\unicode[STIX]{x1D703})$ .

The instability zones for each oscillators are associated with the classical resonance conditions $n\unicode[STIX]{x1D714}/2=N_{0}\sin \unicode[STIX]{x1D703}$ ( $\forall n$ integer). For instance, a sub-harmonic instability is triggered for $n=1$ while $n=2$ corresponds to the harmonic case. The instability regions in the $((N_{0}/\unicode[STIX]{x1D714})^{2},F)$ -plane, the Mathieu stability diagram, are the well-known instability tongues and their boundaries the transition curves as shown in figure 1. In practice, the instability zones are obtained numerically following the method described in Kumar & Tuckerman (Reference Kumar and Tuckerman1994). The linear problem is always stable if the entire segment $[0,(N_{0}/\unicode[STIX]{x1D714})^{2}]$ lies inside the stable region. Otherwise, if at least one angle $\unicode[STIX]{x1D703}$ or $(N_{0}\sin \unicode[STIX]{x1D703}/\unicode[STIX]{x1D714}^{2})^{2}$ lies inside the unstable region, then linear analysis shows there would be a growth of $B$ and hence the mixing zone.

Figure 1. Stability diagram of the Mathieu equation as a function of ratio of frequency $N^{2}/\unicode[STIX]{x1D714}^{2}$ and amplitude of forcing $F$ . Harmonic and sub-harmonic unstable regions are shown. Green solid line: frequencies excited in the RA model corresponding to the different orientation angles $\unicode[STIX]{x1D703}$ of gravity waves. The two graphs show the situation (a) at an initial time and (b) at a later time after the growth and saturation of $L$ due to the parametric instability. Green star: modes corresponding to the maximum growth rate of the instability. Red dashed line: neutral curve corresponding to the expected saturation width of the mixing zone.

From this, one can easily deduce the following sufficient condition for instability:

(3.1) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x1D714}}{N_{0}}\right)^{2}<{\mathcal{G}}(F), & & \displaystyle\end{eqnarray}$$

where ${\mathcal{G}}$ indicates the first resonance boundary region triggering parametric growth of $B$ . This suggests the following scenario: as long as there is instability for some angles $\unicode[STIX]{x1D703}$ , the mixing zone will enlarge, reducing the buoyancy frequency $N(t)$ . We can expect saturation when condition (3.1) is no longer fulfilled, i.e.

(3.2a,b ) $$\begin{eqnarray}\displaystyle \left(\frac{\unicode[STIX]{x1D714}}{N_{sat}}\right)^{2}={\mathcal{G}}(F)\quad \text{or}\quad L_{sat}=\frac{2{\mathcal{A}}G_{0}}{\unicode[STIX]{x1D714}^{2}}{\mathcal{G}}(F). & & \displaystyle\end{eqnarray}$$

The dynamics of the instability is also governed by angles corresponding to maximum growth rates, not necessarily $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ as indicated in figure 1. It is interesting to note that if $(N/\unicode[STIX]{x1D714})^{2}$ is large enough, then both harmonic and sub-harmonic instabilities can be triggered at the same time but for different angles (see § 4.2). If initially $(\unicode[STIX]{x1D714}/N_{0})^{2}\geqslant {\mathcal{G}}(F)$ , the instability is not triggered.

Let us call (3.2) the saturation criterion, and the main objective of this work is to verify its validity. Indeed nonlinear effects corresponding to the evolution of $L$ need also to be considered in the analysis.

3.2 Perturbation analysis

The aim of this section is to confirm the validity of the saturation criterion (3.2) even by taking into account weakly nonlinear effects. We use a classical multiple time scale analysis (see Bender & Orszag Reference Bender and Orszag1978, Godrèche & Manneville Reference Godrèche and Manneville2005) to find a perturbation expansion of the RA model solutions. These techniques are also used to find transition curves of the Mathieu equation under small forcing or to study weakly nonlinear oscillators.

In order to solve (2.12a,b ), we propose the following expansion for $B$ :

(3.3) $$\begin{eqnarray}\displaystyle B(t,\unicode[STIX]{x1D703})=\unicode[STIX]{x1D700}(B^{(0)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})+\unicode[STIX]{x1D700}B^{(1)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})+\unicode[STIX]{x1D700}^{2}B^{(2)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})+\unicode[STIX]{x1D700}^{3}B^{(3)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})+\cdots \,). & & \displaystyle\end{eqnarray}$$

Here, $\unicode[STIX]{x1D700}$ expresses the fact that the perturbation $B$ is small so we can perform a weakly nonlinear analysis. We introduce the slow time variable $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D700}^{2}t$ . This choice differs from multiple scale analysis of the Mathieu equation and comes from the fact that we approximate the integral over $\unicode[STIX]{x1D703}$ in (2.12b ) by a stationary phase expansion.

The multiple scale expansion of (2.12a,b ) without parametric forcing ( $F=0$ ) is not detailed here but given in appendix A to lighten this section. Let us give a short summary of the principal results of this analysis. At leading order, the nonlinear terms introduce a phase shift but do not modify the oscillation amplitude. At next order, the amplitudes for $B$ are modified, giving finite amplitude oscillations for $L$ . However, we shall see below that the leading order is sufficient to take into account parametric effects.

For the parametric forcing analysis, we assume further that the natural frequency of the system corresponding to $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ is close to the first resonance condition and that the forcing amplitude is small. This can be written as

(3.4a,b ) $$\begin{eqnarray}\displaystyle N_{0}^{2}/\unicode[STIX]{x1D714}^{2}-1/4=\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6E5}\quad \text{and}\quad F=\unicode[STIX]{x1D700}^{2}f & & \displaystyle\end{eqnarray}$$

introducing $\unicode[STIX]{x1D6E5}$ and $f$ as the detuning and the forcing amplitude parameter, respectively. In appendix A, it is shown that the nonlinear effects appear at order $\unicode[STIX]{x1D700}^{3}$ for the natural frequency $\unicode[STIX]{x1D6FA}^{2}(t)$ and at order $\unicode[STIX]{x1D700}^{4}$ for the damping term $-2\dot{\unicode[STIX]{x1D6FA}}(t)/\unicode[STIX]{x1D6FA}(t)$ (see (A 1), (A 7)).

Reintroducing the different expansions in (2.12a ) gives the homogeneous equations for the first two leading orders,

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{tt}B^{(0)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})+\sin ^{2}(\unicode[STIX]{x1D703})B^{(0)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})=0\quad \text{at order}~\unicode[STIX]{x1D700}, & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{tt}B^{(1)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})+\sin ^{2}(\unicode[STIX]{x1D703})B^{(1)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})=0\quad \text{at order}~\unicode[STIX]{x1D700}^{2}, & \displaystyle\end{eqnarray}$$

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle B^{(0)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})=a(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})\text{e}^{\text{i}\sin (\unicode[STIX]{x1D703})t}+a^{\ast }(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})\text{e}^{-\text{i}\sin (\unicode[STIX]{x1D703})t}, & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle B^{(1)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})=b(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})\text{e}^{\text{i}\sin (\unicode[STIX]{x1D703})t}+b^{\ast }(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})\text{e}^{-\text{i}\sin (\unicode[STIX]{x1D703})t}. & \displaystyle\end{eqnarray}$$

$B^{(0)}(t,\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})$

$\unicode[STIX]{x1D700}^{2}$

(3.7) $$\begin{eqnarray}\displaystyle \frac{L(t)}{L_{0}}=\frac{1}{\unicode[STIX]{x1D6FA}^{2}(t)}=1+2\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D70F})+O(\unicode[STIX]{x1D700}^{3}), & & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D70F})=\int _{0}^{+\unicode[STIX]{x03C0}}|a(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})|^{2}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}$ . In (3.7), the integral over $\unicode[STIX]{x1D703}$ of the oscillating terms $a^{2}(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})\text{e}^{2\text{i}\sin (\unicode[STIX]{x1D703})t}$ and $a^{\ast 2}(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})\text{e}^{-2\text{i}\sin (\unicode[STIX]{x1D703})t}$ does not contribute to the leading order of the expansion due to the stationary phase approximation and are postponed to order $\unicode[STIX]{x1D700}^{3}$ .

At order $\unicode[STIX]{x1D700}^{3}$ in (2.12a ), we have

(3.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{tt}B^{(2)}+\sin ^{2}(\unicode[STIX]{x1D703})B^{(2)}=-2\unicode[STIX]{x2202}_{t\unicode[STIX]{x1D70F}}B^{(0)}+2\sin ^{2}\unicode[STIX]{x1D703}B^{(0)}\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D70F})-f\sin ^{2}\unicode[STIX]{x1D703}B^{(0)}\sin (\unicode[STIX]{x1D714}/N_{0}t). & & \displaystyle\end{eqnarray}$$

Classically, the equation for the amplitudes $a(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D703})$ in multiple scale analysis are obtained by cancelling secular terms in (3.8). Let us focus on the angle $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ for the leading-order equation (3.6a ):

(3.9) $$\begin{eqnarray}\displaystyle -2\text{i}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}a(\unicode[STIX]{x1D70F})+2a(\unicode[STIX]{x1D70F})\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D70F})-{\textstyle \frac{1}{2}}fa^{\ast }(\unicode[STIX]{x1D70F})\text{e}^{-4\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}=0, & & \displaystyle\end{eqnarray}$$

where the damping still does not appear.

We pose $a(\unicode[STIX]{x1D70F})=g(\unicode[STIX]{x1D70F})\text{e}^{-2\text{i}\unicode[STIX]{x0394}\unicode[STIX]{x1D70F}}$ to render the system autonomous which leads to

(3.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}g(\unicode[STIX]{x1D70F})=\text{i}(2\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6EC}(\unicode[STIX]{x1D70F}))g(\unicode[STIX]{x1D70F})-\text{i}\frac{f}{4}g^{\ast }(\unicode[STIX]{x1D70F}). & & \displaystyle\end{eqnarray}$$

Separating modulus and phase of $g=R(\unicode[STIX]{x1D70F})\text{e}^{\text{i}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70F})}$ leads to

(3.11a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}R(\unicode[STIX]{x1D70F})=-\frac{f}{4}R\sin (2\unicode[STIX]{x1D6FC}), & \displaystyle\end{eqnarray}$$

(3.11b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D70F})=2\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6EC}+\frac{f}{4}\cos (2\unicode[STIX]{x1D6FC}). & \displaystyle\end{eqnarray}$$

$R=0$

$2\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6EC}\pm f/4=0$

$\cos (2\unicode[STIX]{x1D6FC})=1$

$R=0$

$f/4>2|\unicode[STIX]{x1D6E5}|$

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6EC}=2\unicode[STIX]{x1D6E5}+\frac{f}{4}. & & \displaystyle\end{eqnarray}$$

Equation (3.12) is what we seek. Using that $L/L_{0}=1+2\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6EC}$ and $(N_{0}^{2}/\unicode[STIX]{x1D714}^{2})-(1/4)=\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6E5}$ , we obtain easily at leading order:

(3.13) $$\begin{eqnarray}\displaystyle \frac{L_{sat}}{L_{0}}=\frac{N_{0}^{2}}{\unicode[STIX]{x1D714}^{2}}(2F+4). & & \displaystyle\end{eqnarray}$$

The relation equation (3.13) is equivalent to the criterion of equation (3.2); in the limit of small forcing the first transition curve corresponds in this case to ${\mathcal{G}}(F)=2F+4$ .

In conclusion for this part, we have proved by perturbation analysis the validity of the saturation criterion equation (3.2) under the conditions of small forcing and initial perturbation and if the natural frequency satisfies the first sub-harmonic resonance condition.

3.2.1 Numerical solutions of the RA model

In order to explore more broadly the behaviour of the RA model and to assess the validity of the saturation criterion equation (3.2), even with strong nonlinear effects, we compute the solutions of (2.12a,b ) for different $\unicode[STIX]{x1D714}/N_{0}$ ratios and forcing amplitudes $F$ . We choose at $t=0$ the values, $B(0,\unicode[STIX]{x1D703})=0.05$ and $\unicode[STIX]{x2202}_{t}B(0,\unicode[STIX]{x1D703})=0$ . In figure 2, we plot the initial conditions and saturation values in the Mathieu phase diagram. The time evolution of $L$ for a representative case with $F=1.7$ and $N_{0}^{2}/\unicode[STIX]{x1D714}^{2}=3/4$ is also presented.

Figure 2. (a) Representation of the numerical solutions of the RA model in the Mathieu stability diagram. Small symbols: initial conditions. Large symbols: final conditions. Black line: neutral curves of the Mathieu stability diagram. (b) Time evolution of the mixing zone width $L$ in the RA model corresponding to $F=1.7$ and $N_{0}^{2}/\unicode[STIX]{x1D714}^{2}=3/4$ (large circle in the Mathieu stability diagram). Dashed line: expected saturation width.

The parametric instability in the RA model is accompanied by strong oscillations on $L$ and a growth of its mean value. The instability is sub-harmonic as indeed the initial conditions are in the neighbourhood of the first instability tongue. Accordingly, the frequency of oscillation for $L(t)$ with the non-dimensional time is close to $\unicode[STIX]{x1D714}/N_{0}$ during the growth, corresponding to an oscillation period of $\unicode[STIX]{x1D714}/N_{0}/2$ for $B$ (not shown). This recovers the classical Faraday phenomenology.

Before discussing the validity of the saturation criterion, equation (3.2), we need to clarify how $L_{sat}$ is evaluated in our simulations. The RA model has no dissipation mechanisms and the final states of $L$ are strongly oscillating. These oscillations are not damped by the phase mixing between the solutions corresponding to different angles $\unicode[STIX]{x1D703}$ . They remain at late time even if $L/L_{0}\geqslant N_{0}^{2}/\unicode[STIX]{x1D714}^{2}{\mathcal{G}}(F)$ , i.e. the resonance conditions driving the parametric instability are no longer valid. This interesting effect is the imprint of nonlinearities and is observed also when $F=0$ . Using the multiple scale analysis proposed in § A.2, we can explain this phenomenon as the self-excitation of $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ modes induced by the variations of $L$ . In practice, we simply evaluate $L_{sat}$ by averaging $L$ on the interval $\unicode[STIX]{x1D714}t/N_{0}\in [150~200]$ on which the simulations seem to have reached a steady oscillating state.

In figure 2, we see that the final states $L_{sat}$ as obtained by the previous procedure are close to the saturation criterion expressed by (3.2) even if the forcing amplitude $F$ is large. All the final states are located in the stable zone of the Mathieu equation. The prediction works better with initial conditions close to the weakly nonlinear analysis case, i.e. in the neighbourhood of the first resonance condition and for small $F$ . Globally, these results bring strong support to the saturation criterion. The next step of our analysis consists then in checking its validity beyond the RA model.

4.1 Simulation description

We provide the different characteristics of SHT and mixing layer simulations dedicated to Faraday instability in the first and the second part of this section respectively.

4.1.1 SHT simulations

We use for SHT simulations a massively parallel direct numerical simulation (DNS) code based on classical numerical pseudo-spectral methods detailed in Griffond et al. (Reference Griffond, Gréa and Soulard2014). The simulation domain for the turbulent velocity and concentration field is a triply periodic box of size $2\unicode[STIX]{x03C0}$ . For convenience, we work with non-dimensional wavenumbers; in particular, $k_{min}=1$ corresponds to the box size.

Three levels of resolution – small (S) medium (M) and large (L) – with $256^{3}$ , $512^{3}$ , $1024^{3}$ grid points are performed. The characteristics of the different runs at the initial time are detailed in table 1. The viscosity and diffusion coefficients are taken equal (Schmidt number $=1$ ) with values $\unicode[STIX]{x1D707}={\mathcal{D}}=5\times 10^{-3},2\times 10^{-3}$ and $8\times 10^{-4}$ respectively for the S, M and L cases. These values ensure that the different scales of the flow are well resolved even during the instability phase (DNS).

The small and medium simulations (S and M) are principally dedicated to the parametric study, in particular to test the initial positions in the Mathieu diagram (Run A) and the influence of initial perturbation intensity (Run B). There is also a special set of M simulations with $F=0$ (Run A) to ensure that the instability is not triggered when the forcing is switched off. The purpose of large simulations is to investigate turbulent quantities (in particular spectra) during the Faraday instability (Run C) and also to confirm the convergence of results obtained by smaller ones (S and M). Note that large simulations are quite computationally intensive (consuming typically ${>}5\times 10^{5}$ CPU hours spread over $1024$ cores) in order to properly capture the instability development (around $20$ acceleration periods).

In all our SHT simulations, the mean reduced accelerations and frequencies are fixed respectively at $2{\mathcal{A}}G_{0}=1$ and $\unicode[STIX]{x1D714}=2$ . The initial mixing width $L(t=0)$ , amplitude $F$ and turbulent quantities may vary in order to explore the phase space. In total, $220$ SHT simulations are presented in this work.

Table 1. Non-dimensional parameters at $t=0$ and number $n_{s}$ of SHT simulations for the different cases.

The random initial conditions for the concentration $c$ and velocity $u_{i}$ are isotropic. As a consequence we have $\langle u_{i}c\rangle =0$ . A narrow wavenumber band centred around $k_{peak}$ is initially excited in Fourier space (peaked spectra) and we specify further the variances $\langle cc\rangle$ and $\langle u_{i}u_{i}\rangle$ . The choice for $k_{peak}$ is a compromise between the good resolution of the simulations at small scales and keeping a margin for the integral scale to develop freely during the simulations due to nonlinearities. At the end of the simulations, the integral scales usually approach the box size. This effect is difficult to avoid because the Faraday instability needs time to evolve. In Run C, we have tried to reduce these confinement effects by increasing $k_{peak}$ and the resolution. No significant differences concerning the dynamics have been observed compared to less resolved simulations. This is also consistent with results of Thornber (Reference Thornber2016) in the different context of decaying turbulence showing that the effects coming from the integral scale confinement appear very lately in spectral simulations.

The flow corresponding to our initial condition with peaked spectra is not turbulent. However, the spectra broaden quite fast due to nonlinear interactions. In addition, the energy brought to the fluctuations by the instability renders it turbulent at later time. Indeed, the turbulent Reynolds number based on kinetic energy, dissipation and viscosity grows due to the instability and reaches values up to ${\sim}2000$ in L simulations and up to ${\sim}500$ in S simulations.

An important parameter which describes our simulations is the initial Froude number, defined as $Fr_{I}=(1/2\langle u_{i}u_{i}\rangle )^{1/2}k_{peak}/N_{0}$ , giving the ratio between the perturbation and the stratification frequency (see table 1). Here we use a perturbation frequency instead of the classical turbulent one defined from kinetic energy and dissipation since our initial condition is not turbulent. The Froude number expresses whether the simulations are initially buoyancy-driven or ruled by nonlinear effects. This initial Froude number can be easily tuned in our SHT simulations and we show a parametric study of its influence in Run B (discussion postponed to § 4.2).

4.1.2 MZ simulations

The mixing zone simulations (MZ) are performed with the finite difference parallel code, TURBMIX3D, presented in Poujade & Peybernes (Reference Poujade and Peybernes2010), Watteaux (Reference Watteaux2011). This code solves the Navier–Stokes equations for a binary mixture of two incompressible fluids corresponding in the limit of small Atwood number to (2.1a–c ). Here, the Atwood number is chosen as ${\mathcal{A}}=10^{-3}$ . The computational domain is a cubic box of normalized size $L_{box}=1$ with periodicity imposed only in the horizontal direction. The mixing zone width $L$ is classically evaluated from the volume fraction profile defined in (2.5).

Simulations with small (S) $256^{3}$ and medium (M) $512^{3}$ resolutions are presented. The viscosity and the diffusivity coefficients are equal (Schmidt number $=1$ ) and kept constant but the viscous/diffusive scales are not fully resolved. For this reason, our simulations cannot claim to be DNS but are in the category of implicit large eddy simulations. We motivate this choice to reduce viscous or diffusive effects at large scales which may influence the dynamics of large turbulent structures and the mean density gradient. In particular, we try to have a simulation time $t_{sim}\ll L_{sat}^{2}/\unicode[STIX]{x1D708}$ to avoid diffusively evolving mixing layers after the Faraday phase. Note that SHT simulations are not limited by this condition.

Another constraint in MZ simulations which is not present in the SHT framework is that the mixing layer must be smaller than the size of the domain $L\leqslant L_{box}$ . In practice we try to obtain $L_{sat}\sim 0.5{-}0.6\times L_{box}$ , which seems sufficient to avoid confinement effects as shown by Morgan et al. (Reference Morgan, Olson, White and McFarland2017) in the context of Rayleigh–Taylor mixing zones. Also, the use of different resolutions leading to same results guarantees that confinement effects seem very limited in the simulations presented. This leads to the parameters given in table 2 for the $29$ MZ simulations.

Table 2. Non-dimensional parameters at $t=0$ and number $n_{s}$ of mixing zone simulations.

The initial conditions for our MZ simulations correspond to a developed mixing zone generated by a Rayleigh–Taylor (RT) instability with constant acceleration $G_{RT}$ . The RT mixing phase is initialized using a multi-mode perturbation of the interface. The spectrum for the perturbation height is annular as in Dimonte et al. (Reference Dimonte, Youngs, Dimits, Weber, Marinak, Wunsch, Garasi, Robinson, Andrews and Ramaprabhu2004) with wavenumbers $k\in [32~64]$ for S simulations, and $k\in [64~128]$ for M simulations. Once mixing has been created, the RT phase is halted relatively fast by around the time $\unicode[STIX]{x1D70F}=(L_{box}/{\mathcal{A}}G_{RT})^{1/2}$ corresponding to the nonlinear mixing transition. At this time, nonlinear coupling and merging processes do not produce fully turbulent mixing (Cook & Dimotakis Reference Cook and Dimotakis2001), but this choice appears as a compromise because of limitations in the final width of $L_{sat}$ due to the top/bottom boundaries of the box. Also due to this constraint, the Froude number $Fr_{I}$ defined with the integral scale and kinetic energy at the centre of the mixing layer are initially very small.

In summary, we principally dedicate MZ simulations to challenging the validity of the saturation criterion in realistic configurations accounting for inhomogeneous effects neglected by the theory. In complement, the more tractable SHT simulations allow a parametric study of the influence of initial conditions and to investigate the scales of turbulence in more detail.