4.1. Phenomenology of intermittency at large and low $\mathit{Re}$

In figure 6, the streamwise velocity in a horizontal plane close to the upper wall, $y=0.95$ , is shown for three weakly intermittent simulations. The value of $\mathit{Re}_{{\it\tau}}$ increases from (a) to (c), corresponding to $\mathit{Re}_{{\it\tau}}=130,180$ and 520, simulations 8, 10 and 13 respectively. In all three cases, laminar patches appear in the flow, similar to those observed by García-Villalba et al. (Reference García-Villalba, Azagra and Uhlmann2011) in stratified channel flows. It should be noted that in the three cases the dimension of the box varies considerably and the laminar patches at low $\mathit{Re}_{{\it\tau}}$ are significantly larger. Nevertheless, the three flows exhibit similar dynamics with laminar patches cyclically appearing and disappearing in the flow with time scales of the order of hundreds of convective time units. For these simulations, spatially averaged quantities show small fluctuations in time, suggesting that the chosen box sizes are suitable to study the intermittent dynamics. The degree of stratification in these cases is not strong enough to lead to the generation of sustained and steady inclined laminar/turbulent bands as observed by Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012) and Deusebio et al. (Reference Deusebio, Brethouwer, Schlatter and Lindborg2014). As $\mathit{Ri}$ increases, the extent of the laminar regions increases and the dynamics becomes increasingly more intermittent in time. Figure 7(a) shows a snapshot of the streamwise velocity close to the upper wall, at $y=0.95$ , for simulation 14. There is only one band of turbulence spanning the entire streamwise extent, and this band is surrounded by laminar flow. As is apparent from the typical values of the streamwise velocity (the same colour bar is used in figures 6 c and 7 a), the streamwise velocity is substantially closer to its laminar value at $y=0.95$ for simulation 14 than for the more vigorously turbulent simulation 13. In this simulation, $\mathit{Re}_{{\it\tau}}$ varies considerably in time, ranging between 380 and 540, which suggests that the size of the computational domain is only marginally large enough to accommodate sustained intermittent dynamics. The topology found in figure 7 was similarly observed in simulations 15 and 16, confirming that also in these cases the computational domains were close to marginal.

Figure 6.

Streamwise velocity in a horizontal plane close to the upper wall, at $y=0.95$ . It should be noted that the full domain is shown in each of the three cases, which vary in size. (a) Simulation 8; (b) simulation 10; (c) simulation 13.

Figure 7.

Flow field for a highly intermittent simulation (simulation 14). The entire domain is shown in both panels. (a) Streamwise velocity in a horizontal plane close to the upper wall, at $y=0.95$ . The full domain is shown and colours range from blue (low values) to red (high values), using the same colour bar as figure 6(c). (b) Three-dimensional view of the enstrophy in three planes normal to the Cartesian axes: $x$ -normal at $x=12$ (black), $y$ -normal at $y=-0.2$ (red) and $z$ -normal at $z=1$ (blue). The scale logarithmically ranges from 1 $U^{2}/h^{2}$ to 200 $U^{2}/h^{2}$ and it is the same in all of the plots. Colours are used only to highlight which plane the contours belong to.

We only observe bands similar to those observed by Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012) at low $\mathit{Re}_{{\it\tau}}<100$ , for which very large computational domains could be used. In smaller computational domains, the spatial intermittency is replaced by temporal intermittency and the turbulent bands align in the streamwise direction. However, the fact that at $\mathit{Re}_{{\it\tau}}=540$ we find streamwise stripes of turbulence similar to those found in low- $\mathit{Re}$ small-domain simulations suggests that inclined bands might be observed at higher $\mathit{Re}$ if larger computational domains were employed.

