4.1. Development of stably stratified grid turbulence

Figure 2 visualizes the development of grid turbulence in R10F1. The left panels show the colour contour of $u$ on the surfaces of the computational domain while the right panels visualize $u$ on a small part of an $x$–$y$ plane. At $t/t_r=5$ in ($a$), the instantaneous velocity profiles are strongly influenced by the wakes of the vertical plates. The wakes develop with time and interact with each other, resulting in the formation of the fully developed grid turbulence. Figure 3 visualizes $u$ on a vertical plane at $t/t_r=15$ in R10F01 and R10F6. The strongly stable stratification at $Fr_M=0.1$ causes highly anisotropic velocity fluctuations. However, more isotropic distribution of $u$ is found for $Fr_M=6$. For both $Fr_M$, the imprints of the wakes are hardly seen in the velocity profiles, and the grid turbulence has developed by this time.

Figure 2.

Development of temporally evolving grid turbulence in R10F1 at (a) $t/t_r=5$, (b) $t/t_r=10$ and (c) $t/t_r=15$. The streamwise velocity $u$ is visualized on the boundaries of the computational domain on the left panels while two-dimensional contours of $u$ on an $x$–$y$ plane are shown on the right panels.

Figure 3.

Streamwise velocity profiles on a vertical plane in (a) R10F01 and (b) R10F6 at $t/t_r=15$.

Figure 4 compares the temporal evolution of volume-averaged horizontal kinetic energy $\langle u^2+v^2\rangle _V/2$ between the DNS and experiments of towed-grid experiments by Praud et al. (Reference Praud, Fincham and Sommeria2005), where the average is taken in the computational domain of DNS or the measurement volume of scanning particle image velocimetry. The experiments were conducted in a water tank, where the density stratification was generated with saltwater. Therefore, the Schmidt number was much larger than 1 in the experiments while all DNS consider $Pr=1$, and this difference may affect the comparison between the DNS and experiments. Figure 4 shows that the horizontal kinetic energy decreases with time. The present DNS with small $Fr_M$ agrees well with the experiments, which were also conducted with $Fr_M\approx 0.1$. The horizontal kinetic energy contains the contribution from both mean and fluctuating velocity components. When the grid turbulence is generated, the kinetic energy in the mean velocity is converted to that of the velocity fluctuations by the production term of turbulent kinetic energy, which is expressed as the product of the Reynolds stress and the mean velocity gradient. The decay of horizontal velocity fluctuations is caused by the viscous dissipation once grid turbulence has fully developed. Although the mean velocity profile given as the initial condition of DNS is only a simple approximation of the flow behind a grid, the initial decay of $\langle u^2+v^2\rangle _V/2$ is in good agreement between the DNS and the experiments. In the DNS the horizontal kinetic energy tends to be small for large $Fr_M$, while the $Re_M$ dependence is not significant for a fixed value of $Fr_M$. The energy production by the mean velocity gradient contributes to the growth of the streamwise velocity component of turbulent kinetic energy, which is converted to the other components by the pressure-strain correlation (Pope Reference Pope2000). This conversion to the vertical component is suppressed in a stably stratified fluid, and a smaller amount of horizontal kinetic energy is converted to the vertical component at lower $Fr_M$ during the production phase of grid turbulence. This stratification effect explains larger $\langle u^2+v^2\rangle _V/2$ for smaller $Fr_M$.

Figure 4.

Temporal evolutions of volume-averaged horizontal kinetic energy $\langle u^2+v^2\rangle _V/2$ compared with towed-grid experiments by Praud et al. (Reference Praud, Fincham and Sommeria2005).

Figure 5 shows the lateral profiles of mean velocity $\langle u\rangle$ and r.m.s. velocity fluctuations $u_{rms}=\langle u'^2\rangle ^{1/2}$. The turbulent wake develops for $0\leq y/M \leq 0.1$ and $0.9\leq y/M \leq 1$ and the mean velocity difference between the wake and the open region of the grid becomes small with time. The mean velocity profile is homogeneous at $t/t_r=25$. The velocity fluctuations in the $x$ direction are generated by the mean velocity gradient, which is large at $y/M=0.1$ and 0.9 at $t/t_r=2$. Therefore, $u_{rms}$ also reaches peaks at these locations. The wakes spread in the $y$ direction and merge with the adjacent wakes. The profile of $u_{rms}$ also becomes homogeneous at $t/t_r=25$, and $u_{rms}$ decays with time.

Figure 5.

Lateral profiles of (a) mean velocity $\langle u\rangle$ and (b) r.m.s. velocity fluctuations $u_{rms}$ in R10F1. The results are taken at five instances of $t/t_r=2$, 6, 12, 25 and 50.

Figure 6 compares the temporal evolutions of the r.m.s. fluctuations of three velocity components ($u_{rms}$, $v_{rms}$ and $w_{rms}$) at $y/M=0$ and 0.5. Here, the results are also compared among different $Fr_M$ for $Re_M=10\,000$. As the turbulence is generated at an early time, the r.m.s. velocity fluctuations increase with time. This increase is more rapid at $y/M=0$ than at $y/M=0.5$ because the turbulent wake of each plate develops at $y/M=0$. Then, they begin to decrease as the turbulence decays. At a late time, the r.m.s. velocity fluctuations hardly differ at $y/M=0$ and $0.5$ and the flow reaches a statistically homogeneous state. Furthermore, the horizontal velocity fluctuations become isotropic ($u_{rms}\approx v_{rms}$). The time required for stably stratified grid turbulence to be homogeneous and isotropic in the horizontal directions is different depending on $Fr_M$. This time is about $10t_r$, $20t_r$, $25t_r$ and $10t_r$ for $Fr_M=0.1$, 0.4, 1 and 6, respectively. For $Fr=0.1$ and 6, the r.m.s. velocity fluctuations monotonically decrease after they reach the peaks. However, $u_{rms}$ and $v_{rms}$ for $Fr_M=0.4$ and 1 hardly decrease with time for $7\lesssim t/t_r \lesssim 12$ after the rapid decay from the peak at $t/t_r\approx 4$. This tendency is clearer for $y/M=0$ in the wakes of vertical plates. Therefore, the slow decay for $7\lesssim t/t_r \lesssim 12$ might be related to the decay characteristics of the wake.

Figure 6.

Temporal evolutions of r.m.s. velocity fluctuations at $y/M=0$ and 0.5 in the cases of (a) R10F01, (b) R10F04, (c) R10F1 and (d) R10F6.