At low $\mathit{Re}$ , the laminar regions extend over the entire channel depth. However, we find that the vertical structure of the laminar/turbulent regions is inherently different for simulations at low and high $\mathit{Re}$ . As $\mathit{Ri}$ increases, laminar regions appearing at higher $\mathit{Re}$ become more confined to the walls. Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012) found that in stratified channel flows partial relaminarisation occurs at the walls, whereas the dynamics in the interior remains highly turbulent, without any dominant presence of large-scale structures. In figure 7(b) the enstrophy (a useful and simple proxy for turbulence) in planes normal to the three Cartesian axis is shown for simulation 14. The contours projected onto the three Cartesian planes are coloured differently (although with the same intensity scale). The $z$ -normal plane clearly shows the presence of two turbulent layers separated by a quieter laminar region close to the channel mid-plane. This quieter region is confined horizontally and it is bordered in the spanwise direction by a turbulent patch which originates from the bending of the upper vertically confined turbulent layer, as shown by the $x$ -normal plane. Figure 7(b) highlights a rich dynamical picture in which intermittency at the wall connects to intermittent dynamics in the centre of the channel, with streamwise turbulent layers bending and mutually interacting. The dynamical transition between the turbulent bands (filling the entire channel gap) observed at low $\mathit{Re}$ and the layering structure in figure 7(b) is gradual. The bands moderately develop a spanwise inclination which leads to an inhomogeneity in the vertical direction.

It is unclear whether the absence of intermittent dynamics in the interior as observed by Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012) is due to a weaker stratification or to a difference of flow configuration with respect to our simulations. It is nevertheless worth noting that in PCF intermittency at the wall significantly affects the overall dynamics. From an energetic point of view, PCF is forced by shear stresses at the upper and lower walls, and the dynamics in the interior adjusts in order to carry the momentum flux resulting from the near-wall dynamics. Full or partial relaminarisation leads to a significant drop in the magnitude of the wall shear stress ${\it\tau}_{w}$ , thus also strongly affecting the dynamics in the interior. The influence of the partial relaminarisation at the walls might therefore be significantly stronger than in other flow configurations where the driving force is not affected by relaminarisation, such as the channel flows considered by Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012) driven by a constant pressure gradient or the throttling method used by Chung & Matheou (Reference Chung and Matheou2012) in stably stratified homogeneous turbulence.

4.2. Quantifying intermittency

One of our main aims here is to describe the intermittent regime arising due to a stable stratification. In order to quantify the degree of intermittency, an objective metric based on the identification of laminar and turbulent regions is needed. Several choices are possible, including evaluation of shear stress, vertical velocity (Brethouwer et al. Reference Brethouwer, Duguet and Schlatter2012) and vorticity (Corrsin & Kistler Reference Corrsin and Kistler1955; Pope Reference Pope2000).

We base our criterion on the wall enstrophy,

(4.1) $$\begin{eqnarray}{\it\eta}_{\pm }(x,z,t)=\left(\frac{\partial u}{\partial y}\right)^{2}+\left(\frac{\partial w}{\partial y}\right)^{2}\quad \text{at }y=\pm 1,\end{eqnarray}$$

and identify the turbulent fraction ${\it\gamma}$ using moments of the distribution of ${\it\eta}_{\pm }$ as discussed below. We are motivated to use the wall enstrophy by two considerations. First, as observed by Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012), at large $\mathit{Re}$ when intermittent dynamics appears laminar patches are generally confined close to the walls. Second, as discussed in § 4.1, because the driving force is applied at the solid walls, near-wall dynamics plays an important role in maintaining turbulence in PCF. In order to identify ${\it\gamma}$ , we consider the horizontal spatial fluctuation of ${\it\eta}_{\pm }$ at the wall,

(4.2) $$\begin{eqnarray}{\it\eta}_{\pm ,mbf}=\sqrt{\frac{1}{A}\int _{x^{+},z^{+}}{\it\eta}_{\pm }^{2}\,\text{d}x^{+}\,\text{d}z^{+}-\left[\frac{1}{A}\int _{x^{+},z^{+}}{\it\eta}_{\pm }\,\text{d}x^{+}\,\text{d}z^{+}\right]^{2}},\end{eqnarray}$$

over subdomains corresponding to a minimal flow unit (Jiménez & Moin Reference Jiménez and Moin1991), specifically ${\rm\Delta}x^{+}\approx 200$ and ${\rm\Delta}z^{+}\approx 100$ , with $A={\rm\Delta}x^{+}\,{\rm\Delta}z^{+}$ being the subdomain area. The subscript $mbf$ in (4.2) stands for minimal box fluctuation. Although here we use ${\it\eta}_{\pm }$ as a proxy in (4.2), alternative quantities could also be used, such as velocity, temperature or turbulent stresses. When the flow is laminar, the wall enstrophy ${\it\eta}_{\pm }$ is nearly homogeneous within the subdomains and small values of (4.2) are to be expected. On the other hand, in turbulent regions the streaky patterns associated with the near-wall dynamics lead to significantly larger fluctuations. This measure is similar to the measure used by Deusebio et al. (Reference Deusebio, Brethouwer, Schlatter and Lindborg2014), although here there is a greater focus on identifying an appropriate threshold between laminar and turbulent regions. Figure 8(a) shows some examples of the probability density function (p.d.f.) of (4.2). The cases are characterised by different levels of intermittency: fully developed turbulence has a p.d.f. resembling a log-normal shape (the $x$ -axis has a logarithmic scale) with a single well-identified peak (Örlü & Schlatter Reference Örlü and Schlatter2011). As intermittent dynamics occurs in the flow, the p.d.f. skews towards the left part of the plot and, in sufficiently strongly intermittent cases, a double peak appears. The double-peak shape also allows us to identify a threshold $\mathscr{T}$ separating turbulent and laminar regions as the local minimum between the two peaks of the p.d.f. We identify laminar flow at the upper or lower wall as the region for which ${\it\eta}_{\pm ,mbf}<\mathscr{T}$ and turbulent flow as the region for which ${\it\eta}_{\pm ,mbf}>\mathscr{T}$ . In general, the threshold $\mathscr{T}$ may depend on $\mathit{Re}_{{\it\tau}}$ and we have determined its $\mathit{Re}_{{\it\tau}}$ dependence (shown in figure 8 b) by considering the most intermittent simulations at several values of $\mathit{Re}_{{\it\tau}}$ . Although the value of the threshold increases with $\mathit{Re}_{{\it\tau}}$ , it does not vary as strongly as the range spanned by the p.d.f., which in intermittent cases extends over several orders of magnitude. In the following we will consider the intermittency factor ${\it\gamma}$ for all of the simulations in table 1 by taking the average of the fraction of the upper and lower wall area for which ${\it\eta}_{\pm }>\mathscr{T}$ respectively. Typically, these two values are close. The algorithm described above was applied to single flow fields and results were averaged over at least 30 different sample times approximately equally spaced during the simulation.

Figure 8.

(a) Examples of p.d.f.s of small-box fluctuations of ${\it\eta}_{\pm }$ for a typical fully turbulent simulation ( $\cdots \cdots$ , simulation 9), weakly intermittent simulation (- - - -, simulation 10) and strongly intermittent simulation (——, simulation E). (b) Variation of $\mathscr{T}$ with $\mathit{Re}_{{\it\tau}}$ .

Although we mostly focus on the near-wall value of ${\it\gamma}$ , in order to study the mixing properties of laminar and turbulent regions separately (as discussed in § 5), we have further extended the algorithm to horizontal planes away from the wall using the total enstrophy in place of its wall value ${\it\eta}_{\pm }$ . To account for the vertical inhomogeneity, we have adjusted the value of the threshold $\mathscr{T}$ based on the local turbulent kinetic energy dissipation, ${\it\varepsilon}_{K}=\mathit{Re}^{-1}\langle \overline{\partial u_{i}^{\prime }/\partial x_{j}\partial u_{i}^{\prime }/\partial x_{j}}\rangle$ , as

(4.3) $$\begin{eqnarray}\mathscr{T}_{s}(y)=\frac{\mathscr{T}}{{\it\varepsilon}_{K,\pm 1}}{\it\varepsilon}_{K}(y),\end{eqnarray}$$

where the $\pm 1$ subscript denotes upper and lower wall values respectively. Using the threshold $\mathscr{T}_{s}$ , we define the vertically varying intermittency factor ${\it\gamma}_{s}$ as the ratio of the turbulent area to the total area for each horizontal plane. Figure 9 shows an example of the identification of laminar/turbulent regions for the flow field shown in figure 6(c). Laminar and turbulent boundaries are well estimated, both close to the wall, as shown in figure 9(a) and in the centre of the channel, as shown in figure 9(b). The average value of ${\it\gamma}_{s}$ using roughly 30 flow fields (shown in figure 9 c) slightly varies in the vertical direction, with somewhat smaller values close to the wall, where an overall intermittency factor of 82 % is found. This is consistent with the observation of Brethouwer et al. (Reference Brethouwer, Duguet and Schlatter2012), who found laminar patches to be mostly confined at the wall. That the smallest value of ${\it\gamma}_{s}$ generally corresponds to the wall value ${\it\gamma}$ also supports our conjecture that the relaminarisation process in stratified PCF originates in the near-wall region.