Figure 7(a) presents $u_{rms}/U_0$ as a function of time normalized by the buoyancy period $N^{-1}$. A comparison between figures 6 and 7(a) suggests that the initial growth of the r.m.s. velocity fluctuations is better characterized by $t_r$. The time interval corresponding to the slow decay of $u_{rms}$ is found for $14\lesssim tN \lesssim 30$ for $Fr=0.4$ and $7\lesssim tN \lesssim 20$ for $Fr=1$. Experiments and numerical simulations of a turbulent wake of a sphere in a stably stratified fluid have reported that the wake development can be divided into three regimes: a three-dimensional regime for $0\leq tN \leq 2$, a non-equilibrium (NEQ) regime for $2\leq tN \leq 50$ and a quasi-two-dimensional regime for $50\leq tN$ (Spedding Reference Spedding1997; Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011). The decay rates of the mean velocity deficit and r.m.s. velocity fluctuations are different among these regimes. Zhou & Diamessis (Reference Zhou and Diamessis2019) showed that the decay of the horizontal r.m.s. velocity fluctuations becomes slow for $10\lesssim tN\lesssim 50$, which is within the NEQ regime. In the present DNS, the range of $10\leq tN\leq 50$ corresponds to $4\leq t/t_r \leq 20$ for $Fr_M=0.4$ and $10\leq t/t_r \leq 50$ for $Fr_M=1$. The slow decays of $u_{rms}$ and $v_{rms}$ are observed within these intervals. In contrast, the NEQ regime corresponds to $1\leq t/t_r \leq 5$ for $Fr_M=0.1$ and $60\leq t/t_r \leq 300$ for $Fr_M=6$. When the turbulent wake is being generated, $u_{rms}$ at $y=0$ increases with time. Once the turbulent wake has developed, $u_{rms}$ begins to decay and the interaction of the turbulent wakes generates homogeneous turbulence, where $u_{rms}$ does not depend on $y$. In figure 7(a), $u_{rms}$ in R10F01 increases with time for the time corresponding to the possible NEQ regime ($10\lesssim tN\lesssim 50$). Therefore, the NEQ regime is over before the turbulent wake has developed. On the other hand, $u_{rms}$ in R10F6 hardly differs for $y/M=0$ and 0.5 in the NEQ regime and decreases with time, suggesting that homogeneous turbulence has developed for $Fr_M=6$ before the wake reaches the NEQ regime. Thus, the slow decay during the NEQ regime of stratified wakes is not observed for $Fr_M=0.1$ and 6. It should be noted that the discussions are based on the studies on the stratified wake of a sphere, which is different from the wake of a vertical plate considered in this study. The slow decay in the NEQ regime is only one of the possible explanations for the behaviour $u_{rms}$ and $v_{rms}$.

Figure 7.

Temporal evolutions of (a) $u_{rms}/U_0$ at $y/M=0$ and 0.5, and (b) the mean velocity deficit of the wake, $\Delta u$, which are plotted as functions of $tN$.

The wake development of each plate can also be evaluated with the mean velocity deficit $\Delta u$, which is defined as the difference in $\langle u\rangle$ between $y/M=0$ and 1. Figure 7(b) shows $\Delta u$ as a function of $tN$. In a freely developing wake of an object, such as a sphere and a cylinder, $\Delta u$ becomes small with time although $\Delta u$ does not reach zero. The wakes are generated by many plates in the case of grid turbulence. Therefore, $\Delta u$ in grid turbulence rapidly decreases to zero with time. In a stably stratified wake of a sphere, the decay of $\Delta u$ becomes slower once the flow enters the NEQ regime (Spedding Reference Spedding1997). The decay of $\Delta u$ for R10F04 and R10F1 becomes slightly slower at $tN\approx 10$ and 5, respectively, than at an earlier time. However, the change of the decay rate is less clear than that reported for stably stratified wakes. This might be partially due to the interaction with adjacent wakes, which results in $\Delta u\rightarrow 0$ while $\Delta u$ in a conventional wake is greater than 0 even at a late time.

The budget of the turbulent kinetic energy $k_T(\kern0.7pt y,t)= ( \langle u'^2\rangle +\langle v'^2\rangle +\langle w'^2\rangle )/2$ is examined with the transport equation of $k_T$. Once the grid turbulence has developed, the production rate of $k_T$ becomes much smaller than the dissipation rate. Here, an average in the $y$ direction is defined for a statistical quantity $F(\kern0.7pt y,t)$ as $\bar {F}(t)=(1/L_y)\int F(\kern0.7pt y,t)\,{{\rm d}y}$. The spatially averaged turbulent kinetic energy $\overline {k_T}$ evolves following the equation given by

(4.1)\begin{equation} \frac{{\rm d} \overline{k_T}}{{\rm d} t} =\bar{P}-\bar{\varepsilon} + \bar{B}, \end{equation}

where $P$, $\varepsilon$ and $B$ are the production, dissipation and buoyancy flux in the transport equation of $k_T(\kern0.7pt y,t)$,

(4.2a–c)\begin{equation} P={-} \langle u'_{i}u'_{j} \rangle \frac{\partial \langle u_i\rangle}{\partial x_j},\quad \varepsilon= \nu\left\langle\frac{\partial u'_{i}}{\partial x_j} \frac{\partial u'_{i}}{\partial x_j}\right\rangle,\quad B={-} \frac{g}{\rho_0}\langle w'\rho'\rangle. \end{equation}

Here, $\varepsilon$ and $B$ are identical to the first and second terms of (2.6), respectively, when the flow is statistically homogeneous. Figure 8(a) shows the temporal evolutions of $\bar {P}$, $\bar {\varepsilon }$ and $\bar {B}$ for R10F1. Initially, $\bar {P}$ is larger than $\bar {\varepsilon }$ until $t/t_r\approx 3$ because the turbulent kinetic energy is initially small and increases due to the production term. Until grid turbulence has developed by the wake interaction, $\bar {P}$ and $\bar {\varepsilon }$ are comparable for $3\lesssim t/t_r\lesssim 10$. This relationship of $\bar {P}\sim \bar {\varepsilon }$ is widely observed in turbulent free shear flows, where the ratio between the production and dissipation terms of turbulent kinetic energy is close to 1 (Panchapakesan & Lumley Reference Panchapakesan and Lumley1993; Taveira & da Silva Reference Taveira and da Silva2013; Watanabe et al. Reference Watanabe, Sakai, Nagata and Ito2016). Then, $\bar {P}$ rapidly decreases because the mean velocity gradient ${\partial \langle u\rangle }/{\partial y}$ and the Reynolds stress $\langle u'v'\rangle$ become small as grid turbulence develops into a statistically homogeneous state. Because the buoyancy flux $\bar {B}$ is always negative, $-\bar {B}$ is shown in the figure. Positive values of $-\bar {B}$ suggest that the turbulent kinetic energy is converted to the potential energy. However, $-\bar {B}$ is smaller than $\bar {\varepsilon }$ and the decay of the turbulent kinetic energy is dominated by $\bar {\varepsilon }$. In all DNS, $\bar {\varepsilon }$ is at least 10 times larger than $-\bar {B}$ in the decay period of grid turbulence. Figure 8(b) shows the temporal evolutions of $\bar {P}/\bar {\varepsilon }$. The ratio with $\bar {P}/\bar {\varepsilon }\approx 10^{0}$ is observed for $3\lesssim t/t_r\lesssim 10$ except for R10F01, for which $\bar {P}/\bar {\varepsilon }$ begins to decay at an earlier time than for other cases. In all DNS, $\bar {P}/\bar {\varepsilon }$ rapidly decreases because of the development of grid turbulence. However, a comparison between different $Fr_M$ suggests that the rapid decrease of $\bar {P}/\bar {\varepsilon }$ occurs at an earlier time for lower $Fr_M$. Table 1 includes time $t_H$ at which $\bar {P}/\bar {\varepsilon }$ becomes smaller than 0.01. Beyond this time, turbulence is nearly homogeneous and the production of turbulent kinetic energy is almost negligible. The Reynolds number hardly affects $t_H$, which increases with $Fr_M$. Therefore, the grid turbulence reaches a homogeneous state faster under stronger stable stratification. Figures 9(a,b) present the temporal evolutions of $Fr_H=u_{rms}/L_uN$ and $Re_H=u_{rms}L_u/\nu$ at $y/M=0$, where $L_u$ is the longitudinal integral scale defined with $u$. After $Fr_H$ increases by the development of grid turbulence, it monotonically decays. Comparison among all DNS cases suggests that $Fr_H$ depends on $Fr_M$ and does not exhibit a strong dependence on $Re_M$. Similarly, $Re_H$ reaches a peak when turbulence is generated. However, the decay of $Re_H$ is much slower than that of $Fr_H$ because $L_u$ increases with time. Figure 9(c) shows the buoyancy Reynolds number $Re_b$ at $y/M=0$. Except for $Fr_M=0.1$, $Re_b$ exceeds $10$ during the turbulence-production phase. Then, $Re_b$ rapidly decreases to $Re_b<10^0$ for $Fr_M=0.4$ and 1. However, the grid turbulence with $Fr_M=6$ has $Re_b>10^0$ until $t/t_r\approx 70$, and small-scale turbulent motions are active in the early decay period. For $Fr_M=0.1$, $Re_b$ does not exceed $10^0$ even at an early time, and three-dimensional small-scale turbulence is inhibited even when the turbulence is generated.

Figure 8.

Temporal evolutions of (a) the spatially averaged budget of turbulent kinetic energy for R10F1 and (b) the ratio between the production and dissipation terms $\bar {P}/\bar {\varepsilon }$.

Figure 9.

Temporal evolutions of (a) horizontal Froude number $Fr_H$, (b) Reynolds number based on the horizontal integral scale $Re_H$ and (c) buoyancy Reynolds number $Re_b$ at $y/M=0$. (d) The relation between $Re_H$ and $1/Fr_H$ during the decay of grid turbulence at $y/M=0$.

It is useful to diagnose stably stratified turbulence in the parameter space of $(Re_H, Fr_H^{-1})$ (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). Figure 9(d) shows the trajectory of $(Re_H, Fr_H^{-1})$ during the decay at $y/M=0$ after the peak of $Re_b$. The broken lines distinguish four regimes of stratified turbulence proposed in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007): viscosity-affected stratified flow, strongly stratified turbulence, weakly stratified turbulence and classical Kolmogorov turbulence. The two horizontal broken lines represent $Fr_H=1$ and 0.02 while the diagonal broken line is $Re_{H}Fr_{H}^2=1$. As discussed above, $Re_{H}Fr_{H}^2$ is close to $Re_b$ when $\alpha =\varepsilon /(U_{H}^3/L_H)\sim 10^{0}$, and the plots of $(Re_H, Fr_H^{-1})$ complement the temporal evolution of $Re_b$. For $(Re_M, Fr_M)=(10\,000,6)$, the decay of turbulence begins in the classical Kolmogorov turbulence regime, enters the weakly stratified turbulence regime and then reaches the viscosity-affected stratified flow regime. For the other cases, the turbulence is generated in the weakly stratified turbulence regime and the decay proceeds in the viscosity-affected stratified flow regime. The decay of turbulence is examined as a function of time in § 4.2. The time at which the flow enters the viscosity-affected stratified flow regime $(t_V)$ is given in table 1. Time $t_V$ is normalized by the reference time scale of grid turbulence $t_r$ or the buoyancy period $N^{-1}$. For $Re_M=10\,000$, $t_VN$ is 80 for both $Fr_M=0.1$ and 0.4. However, $t_V$ is not determined solely by $t_r$ or $N^{-1}$.