Figure 9.

Identification of laminar and turbulent regions using the automatic algorithm outlined in § 4.2 for a snapshot of simulation 13 at (a) $y=0.95$ and (b) $y=0$ . The thick line represents the boundary between laminar and turbulent regions which has been overlaid on the enstrophy field for reference. The flow corresponds to the same time instant as the flow depicted in figure 6(c) (simulation 13). (c) Wall-normal evolution of the intermittency factor ${\it\gamma}_{s}$ (averaged over roughly 30 snapshots).

4.3. Scaling of intermittent dynamics

As discussed in the introduction, the most appropriate indicator for the onset of intermittent dynamics is still an open issue and several criteria have been proposed. The general disagreement on the relevant quantity generally reflects the fact that it is not yet clear which physical process governs relaminarisation and how to model it. The Miles–Howard criterion (Howard Reference Howard1961; Miles Reference Miles1961), based on linear stability analysis of a stratified shear layer, is sometimes invoked to explain the suppression of turbulence, although there is now a growing consensus that turbulence may still be observed even when $\mathit{Ri}_{g}>1/4$ (see, for example, Zilitinkevich et al. Reference Zilitinkevich, Elperin, Kleeorin, Rogachevskii, Esau, Mauritsen and Miles2008). As shown in figure 3(b), $\mathit{Ri}_{g}$ attains small values at the wall and increases towards the centreline, the location at which higher values are generally found. The maximum values of $\mathit{Ri}_{g}$ before the onset of intermittency and relaminarisation show an approximate dependence with bulk $\mathit{Ri}$ , and no single value describes the appearance of intermittency. Besides, intermittency is observed at the wall where $\mathit{Ri}_{g}$ is small, $\mathit{Ri}_{g}<1/4$ .

García-Villalba & del Álamo (Reference García-Villalba and del Álamo2011) used the Nusselt number $\mathit{Nu}$ , as defined in (2.10), as a criterion to identify the appearance of intermittency in stratified channel flows, arguing that relaminarisation occurs as $\mathit{Nu}$ approaches unity. This criterion is tested in figure 10, where ${\it\gamma}$ is plotted against $\mathit{Nu}-1$ . Although as intermittency increases the data generally tend to lower values of $\mathit{Nu}$ , a general scaling cannot be observed and the data show a large scatter with $\mathit{Re}_{{\it\tau}}$ .

Figure 10.

A plot of ${\it\gamma}$ versus $\mathit{Nu}-1$ for the entire set of simulations listed in table 1. The colour code represents $\mathit{Re}_{{\it\tau}}$ , ranging from 40 (blue) to 1500 (red).

As noted in the introduction, another criterion was proposed by Flores & Riley (Reference Flores and Riley2010), who hypothesised that relaminarisation is observed when $L^{+}\approx 100$ , thus suggesting that relaminarisation is a near-wall process since $L^{+}$ is a criterion only based on inner-scaling quantities. We test this hypothesis in figure 11, where the turbulent fraction ${\it\gamma}$ is plotted versus $L^{+}$ for all of the simulations listed in table 1. We note that, apart from the two lowest values of $\mathit{Re}_{{\it\tau}}$ (plotted with the four darkest blue points in the plot), the data collapse onto a single curve, unlike the data shown in figure 10. This is quite remarkable given the $\mathit{Re}_{{\it\tau}}$ span considered in the study, which covers more than one order of magnitude. The failure to collapse the two lowest values of $\mathit{Re}_{{\it\tau}}$ is also not surprising as the $L^{+}$ criterion is expected to apply only when the onset of intermittency is determined by the imposed stratification as opposed to by viscous effects. For smaller $\mathit{Re}$ values, viscosity becomes crucial and leads to a viscosity-driven intermittency, which in the unstratified limit occurs when $\mathit{Re}\approx 325$ , corresponding to $\mathit{Re}_{{\it\tau}}\approx 35$ . Nevertheless, the asymptotic regime in which intermittency dynamics is buoyancy-driven and the intermittency factor ${\it\gamma}$ becomes independent of $Re$ is already reached at a value of $\mathit{Re}_{{\it\tau}}$ that is only three times larger than the unstratified value of 35. As shown in figure 11, in the asymptotic regime of high $\mathit{Re}$ , $\mathit{Re}>4000$ , intermittent dynamics appears at $L^{+}\approx 200$ , and the intermittency factor ${\it\gamma}$ quickly drops for smaller values.