Figure 10 shows the temporal evolutions of the vertical Froude number $Fr_V=U_H/L_VN$. The theories of strongly stratified turbulence with high $Re_b$ assume that the stable stratification is sufficiently strong, and require large $Re_H$ and small $Fr_H$ in addition to $Re_b\sim Re_H Fr_H^2\gg 1$. The decay laws for high $Re_b$ are obtained with (2.12), which assumes that $Fr_V$ is independent of time. de Bruyn Kops & Riley (Reference de Bruyn Kops and Riley2019) has observed that $Fr_V$ hardly varies with time in decaying turbulence with stable stratification at high $Re_H$ and low $Fr_H$. In figure 10, $Fr_V$ reaches a maximum value at $t/t_r\approx 10^1$ and decays with time. Therefore, constant $Fr_V$ is not observed in the DNS even though $Re_b$ exceeds $10^2$ for large $Fr_M$ cases. This is explained by the weak stratification in these cases, and high $Re_b$ is not a sufficient condition to apply the theories of strongly stratified turbulence.

Figure 10.

Temporal evolutions of vertical Froude number $Fr_V$ at $y/M=0$.

4.2. Decay properties of temporally evolving, stably stratified grid turbulence

The decay of temporally evolving grid turbulence is compared with the theories for Saffman and Batchelor turbulence. The velocity and length scales considered in § 2 are defined as

(4.3a–c)\begin{equation} U_{H} = (\langle u'^2\rangle + \langle v'^2\rangle)/2,\quad L_{H} = (L_{uH} + L_{vH})/2,\quad L_{V} = (L_{uV} + L_{vV})/2, \end{equation}

where the horizontal and vertical integral scales of each velocity component are defined as

(4.4a,b)$$\begin{gather} L_{uH}=\int \frac{\langle u'(x,y,z)u'(x+r_x,y,z)\rangle}{\langle {u'}^2\rangle} {\rm d} r_x,\quad L_{vH} =\int \frac{\langle v'(x,y,z)v'(x,y+r_y,z)\rangle}{\langle {v'}^2\rangle} {\rm d} r_y, \end{gather}$$

(4.5a,b)$$\begin{gather}L_{uV}=\int \frac{\langle u'(x,y,z)u'(x,y,z+r_z)\rangle}{\langle {u'}^2\rangle} {\rm d} r_z,\quad L_{vV} =\int \frac{\langle v'(x,y,z)v'(x,y,z+r_z)\rangle}{\langle {v'}^2\rangle} {\rm d} r_z. \end{gather}$$

These statistics depend on $y$ and $t$. Hereafter, the results are presented only for $y/M=0$ because the statistics are homogeneous during the decay for which the DNS results are compared with the theories.

Time independence of $U_{H}^2 L_H^4 L_V$ and $U_{H}^2 L_H^2 L_V$ considered in § 2 is derived by the assumption of ‘partial self-similarity’, which means the self-similarity for scales of the order of integral scales (Davidson Reference Davidson2010). If the partial self-similarity is applied, the autocorrelation functions in (4.4a,b) and (4.5a,b) are functions of the separation distance normalized by the integral scale, e.g. $r_x/L_{uH}$. The assumption of the partial self-similarity is examined with the autocorrelation functions of $u$, which are denoted by $f_{uH}(r_x,y,t)=\langle u'(x,y,z,t)u'(x+r_x,y,z,t)\rangle /\langle u'^{2}\rangle$ and $f_{uV}(r_z,y,t)=\langle u'(x,y,z,t)u'(x,y,z+r_z,t)\rangle /\langle u'^{2}\rangle$. Figure 11 plots $f_{uH}$ and $f_{uV}$ at $y=0$ as functions of $r_x/L_{uH}$ and $r_z/L_{uV}$, respectively, in R10F01. For $t/t_r=10$ and 20, the profiles of $f_{uH}$ vary with time. Similarly, $f_{uV}$ is slightly different between $t/t_r=10$ and later times. However, the curves collapse well at later times, and $f_{uH}$ and $f_{uV}$ are functions of $r_x/L_{uH}$ and $r_z/L_{uV}$, respectively, which do not depend on time. The results suggest that the partial self-similarity for large scales is valid except for a very early time.