Figure 11.

A plot of ${\it\gamma}$ versus $L^{+}$ for the entire set of simulations. The colour code represents $\mathit{Re}_{{\it\tau}}$ , ranging from 40 (blue) to 1500 (red).

6.1. Monin–Obukhov self-similarity theory

The analytical model is derived from the self-similarity theory proposed by Obukhov (Reference Obukhov1971) and Monin & Obukhov (Reference Monin and Obukhov1954). It provides analytical expressions for velocity and temperature profiles in unstable and stable stratifications. More generally, turbulent statistical moments are assumed to depend only on the vertical momentum transfer $-\overline{u^{\prime }v^{\prime }}$ (at first approximation equal to the wall shear stress ${\it\tau}_{w}$ ), the buoyancy parameter $g{\it\alpha}_{V}$ , the height from the surface $y$ and the turbulent heat flux $-\overline{v^{\prime }T^{\prime }}$ (at first approximation equal to the wall heat flux $q_{w}$ ). Dimensional analysis therefore leads to only one independent dimensionless group, suggesting that any statistical quantity, when properly normalised, is only a function of such a group, which could be (arbitrarily) identified by

(6.2) $$\begin{eqnarray}{\bf\xi}=\frac{y}{L},\end{eqnarray}$$

where $L$ is the Monin–Obukhov length scale, as defined in (2.11). In particular, the mean velocity and temperature gradients (dimensional) can therefore be written as

(6.3a,b ) $$\begin{eqnarray}\frac{\partial u}{\partial y}=\frac{u_{{\it\tau}}}{C_{{\it\kappa}}y}{\it\phi}_{m}({\bf\xi})\quad \text{and}\quad \frac{\partial T}{\partial y}=\frac{T_{{\it\tau}}}{C_{{\it\kappa}}y}{\it\phi}_{h}({\bf\xi}),\end{eqnarray}$$

where ${\it\phi}_{m}$ and ${\it\phi}_{h}$ are dimensionless functions, the von Kármán coefficient, $C_{{\it\kappa}}$ , is historically included in the definitions and $T_{{\it\tau}}$ is the friction temperature, i.e. $T_{{\it\tau}}=q_{w}/u_{{\it\tau}}$ . At first approximation, ${\it\phi}_{m}$ and ${\it\phi}_{h}$ are assumed to be linear functions of ${\bf\xi}$ , i.e.

(6.4a,b ) $$\begin{eqnarray}{\it\phi}_{m}({\bf\xi})=1+{\it\beta}{\bf\xi}\quad \text{and}\quad {\it\phi}_{m}({\bf\xi})=\mathit{Pr}_{T}+{\it\beta}{\bf\xi},\end{eqnarray}$$

with $\mathit{Pr}_{T}$ being the turbulent Prandtl number and ${\it\beta}$ a dimensionless constant. Despite its simplicity, Monin–Obukhov self-similarity theory has been successfully compared with in situ field measurements (Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971; Kaimal et al. Reference Kaimal, Wyngaard, Haugen, Coté, Izumi, Caughey and Readings1976) and numerical simulations (García-Villalba & del Álamo Reference García-Villalba and del Álamo2011; Ansorge & Mellado Reference Ansorge and Mellado2014; Deusebio et al. Reference Deusebio, Brethouwer, Schlatter and Lindborg2014). In the following we develop a model based on it. By integrating the mean gradients given by (6.3) together with definitions (6.4), we arrive at an expression for the velocity and buoyancy jumps between the upper and lower walls,

(6.5) $$\begin{eqnarray}U=2\frac{u_{{\it\tau}}}{C_{{\it\kappa}}}\left(\log \mathit{Re}_{{\it\tau}}+{\it\beta}\frac{h}{L}+C_{{\it\kappa}}C_{1}\right)\end{eqnarray}$$

and

(6.6) $$\begin{eqnarray}T=2\frac{T_{{\it\tau}}}{C_{{\it\kappa}}}\left(\mathit{Pr}_{T}\log \mathit{Re}_{{\it\tau}}+{\it\beta}\frac{h}{L}+C_{{\it\kappa}}C_{2}\right),\end{eqnarray}$$

with $C_{1}$ and $C_{2}$ being constants. We note that Monin–Obukhov theory assumes a very similar shape for the velocity and temperature mean profiles, consistently with our simulations, as shown by the empirically observed scaling $\mathit{Ri}_{f}^{G}\sim \mathit{Ri}$ . Equations (6.5) and (6.6) can be rearranged to provide expressions relating $\mathit{Re}_{{\it\tau}}$ , $\mathit{Re}$ , $L$ and $\mathit{Ri}$ , i.e.

(6.7) $$\begin{eqnarray}\mathit{Re}=\frac{\mathit{Re}_{{\it\tau}}}{C_{{\it\kappa}}}\left(\log \mathit{Re}_{{\it\tau}}+{\it\beta}\frac{\mathit{Re}_{{\it\tau}}}{L^{+}}+C_{1}\right)\end{eqnarray}$$

and

(6.8) $$\begin{eqnarray}\mathit{Ri}=\frac{h}{L}\frac{\displaystyle \mathit{Pr}_{T}\log \mathit{Re}_{{\it\tau}}+{\it\beta}\frac{\mathit{Re}_{{\it\tau}}}{L^{+}}+C_{2}}{\left(\log \mathit{Re}_{{\it\tau}}+{\it\beta}{\displaystyle \frac{\mathit{Re}_{{\it\tau}}}{L^{+}}}+C_{1}\right)^{2}},\end{eqnarray}$$

from which the value of $L^{+}$ can be estimated. Here, we set ${\it\beta}=4.8$ , $\mathit{Pr}_{T}=0.7$ and $C_{1}=C_{2}=5.5$ in (6.7) and (6.8), similar to the values suggested by Wyngaard (Reference Wyngaard2010) and consistent with our simulations, indicating a $\mathit{Pr}_{T}$ slightly smaller than unity.

Figure 16(a) shows a comparison between the velocity profiles from the numerical simulations (already presented in figure 3) and the predictions provided by the Monin–Obukhov theory. Despite its simplicity, we observe a reasonable agreement, particularly in the core region of the flow. On the other hand, the agreement close to the walls is less good (especially for the unstratified simulation 6 for which $\mathit{Re}$ may well be too small to allow the development of a proper logarithmic layer), resulting in an error in the prediction of $\mathit{Re}_{{\it\tau}}$ and $L^{+}$ . In figure 16(b) we plot the percentage error difference,

(6.9) $$\begin{eqnarray}\frac{|L_{DNS}^{+}-L_{M}^{+}|}{L_{DNS}^{+}}\,\%,\end{eqnarray}$$

between the value obtained from the numerical simulations, $L_{DNS}^{+}$ , and the estimated value of the model using (6.7) and (6.8), $L_{M}^{+}$ . It is worth noting that the error in the prediction of $L^{+}$ tends to be larger when intermittency appears in the flow. In fully developed turbulent regimes the predicted value of $L^{+}$ is within a relative error of 20–40 %, which we believe is acceptable.

Figure 16.

(a) Comparison between numerical simulations (——) and the analytical model (- - - -) defined in (6.5) and (6.6). The three profiles correspond to $\langle \overline{u}\rangle$ velocity profiles for the simulations 6, 11 and 13, shown in figure 3. (b) Percentage errors in the estimation of $L^{+}$ using the Monin–Obukhov model, as defined in (6.9). The symbols represent all of the simulations given in table 1 and are coloured by the turbulent fraction ${\it\gamma}$ , ranging from blue (laminar) to red (fully turbulent).

6.2. Explicit algebraic models

The more complex model we consider here is an EAM for turbulent Reynolds stresses and heat fluxes for stably stratified flows. An EAM is generally developed from differential Reynolds-stress models (DRSMs) in which the evolution equations for the Reynolds stresses and heat fluxes take the form (Pope Reference Pope2000)