Figure 11.

Temporal evolutions of the autocorrelation functions of streamwise velocity defined with separation distances in (a) the streamwise direction, $r_x$, and (b) the vertical direction, $r_z$. The results are taken at $y=0$ for R10F01. The distance is normalized by the integral scale, $L_{uH}$ or $L_{uV}$.

Saffman turbulence and Batchelor turbulence have a three-dimensional energy spectrum with $E(k)\sim k^2$ and $k^4$ at large scales, respectively. Axisymmetric turbulence, such as homogeneous turbulence subject to stable stratification or rotation, also has these spectral shapes for Saffman turbulence and Batchelor turbulence. The spectrum in temporally evolving grid turbulence is evaluated with a spherical-shell average following Valente, da Silva & Pinho (Reference Valente, da Silva and Pinho2016). The Fourier transform is applied to the instantaneous velocity $u_{i}(x,y,z,t)$ in the $x$, $y$ and $z$ directions to obtain the Fourier coefficient $\hat {u}_{i}(k_x,k_y,k_z,t)$. The energy spectrum is computed as

(4.6)\begin{equation} E(k,t)=4{\rm \pi} k^2\langle \tfrac{1}{2} \hat{u}_{i}\hat{u}_{i}^{*}\rangle_{k},\end{equation}

where $\hat {u}_{i}^{*}$ is the complex conjugate of $\hat {u}_{i}$, $k$ is the wavenumber and the spherical-shell average $\langle\, f\rangle _{k}$ of $f(k_x, k_y, k_z)$ is defined as