(6.10) $$\begin{eqnarray}\frac{\text{D}\overline{u_{i}u_{j}}}{\text{D}t}-\mathscr{D}_{ij}=\mathscr{P}_{ij}+{\it\Pi}_{ij}-{\it\varepsilon}_{ij}+\mathscr{G}_{ij}\end{eqnarray}$$

and

(6.11) $$\begin{eqnarray}\frac{\text{D}\overline{u_{i}{\it\theta}}}{\text{D}t}-\mathscr{D}_{{\it\theta}i}=\mathscr{P}_{{\it\theta}i}+{\it\Pi}_{{\it\theta}i}-{\it\varepsilon}_{{\it\theta}i}+\mathscr{G}_{{\it\theta}i},\end{eqnarray}$$

where the terms correspond to (from left to right) advection, diffusion, production, pressure redistributions, dissipation and buoyancy. Following Lazeroms et al. (Reference Lazeroms, Brethouwer, Wallin and Johansson2013), equations (6.10) and (6.11) can be recast as evolution equations for the dimensionless counterparts $\unicode[STIX]{x1D622}_{ij}=\overline{u_{i}u_{j}}/K-2/3{\it\delta}_{ij}$ and ${\it\xi}_{i}=\overline{u_{i}{\it\theta}}/\sqrt{K\,K_{{\it\theta}}}$ , where $K$ is the turbulent kinetic energy, $\overline{u_{i}u_{i}}/2$ , and $K_{{\it\theta}}$ is the scalar variance, $\overline{{\it\theta}^{2}}/2$ . In EAMs, the weak-equilibrium assumption (Rodi Reference Rodi1976) is generally made at this stage, neglecting advection and diffusion terms in the evolution equations for $\unicode[STIX]{x1D622}_{ij}$ and ${\it\xi}_{i}$ . This allows one to reduce these equations to an algebraic form in which time derivatives and higher-order spatial derivatives of Reynolds stresses disappear, i.e.

(6.12) $$\begin{eqnarray}\frac{\overline{u_{i}u_{j}}}{K}(\mathscr{P}-{\it\varepsilon}+\mathscr{G})=\mathscr{P}_{ij}+{\it\Pi}_{ij}-{\it\varepsilon}_{ij}+\mathscr{G}_{ij}\end{eqnarray}$$

and

(6.13) $$\begin{eqnarray}\frac{\overline{u_{i}{\it\theta}}}{2}\left(\frac{\mathscr{P}-{\it\varepsilon}+\mathscr{G}}{K}+\frac{\mathscr{P}_{{\it\theta}}-{\it\varepsilon}_{{\it\theta}}}{K_{{\it\theta}}}\right)=\mathscr{P}_{i{\it\theta}}+{\it\Pi}_{i{\it\theta}}-{\it\varepsilon}_{i{\it\theta}}+\mathscr{G}_{ij}.\end{eqnarray}$$

In (6.12) and (6.13), the terms ${\it\Pi}_{ij}$ and ${\it\Pi}_{i{\it\theta}}-{\it\varepsilon}_{i{\it\theta}}$ can be modelled following Launder (Reference Launder1975), with expressions involving $K$ , $K_{{\it\theta}}$ and the other variables, i.e. $U_{i}$ , $T$ , $\unicode[STIX]{x1D622}_{ij}$ and ${\it\xi}_{i}$ . The dissipation term ${\it\varepsilon}_{ij}$ is assumed to be isotropic, i.e. ${\it\varepsilon}_{ij}=(2/3){\it\varepsilon}{\it\delta}_{ij}$ , where ${\it\varepsilon}$ is the turbulent kinetic energy dissipation. Equations (6.12) and (6.13) are implicit in $\unicode[STIX]{x1D622}_{ij}$ and ${\bf\xi}$ and highly nonlinear. However, explicit algebraic equations relating $\unicode[STIX]{x1D622}_{ij}$ and ${\bf\xi}$ to the other mean quantities can be derived by expanding $\unicode[STIX]{x1D622}_{ij}$ and ${\bf\xi}$ in tensor groups (Pope Reference Pope1975) and using an approximate equation for the kinetic energy production-to-dissipation ratio (Lazeroms et al. Reference Lazeroms, Brethouwer, Wallin and Johansson2015). These relations together with suitable expressions for $K$ , $K_{{\it\theta}}$ and ${\it\varepsilon}$ such as the $K{-}{\it\omega}$ model (Menter Reference Menter1992) form a closed set of equations which can be numerically solved in time. For further details on the derivation of the model we refer the interested reader to Lazeroms et al. (Reference Lazeroms, Brethouwer, Wallin and Johansson2013, Reference Lazeroms, Brethouwer, Wallin and Johansson2015).