(4.7)\begin{equation} \left\langle\, f \right\rangle_{k}(k)=\frac{1}{N_k} \sum_{|k'-k|<\frac{\Delta k}{2}} f(k_x, k_y, k_z),\quad \text{with } k'=\sqrt{k_x^2+ k_y^2+k_z^2}. \end{equation}

Here, the summation is taken for $N_k$ Fourier coefficients whose wavenumber $k'$ satisfies $|k'-k|<\Delta k/2$. The shell thickness is set to be $\Delta k = 0.296M^{-1}$, which is about 1.4 times larger than the wavenumber of the horizontal domain size $2{\rm \pi} /L_x=0.209M^{-1}$.

Figure 12 shows $E(k)$ of the initial velocity field, where the mean velocity profile is given by (3.5). The wavenumber corresponding to the mesh size $M$ is $k_M=2{\rm \pi} /M$. The inset shows $E(k)$ at a low-wavenumber range. Here, $E(k)$ has distinct peaks at wavenumbers related to $M$, namely $k=k_M$, $2k_M$, $3k_M$ and $4k_M$. The initial velocity field is laminar and does not have the energy spectrum of turbulence, such as $E(k)\sim k^2$ and $k^4$, which can be formed only after the grid turbulence is generated.

Figure 12.

A three-dimensional energy spectrum $E(k)$ at $t=0$.

The analysis of the turbulent kinetic energy budget suggests that homogeneous turbulence has developed at $t_H/t_r=12$–33, as summarized in table 1. The energy spectrum after $t/t_r=30$ is examined for the lowest and highest $Fr_M$ cases. Figure 13 shows the temporal evolutions of $E(k)$ for R10F01 and R10F6, which are compared with $E(k)\sim k^2$ and $k^4$. The oscillation of $E(k)$ is due to statistical convergence because $E(k)$ is evaluated from a single snapshot and small $\Delta k$ is used to examine $E(k)$ at large scales. For both cases, $E(k)$ for $kM\lesssim 3$ tends to follow $E(k)\sim k^2$ rather than $E(k)\sim k^4$. As turbulence decays with time, $E(k)$ at high wavenumbers becomes small. However, $E(k)$ at low wavenumbers hardly decay with time, suggesting that a coefficient $A$ in $E(k)\approx Ak^2$ is almost independent of time. The time independence of $A$ is related to the linear momentum conservation, which leads to the invariance of the Saffman integral (Davidson Reference Davidson2013). Therefore, temporally evolving grid turbulence in a stably stratified fluid evolves following the theory of axisymmetric Saffman turbulence. These results suggest that the mean velocity profile (3.5) results in the formation of Saffman turbulence with $E(k)\sim k^2$ in a stably stratified fluid.

Figure 13.

Temporal evolutions of the three-dimensional energy spectrum $E(k)$ for (a) R10F01 and (b) R10F6. The thin black lines represent $E(k)\sim k^2$ and $k^4$.

The invariants of Saffman and Batchelor turbulence are examined in figure 14, which presents the temporal evolutions of $U_{H}^{2}L_{H}^{2}L_V$ and $U_{H}^{2}L_{H}^{4}L_V$ at $y/M=0$. Once the turbulence is generated, $U_{H}^{2}L_{H}^{2}L_V$ hardly varies while $U_{H}^{2}L_{H}^{4}L_V$ increases with time. Experiments of non-stratified grid turbulence also observed that the invariant of Saffman turbulence, $U^{2}L^{3}$, hardly varies with the distance from the grid while the invariant of the Batchelor turbulence, $U^{2}L^{5}$, monotonically increases as the turbulence decays (Krogstad & Davidson Reference Krogstad and Davidson2010; Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014). These tendencies agree with the results for $U_{H}^{2}L_{H}^{2}L_V$ and $U_{H}^{2}L_{H}^{4}L_V$ defined for axisymmetric Saffman and Batchelor turbulence. Most wind-tunnel experiments of grid turbulence used conventional square grids, which consist of vertical and horizontal bars while this study considers the rake of vertical plates (Fincham et al. Reference Fincham, Maxworthy and Spedding1996; Praud et al. Reference Praud, Fincham and Sommeria2005). The time independence of $U_{H}^{2}L_{H}^{2}L_V$ indicates that grid turbulence generated by the vertical grid behaves as Saffman turbulence similarly to square-grid turbulence. Comparison of different DNS indicates that $U_{H}^{2}L_{H}^{2}L_V$ becomes independent of time sooner for smaller $Fr_M$. The turbulent kinetic energy budget in figure 8 also suggests a faster transition to homogeneous turbulence for smaller $Fr_M$. These results imply that $U_{H}^{2}L_{H}^{2}L_V=\mathrm {Const}.$ is achieved after the grid turbulence begins to freely decay without the production of turbulent kinetic energy.

Figure 14.

Temporal evolutions of (a) $U_{H}^{2}L_{H}^{2}L_V$ and (b) $U_{H}^{2}L_{H}^{4}L_V$ at $y/M=0$.

Figure 15(a) presents the temporal variation of the non-dimensional dissipation rate $\alpha =\varepsilon /(U_H^3/L_H)$, which is assumed to be a constant of order 1 in the derivation of the decay laws in § 2. After $t/t_r\gtrsim 30$, $\alpha$ hardly varies with time except for R10F6. Although $\alpha$ for R10F6 decreases after it attains the peak, $\alpha \approx 0.85$ tends to be constant for $t/t_r\gtrsim 100$. Although $\alpha$ hardly varies with time in each condition, the values of $\alpha$ depend on $Re_M$ and $Fr_M$. Therefore, $\alpha$ is not a universal constant in stratified grid turbulence at low $Re_b$. The values of $\alpha$ for $Fr_M=0.4$ and 1 slightly depend on $Re_M$, and $\alpha$ is slightly larger for $Re_M=5000$. This tendency agrees with DNS of decaying, stably stratified, homogeneous turbulence initialized with a prescribed energy spectrum, where the stationary values of $\alpha$ are larger for lower-Reynolds-number cases (Maffioli & Davidson Reference Maffioli and Davidson2016). As explained in § 2, $\varepsilon =\alpha (U_H^3/L_H)$ for low $Re_b$ is derived with $\varepsilon \sim \nu U_H^2/L_V^2$ and $L_V \sim (\nu L_H/U_H)^{1/2}$. As there are no reasons for prefactors $A$ and $B$ in $\varepsilon =A\nu U_H^2/L_V^2$ and $L_V =B(\nu L_H/U_H)^{1/2}$ to be universal constants which do not depend on flows, $\alpha =A/B^2$ should not be considered as a universal constant. Small temporal variations of $\alpha$ imply that these scalings, $\varepsilon \sim \nu U_H^2/L_V^2$ and $L_V \sim (\nu L_H/U_H)^{1/2}$, are valid during the decay. The decay laws in § 2 are derived based on these scalings, which are satisfied by the power laws in (2.19a–d) and (2.20a–d). The results presented below will also show that the temporal evolutions of $U_H$, $L_H$, $L_V$ and $\varepsilon$ agree with the decay laws of low-$Re_b$ Saffman turbulence except for R10F6, for which $\alpha$ becomes time independent only at a late time. Figure 15(b) shows the relationship between $\alpha$ and $Re_b$. Although $Re_b$ decreases with time, $\alpha$ weakly depends on time. As the flow is not in a strongly stratified turbulence regime, the scaling $\varepsilon \sim U_H^3/L_H$ requires the low-$Re_b$ scaling for $L_V$ and $\varepsilon$. Therefore, $\alpha$ tends to be time independent only when $Re_b$ is small. However, $\alpha$ in R10F01 varies with time even for $0.1\lesssim Re_b \lesssim 0.7$. In this range of $Re_b$, $\alpha$ for other cases varies slowly with $Re_b$. Filled circles in the figure mark $t=t_H$, before which the flow is statistically inhomogeneous. The time independence of $\alpha$ for R10F01 is achieved only after $t=t_H$. Therefore, $\alpha \approx$ Const. also requires the development of grid turbulence, and $\varepsilon \sim U_H^3/L_H$ seems invalid in the wake of each vertical plate.

Figure 15.

(a) Temporal evolutions of non-dimensional energy dissipation rate $\alpha =\varepsilon / (U_H^3/L_H)$ at $y/M=0$; (b) $\alpha$ plotted as a function of $Re_b$ at $y/M=0$. The time after $Re_b$ reaches a peak is shown in (b).

The dissipation scaling in non-stratified turbulence can be written as $\varepsilon =C_\varepsilon (U^3/L)$ with the r.m.s. velocity fluctuations $U$ and integral length scale $L$. The DNS of statistically steady isotropic turbulence has suggested that $C_\varepsilon =\varepsilon /(U^3/L)$ is a universal constant at a high Reynolds number (Ishihara, Gotoh & Kaneda Reference Ishihara, Gotoh and Kaneda2009). However, recent studies have confirmed that $C_\varepsilon$ differs depending on flows even at high $Re$ (Vassilicos Reference Vassilicos2015). Direct numerical simulations and experiments suggest that $C_\varepsilon$ may vary depending on large-scale flow characteristics (Mazellier & Vassilicos Reference Mazellier and Vassilicos2008; Goto & Vassilicos Reference Goto and Vassilicos2016; Takamure et al. Reference Takamure, Sakai, Ito, Iwano and Hayase2019). Large scales of stratified turbulence should be strongly influenced by the Froude number. For example, pancake-like large-scale vortices are often observed in stratified turbulence (Godoy-Diana et al. Reference Godoy-Diana, Chomaz and Billant2004). The dependence of large scales on $Re_M$ and $Fr_M$ can result in variations of $\alpha$. Furthermore, wind-tunnel experiments of active-grid turbulence indicate that $C_\varepsilon$ are different depending on the grid geometry, the operation mode of the active grid and the Reynolds number, although $C_\varepsilon$ in each experiment hardly varies with the decay of turbulence (Mora et al. Reference Mora, Pladellorens, Turró, Lagauzere and Obligado2019; Zheng, Nagata & Watanabe Reference Zheng, Nagata and Watanabe2021b). This behaviour is similar to the evolution of $\alpha$ in figure 15(a) since $\alpha$ in a late time weakly depends on time but differs depending on $Re_M$ and $Fr_M$. Although the underlying assumption is different for $\varepsilon = C_\varepsilon U^3/L$ in non-stratified turbulence and $\varepsilon = \alpha U_H^3/L_H$ in stratified turbulence, these previous studies of non-stratified turbulence imply that $\alpha$ does not have to be a universal constant.

In all DNS the flow reaches the viscosity-affected stratified regime during the decay and $\alpha$ hardly varies with time at a late time. Furthermore, $U_{H}^{2}L_{H}^{2}L_V$ also hardly depends on time as suggested from the theory of Saffman turbulence. In the rest of the paper, the decay properties are compared with the theoretical results of low-$Re_b$ Saffman turbulence presented in § 2. First, we briefly discuss the decay of main variables and then provide detailed discussions on the decay exponents based on least square estimation.

Figure 16(a) examines the decay of $U_H^2$, which is compared with the power law $U_H^2\sim t^{-5/4}$. The slope of $U_H^2$ is close to the $-5/4$th law except for R10F6. The decay of $U_H^2$ for R10F6 is slightly slower than $U_H^2\sim t^{-5/4}$ even for $t/t_r\geq 100$, for which the assumptions of constant values of $\alpha$ and $U_H^2L_H^2L_V$ are valid. Then, the deviation from $U_H^2\sim t^{-5/4}$ should be attributed to the validity of (2.18), which may not be valid when $Fr_M$ is not sufficiently low. It is also clear that the decay exponent of R10F6 is not constant for $30\lesssim t/t_r\lesssim 100$ because $\alpha$ also varies with time, although the grid turbulence has fully developed by $t/t_r=30$. Figure 16(b) shows the variation of the kinetic energy dissipation rate $\varepsilon$. As expected for Saffman turbulence with low $Re_b$, the present DNS results with $Fr_M\leq 1$ follow $\varepsilon \sim t^{-9/4}$ during the decay.

Figure 16.

Decay of (a) horizontal velocity fluctuations $U_H^2$ and (b) turbulent kinetic energy dissipation rate $\varepsilon$ at $y/M=0$.

Figure 17(a,b) shows the variations of integral scales, $L_H$ and $L_V$. The $Re_M$ dependence is weak for $L_H$, which is more dependent on $Fr_M$. The increase of $L_H$ follows $L_H\sim t^{3/8}$ from $t/t_r\approx 20$ except for R10F6, for which $L_H\sim t^{3/8}$ is valid only after $t/t_r\approx 100$. Although $L_V$ depends on both $Re_M$ and $Fr_M$ at an early time, it hardly depends on $Re_M$ at a late time as the lines for different $Fr_M$ converge to a single line. The DNS results with $Fr_M\leq 1$ follow $L_V\sim t^{1/2}$. However, $L_V\sim t^{1/2}$ tends to prevail at a later time than $L_H\sim t^{3/8}$ especially for higher $Fr_M$. For R10F6, the exponent is different from 1/2 even at the end of the simulation.

Figure 17.

Temporal evolutions of (a) horizontal and (b) vertical integral length scales at $y/M=0$.

Turbulence is initially generated as turbulent wakes, whose interaction results in the formation of grid turbulence. Therefore, the flow behaviour at an early time can be related to that of stably stratified turbulent wakes. For an axisymmetric wake generated by a sphere or a disk, the characteristic length of large scales is often defined with the spatial distribution of mean velocity deficit or turbulent kinetic energy (Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011; Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020). The DNS of the stably stratified wake of a disk shows that both horizontal and vertical length scales slightly decrease when the turbulent wake is being generated (Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020). Once the wake becomes turbulent, these length scales increase as the wake evolves. These evolutions of the length scales are similar in figure 17, and the wake length scales defined with mean velocity or turbulent kinetic energy can be related to the integral scales defined with the autocorrelation functions of horizontal velocity.

Figure 18(a) examines the decay of potential energy $E_P$. The present DNS results exhibit wave-like oscillations of $E_P$, as also found in previous DNS (de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019). This can be confirmed from the inset of figure 18(a), where $E_P$ between $t/t_r=105$ and 107 is shown for R10F01 in a linear plot. However, the amplitude of the oscillation is much smaller than the variation of $E_P$ during the decay, for which $E_P$ decreases from $E_P/U_0^2={{O}}(10^{-3})$ to ${{O}}(10^{-8})$. Therefore, the wave-like oscillation of $E_P$ is not visible in the logarithmic plot and does not affect the overall decay rate of $E_P$. The decay for $Fr_M=0.1$ and 0.4 tends to follow the power law of $E_P\sim t^{-7/2}$. However, a deviation from $E_P\sim t^{-7/2}$ is found for $Fr_M=0.4$ after $t/t_r\approx 150$. The power law $E_P\sim t^{-7/2}$ assumes (2.15) and (2.18), which can be used to derive $E_P/k_H \sim Re_b$ for the ratio of $E_P$ to the horizontal kinetic energy $k_H=U_H^2$. Figure 18(b) plots the relation between $Re_b$ and $E_P/k_H$ after $Re_b$ reaches a peak. The present DNS results are compared with the DNS of forced, stably stratified homogeneous turbulence in Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007). Below $Re_b\approx 10$, $E_P/k_H$ decreases with $Re_b$ being small. The ratio $E_P/k_H$ for $Fr_M=0.1$ tends to follow $E_P/k_H \sim Re_b$ for $Re_b\lesssim 10^{-1}$. Although this relation is found for $Re_b\approx 10^{-2}$ in the cases of $Fr_M=0.4$, the DNS data deviate from $E_P/k_H\sim Re_b$ at a lower $Re_b$. For R10F01, $E_P/k_H\sim Re_b$ is valid for a wide range of $Re_b$ and, therefore, the decay of $E_P$ approximately follows $E_P\sim t^{-7/2}$.

Figure 18.

(a) Decay of potential energy $E_P$ at $y/M=0$. (b) The buoyancy Reynolds number dependence of the ratio between the potential and horizontal kinetic energies $E_P/k_H$. The trajectory in $(Re_b,E_P/k_H)$ space at $y/M=0$ is shown after $Re_b$ reaches a peak. Direct numerical simulation results of forced, stably stratified homogeneous turbulence are also shown for comparison (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007).

Although the discussion in § 2 does not provide the decay law of the potential energy dissipation rate $\varepsilon _P$, the decay law in a form of $\varepsilon _P\sim t^{n}$ may be empirically established because the ratio $\varepsilon _P/\varepsilon$ often depends on $Re_b$ (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). Figure 19(a) shows the relation between $Re_b$ and $\varepsilon _P/\varepsilon$, which is compared again with DNS by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007). For $Re_b\gtrsim 10^{-1}$, the variation of $\varepsilon _P/\varepsilon$ with $Re_b$ is consistent between decaying grid turbulence and forced homogeneous turbulence. For $Re_b\lesssim 10^{-1}$, $\varepsilon _P/\varepsilon$ decreases as $Re_b$ becomes small. The present data roughly follows $\varepsilon _P/\varepsilon \sim Re_b$ for $Re_b\lesssim 10^{-2}$. This relation $\varepsilon _P\sim \varepsilon Re_b=\varepsilon ^{2}/N^2\nu$ combined with the decay laws of $\varepsilon$ in (2.19a–d) and (2.20a–d) yields

(4.8)$$\begin{gather} \varepsilon_P\sim t^{{-}9/2}\quad \text{for low-}Re_b\ \text{Saffman turbulence}, \end{gather}$$

(4.9)$$\begin{gather}\varepsilon_P\sim t^{{-}5}\quad \text{for low-}Re_b\ \text{Batchelor turbulence}. \end{gather}$$

An alternative way to derive these decay laws is neglecting the buoyancy flux term in (2.7), and ${\rm d} E_P/{\rm d} t=-\varepsilon _P$ with (2.21) and (2.22) yields (4.8) and (4.9), respectively. Indeed, the wave-like oscillation of $E_P$ arising from the buoyancy flux was shown to be smaller than the overall decay of $E_P$ in figure 18(a). Figure 19(b) compares the decay of $\varepsilon _P$ with $\varepsilon _P\sim t^{-9/2}$. The present DNS data except for R10F6 follows $\varepsilon _P\sim t^{-9/2}$ derived based on the empirical relation $\varepsilon _P\sim \varepsilon Re_b$ combined with the decay law of Saffman turbulence.

Figure 19.

(a) The buoyancy Reynolds number dependence of the ratio between the dissipation rates of potential energy and kinetic energy $\varepsilon _P/\varepsilon$. The results are shown after the peak of $Re_b$. Direct numerical simulation results of forced, stably stratified homogeneous turbulence are also shown for comparison (Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007). (b) Decay of potential energy dissipation rate $\varepsilon _P$. All results are taken at $y/M=0$.