Explicit algebraic models have been successfully applied to stably stratified shear parallel flows and stably stratified channel flows, and provide fair agreement with simulations and significant improvements with respect to standard eddy-viscosity and eddy-diffusivity models (Lazeroms et al. Reference Lazeroms, Brethouwer, Wallin and Johansson2013). In the following, we extend this comparison to PCF. It is worth emphasising here that all of the modelling constants are the same as in Lazeroms et al. (Reference Lazeroms, Brethouwer, Wallin and Johansson2013). The velocity (figure 17 a) and density profiles (not shown) agree reasonably well with the numerical simulations over the entire channel depth. Figure 17(b) shows the percentage error, as defined in (6.9) using EAMs for estimating the model value $L_{M}^{+}$ . The agreement is somewhat poor at low $\mathit{Re}$ but improves for increasing $\mathit{Re}$ , providing good predictions for $\mathit{Re}>5\times 10^{3}$ . That models and simulations do not agree at low $\mathit{Re}<5\times 10^{3}$ is not surprising since at such low $\mathit{Re}$ diffusive effects may well be important and the weak-equilibrium assumption may be inappropriate. Nevertheless, as long as $\mathit{Re}$ is large enough, the relative errors for estimating $\mathit{Re}_{{\it\tau}}$ and $L^{+}$ are approximately $10\,\%$ and $20\,\%$ respectively, slightly better than for the Monin–Obukhov self-similarity theory.

Figure 17.

(a) Comparison between numerical simulations (——) and the EAM (- - - -) from Lazeroms et al. (Reference Lazeroms, Brethouwer, Wallin and Johansson2013). The three profiles correspond to $\langle \overline{u}\rangle$ for the simulations 6, 11 and 13, shown in figure 3. (b) Percentage errors in the estimation of $L^{+}$ using the analytical model. The symbols represent the simulations given in table 1 and are coloured by the turbulent fraction ${\it\gamma}$ , ranging from blue (laminar) to red (fully turbulent).

6.3. The intermittency boundary

We now use the models outlined in §§ 6.1 and 6.2 to predict the intermittency boundary in the $\mathit{Re}{-}\mathit{Ri}$ plane between turbulent and intermittent dynamics. According to figure 11, we predict that intermittency onsets at $L^{+}\approx 200$ . Using the two models, we have estimated $L^{+}$ as a function of $\mathit{Re}$ and $\mathit{Ri}$ , with $\mathit{Re}$ between 10³ and 10⁶ and $\mathit{Ri}$ between 0 and 0.2. A summary of the simulations together with the intermittency boundaries predicted by the two models is shown in figure 18. The colours represent the turbulent fraction ${\it\gamma}$ . The two models show an overall similar behaviour, and they are indistinguishable for most of the $\mathit{Re}$ range. Remarkably, the model based on the Monin–Obukhov self-similarity theory is able to provide good agreement with the numerical simulations and with the more complex EAM. It is nevertheless worth noting that PCF is a simple and idealised geometry in which, in statistically steady states, momentum and heat fluxes are constant across the channel. Analytical models may provide somewhat poorer performance in more complex flow geometries, where more refined models, such as EAM, have already been shown to provide better results. As shown in figure 18, at least at smaller $\mathit{Re}$ and $\mathit{Ri}$ , the intermittency and relaminarisation boundaries are quite close. It is still an open question whether this observation carries over to larger values of $\mathit{Re}$ and $\mathit{Ri}$ , particularly in sufficiently large computational domains.

Figure 18.

Summary of the direct numerical simulations (circles) and predictions of the intermittency boundary given by the Monin–Obukhov theory (solid black line) and the EAM (dashed grey line). The symbols are coloured by the turbulent fraction ${\it\gamma}$ , ranging from blue (laminar) to red (fully turbulent). Triangles representing simulations that fully relaminarised are also shown for reference.