It should be noted that $\varepsilon _P\sim \varepsilon Re_b$ is valid only when $Re_b$ is sufficiently small. Figure 19(a) suggest that $\varepsilon _P/\varepsilon$ weakly depends on $Re_b$ for $1\lesssim Re_b \lesssim 10$. Strong $Re_b$-dependence of $\varepsilon _P/\varepsilon$ is observed for $Re_b\gtrsim 10^1$. This range of $Re_b$ is in the production phase of grid turbulence and the grid turbulence has not fully developed yet. Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) also observed that $\varepsilon _P/\varepsilon =0.4$–$0.5$ hardly depends on $Re_b$ when $Re_b\sim Re_HFr_H^2>1$. Similar values of $\varepsilon _P/\varepsilon$ were also reported by Riley & de Bruyn Kops (Reference Riley and de Bruyn Kops2003), Lindborg (Reference Lindborg2006) and Maffioli & Davidson (Reference Maffioli and Davidson2016). If $\varepsilon _P/\varepsilon$ is constant during the decay, the decay exponent of $\varepsilon _P$ is the same as that of $\varepsilon$. The relation $\varepsilon _P\sim \varepsilon$ might be useful at a high-$Re_b$ regime, for which one can derive, from (2.13a–d) and (2.14a–d), the following decay laws of $\varepsilon _P$:

(4.10)$$\begin{gather} \varepsilon_P\sim t^{{-}9/5}\quad \text{for high-}Re_b\ \text{Saffman turbulence}, \end{gather}$$

(4.11)$$\begin{gather}\varepsilon_P\sim t^{{-}15/7}\quad \text{for high-}Re_b \text{Batchelor turbulence}. \end{gather}$$

However, the Reynolds number is not large enough in the present DNS to verify these decay laws.

One of the important parameters in oceanography is mixing efficiency that is often defined as $\eta _{mix}=\varepsilon _P/(\varepsilon +\varepsilon _P)=(1+ \varepsilon /\varepsilon _P)^{-1}$ (Smyth, Moum & Caldwell Reference Smyth, Moum and Caldwell2001; Gregg et al. Reference Gregg, D'Asaro, Riley and Kunze2018). The DNS of stably stratified turbulence has been utilized to parametrize $\eta _{mix}$ (Salehipour, Peltier & Mashayek Reference Salehipour, Peltier and Mashayek2015; Maffioli, Brethouwer & Lindborg Reference Maffioli, Brethouwer and Lindborg2016; de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019; VanDine, Pham & Sarkar Reference VanDine, Pham and Sarkar2021). If $\varepsilon _P/\varepsilon$ is constant during the decay, $\eta _{mix}$ does not vary with time. However, the decay laws derived for low-$Re_b$ Saffman and Bachelor turbulence suggest that $\eta _{mix}$ decays with time as $\eta _{mix}\sim (1+A_\eta t^{9/4})^{-1}$ for Saffman turbulence and $\eta _{mix}\sim (1+A_\eta t^{5/2})^{-1}$ for Batchelor turbulence. As long as the late time behaviour of grid turbulence is concerned, we may approximate these decay laws as

(4.12)$$\begin{gather} \eta_{mix} \sim t^{{-}9/4}\quad \text{for low-}Re_b\ \text{Saffman turbulence}, \end{gather}$$

(4.13)$$\begin{gather}\eta_{mix} \sim t^{{-}5/2}\quad \text{for low-}Re_b\ \text{Batchelor turbulence}. \end{gather}$$

Figure 20 shows the temporal variation of $\eta _{mix}$ compared with (4.12). As expected, $\eta _{mix}$ decreases as turbulence decays. Before the decay of $\eta _{mix}$ begins, $\eta _{mix}$ is about $0.33$, which is close to the values corresponding to $\varepsilon _P/\varepsilon =0.4$–$0.5$ and is consistent with other DNS studies of stratified homogeneous turbulence (de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019). The decay of $\eta _{mix}$ tends to follow (4.12), although $\eta _{mix} \sim t^{-9/4}$ is observed at a later time for higher $Fr_M$. de Bruyn Kops & Riley (Reference de Bruyn Kops and Riley2019) also found in DNS of decaying, stably stratified homogeneous turbulence that $\eta _{mix}$ stays at an almost constant value at an early time of the decay with high $Re_b$ and then decreases with time once $Re_b$ becomes small, although this decay was not emphasized in their paper.

Figure 20.

Temporal variation of mixing efficiency $\eta _{mix}=\varepsilon _P/(\varepsilon +\varepsilon _P)$ at $y/M=0$. Horizontal broken lines represent $\eta _{mix}$ with $\varepsilon _P/\varepsilon =0.4$ and 0.5.

Figure 21 shows the decay of vertical velocity variance $\langle w'^2\rangle$. The theories for stably stratified turbulence in § 2 concerns large-scale characteristics of the flow and do not provide the decay law of $\langle w'^2\rangle$, which is dominated by small-scale fluctuations in stratified turbulence. However, $\langle w'^2\rangle$ tends to follow $\langle w'^2\rangle \sim t^{-1}$ except for R10F6. Towed-grid experiments by Praud et al. (Reference Praud, Fincham and Sommeria2005) also observed $\langle w'^2\rangle \sim t^{-1}$ at low $Fr_M$.

Figure 21.

Decay of vertical velocity variance $\langle w'^2\rangle$ at $y/M=0$.

These results suggest that temporally evolving grid turbulence in a stratified fluid behaves similarly to Saffman turbulence. Previous numerical studies of Saffman turbulence and Batchelor turbulence consider decaying turbulence initialized with a velocity field that satisfies $E(k)\sim k^2$ or $k^4$ at large scales (Ishida, Davidson & Kaneda Reference Ishida, Davidson and Kaneda2006; Teitelbaum & Mininni Reference Teitelbaum and Mininni2012; Yoshimatsu & Kaneda Reference Yoshimatsu and Kaneda2018; de Bruyn Kops & Riley Reference de Bruyn Kops and Riley2019; Anas et al. Reference Anas, Joshi and Verma2020). It has been shown that the decay of turbulence with the $k^2$ and $k^4$ spectra is consistent with the theories of Saffman and Batchelor turbulence, respectively. However, these studies do not reveal how turbulence with $E(k)\sim k^2$ or $k^4$ is generated because the energy spectrum is imposed as the initial condition. The temporally evolving grid turbulence is initialized with a laminar flow that approximates the mean velocity deficit of the wakes. Therefore, the production of turbulence occurs simply due to the mean velocity gradient imposed as the initial condition. The results presented above imply that the mean velocity deficit caused by the grid is essential for grid turbulence to behave like Saffman turbulence.

4.3. Estimation of decay exponents

The results shown above indicate that most quantities vary with time in accord with a power law of $t$. From experiments and numerical simulations of grid turbulence with and without stable stratification, it is widely accepted that the decay of horizontal velocity fluctuations follows

(4.14)\begin{equation} \frac{U_H^2}{U_0^2}=A\left(\frac{t}{t_r}-\frac{t_0}{t_r}\right)^{n}, \end{equation}

where $t_0$ is the virtual origin. Various methods have been proposed to determine the unknown parameters $A$, $t_0$ and $n$ from data taken in laboratory experiments of non-stratified grid turbulence (Lavoie et al. Reference Lavoie, Djenidi and Antonia2007; Krogstad & Davidson Reference Krogstad and Davidson2010; Valente & Vassilicos Reference Valente and Vassilicos2011). Most methods provide a similar decay exponent $n$, although small variations may exist because the virtual origin is determined in a different manner. These methods are based on linear or nonlinear least square estimation applied to a certain range of decay. It is known that the fitting range $t_{s} \leq t \leq t_{e}$ may significantly affect $n$ in non-stratified grid turbulence, where $t_{s}$ and $t_{e}$ are the start and end points (time) of the fitting range. As discussed in the introduction, a near-grid region has a faster decay rate than a fully developed turbulence region (Komori et al. Reference Komori, Nagata, Kanzaki and Murakami1993; Krogstad & Davidson Reference Krogstad and Davidson2012). Furthermore, $n$ also increases as turbulence decays to a low-Reynolds-number regime. Here, the variation of $n$ is often explained by that of $C_\varepsilon$. Except for R10Fr6 the temporal variation of $\alpha =\varepsilon /(U_H^3/L_H)$ in figure 15 is large for $t/t_r\lesssim 20$, which may correspond to the non-equilibrium decay region in the case of non-stratified grid turbulence generated by a square grid (Isaza, Salazar & Warhaft Reference Isaza, Salazar and Warhaft2014; Nagata et al. Reference Nagata, Saiki, Sakai, Ito and Iwano2017). Therefore, $t/t_r\lesssim 20$ should be eliminated from the data used to estimate $n$. Furthermore, the power law of $L_V$ predicted for low $Re_b$ seems valid for $t/t_r\gtrsim 20$ for $Fr=0.1$, $t/t_r\gtrsim 30$ for $Fr=0.4$ and $t/t_r\gtrsim 60$ for $Fr=1$ in figure 17(b). By taking these points into consideration, $t_{s}/t_r=70$ is adapted except for R10F6. We have also confirmed that the decay exponents obtained with larger $t_{s}$ hardly differ from those with $t_s/t_r=70$. The assumption of constant $\alpha$ becomes valid only after $t\approx 100$ for R10F6, and $n$ is estimated with $t_s/t_r=120$. In all cases, the fitting range is within the viscosity-affected stratified flow regime. Experiments at very large $t/t_r$ $(x/M)$ often have a difficulty in dealing with background velocity fluctuations that can be comparable to turbulent fluctuations and with measurement errors arising from, for example, electronic noise of hot-wire anemometry. Therefore, the end point of the fitting range has to be carefully determined in experiments. However, these issues do not affect the estimation of $n$ in numerical simulations. Furthermore, the decay of turbulence in a stably stratified fluid does not result in the increase of $\alpha$, which is almost constant even in a late time. The estimation of $n$ in stratified grid turbulence is not as sensitive to $t_e$ as in experiments of non-stratified grid turbulence. Therefore, $t_e$ is set as the end of the computations, i.e. $t_e/t_r=250$.

This study compares $n$ obtained by three methods explained here. Method 1 predetermines $t_0/t_r$ and applies a linear least square method to the DNS data of $(U_H^2/U_0^2, t/t_r-t_0/t_r)$ in a form of $\ln (U_H^2/U_0^2)=\ln A+ n\ln (t/t_r-t_0/t_r)$, for which $\ln A$ and $n$ are analytically determined. The fitting is repeatedly applied by changing $t_0/t_r$ from $-100$ to $100$ with an increment of $0.1$. The virtual origin $t_0/t_r$ is determined such that the standard deviation of the DNS data from the fitting function is the smallest for all tested virtual origins. This virtual origin is used to determine $\ln A$ and $n$ with the linear least square method. The condition of the minimum standard deviation is often used to determine one of the fitting parameters, e.g. a range of $t$ (or $x$ in wind-tunnel experiments) where the least square method is applied to calculate $n$ (Krogstad & Davidson Reference Krogstad and Davidson2010; Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014). Method 2 assumes that the virtual origin is equal to zero and the linear least square method is applied for $\ln (U_H^2/U_0^2)=\ln A+ n\ln (t/t_r)$ to determine $\ln A$ and $n$. This method yields the decay exponent equivalent to those considered in some previous studies on stably stratified grid turbulence, where the decay exponent was obtained for $(U_H^2/U_0^2)=A(t/t_r)^{n}$ with the assumption of $t_0=0$ (Liu Reference Liu1995; Praud et al. Reference Praud, Fincham and Sommeria2005; Espa et al. Reference Espa, Avallone and Cenedese2018). Method 3 directly applies a nonlinear least square method (Levenberg–Marquardt method) to simultaneously determine $A$, $t_0/t_r$ and $n$ of (4.14). The nonlinear least square method has been widely used for experimental data of grid turbulence and was shown to provide $n$ consistent with other methods (Valente & Vassilicos Reference Valente and Vassilicos2011; Isaza et al. Reference Isaza, Salazar and Warhaft2014; Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014; Zheng, Nagata & Watanabe Reference Zheng, Nagata and Watanabe2021a).

Table 2 summarizes the decay exponent and virtual origin obtained by the three methods. Exponent $n$ for R10F6 is smaller than for the other cases, as also expected from the slower decay of R10F6 in figure 16(a). However, $n$ for $Fr=0.1$, 0.4 and 1 is close to $n=-1.25$ derived for low-$Re_b$ Saffman turbulence. Although each method yields different virtual origins, the estimated values of $n$ differ only slightly. This weak sensitivity of $n$ to the choice of $t_0$ is related to the fitting range that includes large $t/t_r$ up to 250 because the decay exponent estimated from data at large $t/t_r$ is known to be insensitive to small variations of the virtual origin (Zheng et al. Reference Zheng, Nagata and Watanabe2021a). The variations of $n$ found among the three methods are at a similar level to experimental results of non-stratified grid turbulence (Valente & Vassilicos Reference Valente and Vassilicos2011). Averages of $n$ taken except for R10F6 are $-1.20$, $-1.25$ and $-1.23$ for methods 1, 2 and 3, respectively. These values are fairly close to the prediction $n=-1.25$ for low-$Re_b$ Saffman turbulence.

Table 2.

Decay exponents and virtual origins estimated by three different methods.

The DNS results suggest that a power law $f=A(t-t_0)^{n}$ is also valid for various quantities shown in figures 16–19. Exponent $n$ is estimated by applying the linear least square method to $\ln f=\ln A+ n\ln (t-t_0)$ with the virtual origin estimated from $U_H^2$ by method 1. Table 3 summarizes the estimated values of $n$ for the integral scales, potential energy and dissipation rates of kinetic and potential energies. The exponents for $L_H$, $L_V$ and $\varepsilon$ are close to the values for low-$Re_b$ Saffman turbulence except for R10F6. The faster growth of $L_V$ than $L_H$ is also found for DNS data as predicted for low-$Re_b$ Saffman turbulence. The decay of $E_P$ and $\varepsilon _P$ of the DNS data tends to be slower than low-$Re_b$ Saffman turbulence, although the difference from the theory is not large for some cases, e.g. R10F04 and R10F01. This is because the exponents of $E_P/k_H \sim Re_b^{p}$ and $\varepsilon _P/\varepsilon \sim Re_b^{p}$ deviate from $p=1$ especially for large $Re_b$ while $n=-3.5$ for $E_P$ and $-4.5$ for $\varepsilon _P$ require $p=1$.

Table 3.

Power-law exponents of integral scales, potential energy and dissipation rates of kinetic and potential energies.