Direct numerical simulation of temporally evolving stratified wakes with ensemble average

Abstract

Direct numerical simulations are conducted for temporally evolving stratified wake flows at Reynolds numbers from $10\,000$ to $50\,000$ and Froude numbers from $2$ to 50. Unlike previous studies that obtained statistics from a single realization, we take ensemble averages among 80–100 realizations. Our analysis shows that data from one realization incur large convergence errors. These errors reduce quickly as the number of statistical samples increases, with the benefit of ensemble average diminishing beyond 40–60 realizations. The data with ensemble average allow us to test the previously established scalings and arrive at new scaling estimates. Specifically, the data do not support power-law scaling in the centreline velocity deficit $U_0$ beyond the near wake. Its decay rate increases continuously from 0.1 at the onset of the non-equilibrium regime until the end of our calculations without reaching any asymptote. Additionally, while no power-law scalings could be found in the wake width ($L_H$) and wake height ($L_V$) in the late wake, $L_H\sim (Nt)^{1/3}$ is a good working approximation of the wake's horizontal size, where $N$ is buoyancy frequency and $t$ is time. Besides the low-order statistics, we also report the transverse integrated terms and the vertically integrated terms in the turbulent kinetic energy budget equation as a function of the vertical and transverse coordinates. The data indicate that there are two peaks in the vertically integrated production and transport terms, and one peak when the two terms are integrated horizontally.

1. Introduction

The evolution of a wake in a stratified environment has garnered much attention since the early work by Lin & Pao (1979), and reviews of the literature can be found in Spedding (2014) and de Bruyn Kops & Riley (2019). Two model problems that have been studied extensively numerically are temporally evolving wakes and spatially developing wakes in a linearly stratified environment, which are sketched in figures 1 and 2. The spatial coordinate $x$ in figure 1(a) and the time coordinate $t$ in figure 1(b) are interchangeable according to $x=U_bt$, where $U_b$ is the convective velocity. Simulations of spatially evolving wakes necessitate the concurrent resolution of both the near wake and the far wake. On the other hand, simulations of temporally evolving wakes require only the wake characterization at a specific time. Consequently, when investigating the late wake, simulations of temporally evolving wakes are often more straightforward than simulations of spatially developing wakes (Chongsiripinyo & Sarkar 2020). Nevertheless, the initialization of temporally evolving wake simulations remains challenging practically, unless calculations of a spatially evolving wake at the same flow condition are readily available (Pasquetti 2011; VanDine, Chongsiripinyo & Sarkar 2018).

Figure 1. Schematics of stratified wake flows. (a) A spatially developing wake in (i) the streamwise vertical plane and (ii) the streamwise transverse plane. (b) A temporally evolving wake in (i) the streamwise vertical plane and (ii) the streamwise transverse plane. The fields are from case R20F02 at $Nt=4.5$, 52, 110 and 400, which will be detailed in § 2.

Figure 2. (a) Visualization of the instantaneous vertical vorticity field $\omega _z$ on the horizontal plane $z=0$. (b) Instantaneous spanwise vorticity field $\omega _y$ on the vertical plane $y=0$. Four snapshots are at $Nt=4.5$, 52, 110 and 400, respectively. The fields are from case R20F02.

Significant insights have been gained into the flow physics of stratified wakes since the seminal work of Lin & Pao (1979). Here, we provide a brief overview of the basic flow phenomenology. First, we define the local Froude number $Fr=U_0/(NL)$, where $L$ is an integral length scale, $U_0$ is the centreline velocity deficit, and $N$ is the Brunt–Väisälä frequency. The local Froude number describes the competition between the inertial and buoyancy forces. In many real-world scenarios, such as the wake behind an underwater vehicle in the deep sea, Froude numbers start off large, gradually decrease as the wake decays, and eventually reach order one. From here, buoyancy becomes increasingly important, and the wake loses its memory of the wake-generating body (Meunier & Spedding 2004). Following Spedding (1997), the wake can be divided into three distinct flow regimes: the near-wake regime (NW), the non-equilibrium regime (NEQ), and the quasi-two-dimensional regime (Q2D). Buoyancy effects are negligible in the NW regime, but are important in the NEQ and Q2D regimes. As different physics are at play in the different flow regimes, important flow quantities like the centreline velocity deficit and the wake width and height exhibit different scaling behaviours in these regimes. In the following, we provide a review of existing data, with a focus on direct numerical simulations (DNS) of temporally evolving wakes. Experimental data will be referenced as needed.

The early experiments conducted by Spedding, Browand & Fincham (1996) and Spedding (1997) indicated that the centreline velocity deficit decays as $x^{-2/3}$, $x^{-0.25}$ and $x^{-0.76}$ in the NW, NEQ and Q2D regimes, respectively. Spedding (1997) proposed defining the flow regimes based on the scaling behaviour of the centreline velocity deficit $U_0$. According to that study, the NW regime ends at approximately $Nt=1.7\pm 0.3$, and NEQ ends at approximately $Nt=50\pm 15$, beyond which the flow transitions to the Q2D regime. The above power-law decay of the centreline velocity deficit has also been observed in other studies, such as Diamessis, Domaradzki & Hesthaven (2005), Brucker & Sarkar (2010) and Diamessis, Spedding & Domaradzki (2011), among others. However, the specific time or downstream location at which the NW and NEQ regimes end varies across different investigations. For example, Brucker & Sarkar (2010) identified the termination of the NW regime at $Nt=5$ and the NEQ regime at $Nt=100$. In contrast, Diamessis et al. (2011) and Spedding (1997) reported an early transition to the Q2D regime at $Nt=60$ and $Nt=50$, respectively.

Like the extent of the three flow regimes, the scalings of the centreline velocity deficit and wake sizes in each regime remain open research questions. Here, we summarize the current understanding of these scalings in the NW, NEQ and Q2D regimes. In the NW regime, the behaviour of a stratified wake is similar to that of an unstratified wake. The centreline velocity deficit decays as $U_0\sim x^{-2/3}$, while the vertical and horizontal sizes of the wake grow as $L_v\approx L_h\sim x^{1/3}$ (Tennekes & Lumley 1972; Townsend 1976). Some studies have observed a $x^{-1}$ scaling in $U_0$, which arises from the anisotropic dissipation of turbulent kinetic energy (TKE) (Nedić, Vassilicos & Ganapathisubramani 2013; Dairay, Obligado & Vassilicos 2015; Pal et al. 2017; Chongsiripinyo & Sarkar 2020; Ortiz-Tarin, Nidhan & Sarkar 2021).

Moving on to the NEQ regime, buoyancy becomes significant. The streamwise and vertical velocity fluctuations $u'$ and $w'$ are less correlated compared to the NW regime, resulting in reduced energy transfer from the mean flow to turbulence. As a consequence, the decay of the centreline velocity deficit in the NEQ regime is slower than in the NW regime (Brucker & Sarkar 2010). Meunier, Diamessis & Spedding (2006) argued that the NEQ regime can be divided further into an early NEQ regime and a late NEQ regime. In the early NEQ regime, the centreline velocity deficit decays as $x^{-0.24}$, while in the late NEQ regime, it decays as $x^{-0.37}$ (Brucker & Sarkar 2010). Regarding wake sizes, $L_V$ remains relatively constant in the NEQ regime, while $L_H$ shows approximately $x^{1/3}$ growth (Spedding 1997; Gourlay et al. 2001; Diamessis et al. 2005, 2011; Brucker & Sarkar 2010).

Finally, in the Q2D regime, as we can see in figure 2, large-scale pancake-like eddies emerge (Chongsiripinyo & Sarkar 2020). The centreline velocity deficit decays as $U_0\sim x^{-0.76}$ (Spedding 1997; Gourlay et al. 2001; Diamessis et al. 2005, 2011; Brucker & Sarkar 2010). The scalings of wake sizes in this regime are subject to some controversy. Spedding (1997), Gourlay et al. (2001), Dommermuth et al. (2002) and Diamessis et al. (2011) reported approximately $x^{1/3}$ growth in $L_H$. However, Brucker & Sarkar (2010) found a varying growth rate of the wake's horizontal size, with $L_H\sim x^{0.46}$ and $L_H\sim x^{0.23}$ for $100< Nt<250$ and $Nt>350$, respectively. As for the vertical size, Gourlay et al. (2001), Diamessis et al. (2005), Brucker & Sarkar (2010) and Diamessis et al. (2011) reported $x^{0.32\pm 0.1}$, $x^{0.6}$, $x^{0.5}$ and $x^{0.36-0.48}$, respectively. Tables 1–3 summarize the data. In summary, we see that the available data are scattered and exhibit variation across different studies.

Table 1. Decay rate of the centreline velocity deficit. Here, S97 is Spedding (1997), G01 is Gourlay et al. (2001), D02 is Dommermuth et al. (2002), BS10 is Brucker & Sarkar (2010) and D11 is Diamessis et al. (2011). Diamessis et al. (2005) (D05) did not report the decay rate, leading to the NA in the table. The values of the Froude number are adjusted so that they are consistent with our definition.

Table 2. Growth rate of $L_h$. The flow parameters and the references are the same as in table 1. D02 did not report the exact growth rate, leading to the NA in the table.

Table 3. Growth rate of the wake's vertical size $L_v$. The flow parameters and the references are the same as in table 1. In BS10, there is no numerical number reported, leading to NA in the table. A caveat is that one may get vastly different values depending on the definition of $L_v$.

Next, we try to understand the scattering in the data. In figure 3, we have collected some of the centreline velocity deficit data in the literature (Dommermuth et al. 2002; Diamessis et al. 2005, 2011; Brucker & Sarkar 2010). The data cover the Reynolds number range between 5000 and 400 000, and the Froude number range between 2 and 50. We note that the Froude number in Diamessis et al. (2011) is defined as $Fr \equiv 2U/ND$, which is different from the definition here, $Fr=U/ND$. The numbers in Diamessis et al. (2011) are adjusted so that they are consistent with the definition here. Our data are also included in the figure for comparison purposes. The details of our DNS will be presented in § 2. We observe the following. First, wiggles can be observed in the data, particularly in the late NEQ and Q2D regimes. These fluctuations make it challenging to discern clear scaling trends for $U_0$. Second, it is difficult to identify a specific range where the centreline velocity deficit follows a consistent $-0.25$ or $-0.76$ scaling. In fact, the decay rate does not appear to remain constant over an extended range. These factors contribute to the scattering of data in table 1. In addition, a lack of good theories and clear definitions is also responsible. In comparison to flow statistics in isotropic turbulence and boundary-layer flows, where Kolmogorov's theory of small-scale turbulence (Frisch 1995) and Townsend's attached eddy hypothesis (Marusic & Monty 2019; Yang & Meneveau 2019) provide good baseline estimates, flow statistics in stratified wake flows lack well-established baseline estimates. Meunier et al. (2006) provide estimates of the scalings of the mean flow statistics in a stratified wake, but the work is much less well-established compared to Kolmogorov's and Townsend's theories.

Figure 3. A not comprehensive collection of the data in the literature showing the decay of the centreline velocity deficit. Flow parameters are identified in the legend as R[$Re/1000$]F[$Fr$]. Here, solid lines are our DNS (detailed in § 2); D11 is by Diamessis et al. (2011); D05 is by Diamessis et al. (2005); D02 is by Dommermuth et al. (2002); BS10 is by Brucker & Sarkar (2010). The values of the Froude number are adjusted so that they are consistent with our definition.

Furthermore, there exist different definitions of wake width and height. They may be defined based on the velocity (Spedding 1997, 2001), the TKE (Pal et al. 2017; Ortiz-Tarin et al. 2021) or other statistical objects (Brucker & Sarkar 2010; de Stadler, Sarkar & Brucker 2010). The behaviours of the heights and widths depend on their definitions, which also contribute to the scattering of data in tables 2 and 3. Here, we distinguish $L_H$ and $L_V$ from the integral length scales in Zhou & Diamessis (2019) and de Bruyn Kops & Riley (2019). The integral length scale in Zhou & Diamessis (2019) characterizes the scale of the underlying turbulence, and the integral length scale in de Bruyn Kops & Riley (2019) characterizes the scale of the mean flow.

The discussion above focuses on the classification of stratified wakes into the NW, NEQ and Q2D regimes. Alternatively, we may divide stratified wakes into weakly stratified turbulence (WST), intermediately stratified turbulence (IST), strongly stratified turbulence (SST) and viscous regimes (Billant & Chomaz 2001; Riley & de Bruyn Kops 2003; Lindborg 2006; Brethouwer et al. 2007; Chongsiripinyo & Sarkar 2020). This classification is based on the local horizontal Reynolds number $Re_h$ and the local horizontal Froude number $Fr_h$, which are defined as

(1.1a,b)\begin{equation} Re_h = \frac{u'_hL_{Hk}}{\nu}, \quad Fr_h=\frac{u'_h}{NL_{Hk}}, \end{equation}

where $u'_h$ is the root mean square of the horizontal velocity fluctuation $u'_h=(u'^2+v'^2)^{1/2}$, and $L_{Hk}$ is the distance from the centreline to where the TKE is half of its centreline value in the transverse direction. The parameter $Fr_h$ characterizes the importance of buoyancy, while $Re_h\,Fr_h^2$ characterizes the importance of inertia. Figure 4 illustrates the $Fr_h$–$Re_h\,Fr_h^2$ phase space and the evolution of stratified wakes within it. The WST regime occupies the top right portion of the phase space, where $0.1< Fr_h<1$. Below the WST regime lie the IST and SST regimes, characterized by $0.03< Fr_h<0.1$ and $Fr_h<0.03$, respectively (Chongsiripinyo & Sarkar 2020). The viscous regime is located to the left of the stratified turbulence regimes and is defined by $Re_h\,Fr_h^2<1$. As a stratified wake at high Reynolds and Froude numbers evolves, it passes through the WST regime, IST regime, SST regime, and eventually the viscous regime. Besides $Fr_h$–$Re_h\,Fr_h^2$, another relevant space is the $Fr_h$–$\epsilon /(\nu N^2)$ space. Figure 5 shows the evolution of stratified wakes in this space. Nonetheless, since the statistic $\epsilon /(\nu N^2)$ is rarely reported in DNS studies of stratified wake flows, we could show our data only. The $Fr_h$–$\epsilon /(\nu N^2)$ space is not very different from the $Fr_h$–$Re_h\,Fr_h^2$ space in terms of the positions of the WST, IST, SST and viscous regimes therein. The dissipation in our DNS first increases before it decreases, as in Brucker & Sarkar (2010) and Chongsiripinyo & Sarkar (2020). This gives rise to the arc at the initial stage. The trajectories of the DNS in the $Fr_h$–$\epsilon /(\nu N^2)$ space is otherwise similar to that in the $Fr_h$–$Re_h\,Fr_h^2$ space. A detailed discussion of the two spaces and their physical bearings falls out of the scope of this introduction, and the reader is directed to Riley & de Bruyn Kops (2003) and de Bruyn Kops & Riley (2019) for more detailed discussions. Regime classifications in figures 4 and 5 remove the uncertainties in the spans of the flow regimes. However, the question regarding the scaling of the centreline velocity deficit, width and height of wakes in these stratified turbulence regimes remains open.

Figure 4. A not comprehensive collection of the data in the literature showing the phase diagram of $Re_h\,Fr_h^2$–$Fr_h$. Nomenclature is the same as in figure 3. Here, ZD19 is the LES by Zhou & Diamessis (2019), and CS19 is the LES by Chongsiripinyo & Sarkar (2020). The values of the Froude number are adjusted so that they are consistent with our definition.

Figure 5. The phase diagram of $\epsilon /(\nu N^2)$–$Fr_h$ calculated from the wake-centre value. The legend is the same as in figure 3.

A limitation of the present large-eddy simulations (LES) and DNS data is the lack of statistical samples. A calculation of a temporally evolving wake provides one statistical sample. While one can get reasonably well-converged low-order statistics by averaging in the streamwise direction in one realization, obtaining statistically converged high-order statistics such as the TKE or the terms in the TKE budget equation is challenging with only one sample. In this study, we aim to address this limitation by collecting a larger number of statistical samples at a given simulation parameter. By doing so, we can achieve a better-converged centreline velocity deficit and reasonably well-converged terms in the TKE budget equation. These additional samples will enable us to answer two outstanding questions: the scaling of the centreline velocity deficit and the impact of buoyancy on the balance of terms in the TKE budget equation. It is worth noting that the issue of statistical convergence exists in experimental studies as well (Spedding et al. 1996; Spedding 1997, 2002). In those cases, the concern stems from the limited observational window and insufficient sampling of large-scale structures.

In anticipation of the discussion in § 3, here we review the model in Meunier et al. (2006). In the late wake, the length scales of the mean flow are large in the streamwise direction and small in the vertical and transverse directions; the vertical motions are suppressed. Furthermore, the wake's convective velocity $U_b$ is much larger than the wake's centreline velocity deficit $U_0$. Thus turbulence is assumed to act in the horizontal direction only, and vertical growth is dominated by viscous diffusion. These assumptions give rise to the following simplified form for the streamwise momentum equation:

(1.2)\begin{equation} U_b\,\frac{\partial U}{\partial x}=(\nu+\nu_t)\,\frac{\partial^2 U}{\partial y^2}+\nu\,\frac{\partial^2 U}{\partial z^2}. \end{equation}

Here, the eddy viscosity assumption is invoked, resulting in $-\left \langle u'v'\right \rangle =\nu _t\,\partial U/\partial y$. In Meunier et al. (2006), the eddy viscosity was assumed to be a constant in the wake region. Meunier et al. (2006) argued that the wake profile is approximately Gaussian:

(1.3)\begin{equation} U=U_b-U_0\exp\bigg(-\frac{y^2}{L_y^2}-\frac{z^2}{2L_z^2}\bigg). \end{equation}

Equations (1.2) and (1.3) give estimates for $L_z$, $L_y$ and $U_0$ for the late wake. These estimates can be found in equations (20), (22), and (24) in Meunier et al. (2006), and are not repeated here for brevity.

The remaining sections of the paper are organized as follows. We provide the details of our DNS in § 2. The results are then presented in § 3. We first discuss the number of samples needed to get converged statistics. Subsequently, we analyse the scalings of the centreline velocity deficit, the velocity fluctuations, and the horizontal and vertical sizes of the wake. Finally, we examine the transverse and vertically integrated TKE budget at both high and low Froude numbers. The paper concludes in § 4 with a summary of the findings.

2. Direct numerical simulations

We conduct DNS of temporally evolving wakes. The DNS set-up is similar to that in Brucker & Sarkar (2010) and de Stadler et al. (2010), but we repeat our DNS 80–100 times for ensemble average.

2.1. Governing equations and normalization

The governing equations are the three-dimensional incompressible Navier–Strokes equations with the Boussinesq approximation. In their dimensional forms, the equations are

(2.1)\begin{gather} \frac{\partial{u_{i}^*}}{\partial{x_{i}^*}} = 0, \end{gather}

(2.2)\begin{gather} \frac{\partial{u_{i}^*}}{\partial{t^*}}+\frac{\partial{(u_{j}^*u_{i}^*)}}{\partial{x_{j}^*}} ={-}\frac{1}{\rho_{0}}\,\frac{\partial{{p'^*}}}{\partial{x_{i}^*}}+\nu\,\frac{\partial^2{u_{i}^*}}{\partial{x_{j}^*}\,\partial{x_{j}^*}}-\frac{\rho'^*}{\rho_{0}}\,g\delta_{i3}, \end{gather}

(2.3)\begin{gather} \frac{\partial{\rho^*}}{\partial{t^*}}+\frac{\partial{(u_{i}^*\rho^*)}}{\partial{x_{i}^*}} = \alpha\,\frac{\partial^2{\rho^*}}{\partial{x_{i}^*}\,\partial{x_{i}^*}}, \end{gather}

where (2.1) is the dimensional continuity equation, (2.2) is the dimensional momentum equation, and (2.3) is the dimensional density equation; $u_i^*$ is the dimensional velocity in the $i$th Cartesian direction, and $\rho '^{*}$ is the instantaneous density fluctuation from the background mean field. Here, instead of solving the density equation directly, we solve the transport equation for the temperature,

(2.4)\begin{equation} \frac{\partial{T^*}}{\partial{t^*}}+u_i^*\,\frac{\partial{(\bar{T}^*+T^*})}{\partial{x_{i}^*}}=\kappa\,\frac{\partial^2{T^*}}{\partial{x_{i}^*}\,\partial{x_{i}^*}}, \end{equation}

and then compute the density field using the equation

(2.5)\begin{equation} \frac{\rho^*-\rho_0}{\rho_0}={-}\beta(\bar{T}^*+T^*-T_0),\quad \beta={-}\frac{1}{\rho_0}\left(\frac{\partial{\rho}}{\partial{T}}\right)_P, \end{equation}

where $\bar {T}^*$ is the background temperature, $T^*$ is the dimensional temperature fluctuation, $\rho _0$ is the reference density, and $\beta$ is the thermal expansion coefficient.

We use the convective velocity, the reference temperature $T_0$, the reference density $\rho _0$, and the initial wake diameter $D$ for non-dimensionalization. The non-dimensional flow variables are

(2.6a–f)\begin{equation} t=\frac{t^*U_b}{D},\quad x_{i}=\frac{x_{i}^*}{D},\quad u_{i}=\frac{u_{i}^*}{U_b},\quad T=\frac{T^*}{T_0},\quad \rho =\frac{\rho^*}{\rho_0}, \quad p=\frac{{p'^*}}{\rho_{0}U_b^2}. \end{equation}

The non-dimensional forms of the governing equations are

(2.7)\begin{gather} \frac{\partial{u_{i}}}{\partial{x_{i}}}=0, \end{gather}

(2.8)\begin{gather} \frac{\partial{u_{i}}}{\partial{t}}+\frac{\partial{(u_{j}u_{i})}}{\partial{x_{j}}}={-}\frac{\partial{p}}{\partial{x_{i}}}+\frac{1}{Re}\,\frac{\partial^2{u_{i}}}{\partial{x_{j}}\,\partial{x_{j}}}-\frac{1}{Fr^2}\,{\rho'}\delta_{i3}, \end{gather}

(2.9)\begin{gather} \frac{\partial{T}}{\partial{t}}+u_{i}\,\frac{\partial{(\bar{T}+T)}}{\partial{x_{i}}}=\frac{1}{Re\,Pr}\,\frac{\partial^2{T}}{\partial{x_{i}}\,\partial{x_{i}}}, \end{gather}

where $Re = U_bD/\nu$ is the (initial) Reynolds number, $Pr = \nu /\kappa$ is the Prandtl number, $Fr = U_b/(ND)$ is the Froude number, and $\bar {T}$ is the imposed background temperature. Stable stratification is imposed through $\partial {\bar {T}}/{\partial {t}}=0$, $\partial {\bar {T}}/{\partial {z}}= {\rm const.}$ and $\partial ^2{\bar {T}}/{\partial {z}^2}=0$. From here on, all variables are non-dimensional unless noted otherwise.

2.2. Initial conditions

The initial velocity field is decomposed into a mean velocity field $\langle {u}\rangle$ and a fluctuating field $u'$. The mean field is given by

(2.10)\begin{equation} \langle u_{1} \rangle(r)=U_{0}\exp\Big({-\frac{1}{2}\left({r}/{r_0}\right)^2}\Big), \end{equation}

where $U_{0}$ here is the initial centreline velocity deficit, $r=\sqrt {x_{2}^2+x_{3}^2}$ is the radial distance from the wake centre, and $r_0$ is the initial wake radius. Here, $U_0 = 0.11$ and $r_0 = 1/2$, following Dommermuth et al. (2002) and Brucker & Sarkar (2010). The fluctuating field is generated as follows. First, a random field is generated from a uniform distribution between 0 and 1. Different realizations at given $Re$ and $Fr$ differ in these randomly generated initial fluctuation fields. Second, the field is made divergence-free and digitally filtered to match the energy spectrum

(2.11)\begin{equation} E(k) \sim (k/k_0)^4\,{\rm e}^{{-}2(k/k_0)^2}, \end{equation}

where $k_0 = 4$ (de Stadler et al. 2010). The initial turbulence intensity is denoted as $\alpha$ and is set to $0.055$. Third, the fluctuations are damped exponentially away from the centreline according to the damping function

(2.12)\begin{equation} g(r)=\alpha\Big(1+\frac{r}{r_0}\Big)^2\exp\Big(-\frac{1}{2}(r/r_0)^2\Big). \end{equation}

Fourth, we hold the mean field and allow the fluctuations to evolve until the cross-stream velocity fluctuation correlation is $0.15 < \langle {u_{1}'u_{r}'} \rangle /K < 0.25$, where $K$ is the TKE, $u_1'$ is the root mean square of streamwise velocity fluctuation, and $u_r'$ is the root mean square of radial velocity fluctuation, $u_r'=\sqrt {u_2'^2+u_3'^2}$. This period lasts for approximately $t=7$. During the initial period, the gravitational force ($1/Fr^2$) is turned off for $0< t<6$ and increases from 0 to the correct value in $6< t<8$ according to the ramping function

(2.13)\begin{equation} g(t)=g_{0}\tanh(t-t_{i}/t_{ramp}), \end{equation}

where $t_{i}=6$ and $t_{ramp}=2$. The purpose of gravitational ramping is to compensate for the lack of temperature perturbation in the initial field. The set-up described here is based on the rationale presented in Brucker & Sarkar (2007). For more detailed information on the initial condition, we refer the reader to Bevilaqua & Lykoudis (1978), Uberoi & Freymuth (1970), Dommermuth et al. (2002) and Brucker & Sarkar (2010). According to Meunier & Spedding (2004) and Redford, Lund & Coleman (2015), among others, the initialization should affect only the NW, and the wake should be universal further downstream.

2.3. Further details

Sponge layers are added at the vertical and transverse boundaries of the computational domain to absorb internal gravity waves, following the approach by Brucker & Sarkar (2010). The size of the sponge layers is $1$. Within the sponge layer, the velocity and temperature are damped using a Rayleigh damping function,

(2.14a,b)\begin{equation} -\!\phi(x_i) [U_{i}(x_{i},t)-U_{i,\infty}], \quad -\phi(x_i) \left[ T(x_{i},t)-T_{\infty}\right], \end{equation}

which are added to the right-hand side of the momentum equations and the temperature transport equation, respectively. Here, $U_{i,\infty }$ is the free stream velocity, and $T_{\infty }$ is the background temperature.

The grid used in the simulation is uniform in each of the three Cartesian directions, but the grid spacing varies in different directions. The resolution is chosen such that $\Delta x/\eta < 4$ and $\Delta y/\eta, \Delta z/\eta < 2$, where $\eta$ represents the Kolmogorov scale and is taken as the minimum value of $\eta (x_2, x_3, t)$. The computational domain has size $60D$ in the streamwise direction, and the transverse and vertical sizes are chosen such that the boundaries are sufficiently far from the wake, following the guidelines of Diamessis et al. (2011). Further details of the grid resolution and domain dimensions can be found in table 4.

Table 4. DNS details. The superscripts $1, 2, 3$ denote the different stages of a case, and $t$ is the physical time simulated. ‘No. ens.’ is the number of ensembles computed.

The DNS are run until they reach the viscous regime, which is characterized by $Re_{h}\,Fr_{h}^2<1$. This requires a long simulation time. To save computational resources, the domain is re-gridded during the simulation according to the method proposed by Gourlay et al. (2001), Redford, Castro & Coleman (2012) and Redford et al. (2015). The re-gridding is performed once or twice for each case. Specifically, the domain is re-gridded one time at $t = 110$ for the cases R20F02 and R20F10, two times at $t = 110$ and $t = 570$ for the case R20F50, and two times at $t = 68$ and $t = 265$ for the case R50F50. This leads to two or three stages for one case. This re-gridding process should have little to no influence on the wake according to Brucker & Sarkar (2010).

The DNS are performed using the in-house finite-volume solver AFiD. The code utilizes second-order-accurate spatial discretization with a staggered grid and a low-storage third-order Runge–Kutta method for temporal discretization. Further details of the code can be found in Van Der Poel et al. (2015), Ostilla-Mónico et al. (2015) and Kim & Moin (1985). The code has been used for stratified flows in Jain et al. (2022), Yang et al. (2015, 2020), Du, Zhang & Yang (2021), Bin et al. (2022) and Li & Yang (2021), among others. Here, a single R20F02 or R20F10 case requires approximately 2000 CPU hours; a single R20F50 realization requires 3000 CPU hours due to the increased simulation time; and a single R50F50 realization requires approximately 9000 CPU hours.

Finally, it is worth noting that the computational efficiency of computing multiple ensembles in a typically sized domain is superior to using an extremely long domain in the streamwise direction for our code. This is due to better parallelization and the convergence rate of the pressure Poisson equation. First, parallelization is generally more efficient for many small calculations than for a single large calculation for our code. Second, the convergence rate of the pressure Poisson equation depends on the mode with the longest wavelength. A long computational domain leads to a long wavelength mode, whose convergence rate is slow.

3. Results

3.1. Definitions

We define

(3.1)\begin{equation} {u_i}={\langle{u_i}\rangle} + u_i', \end{equation}

where $u_i$ is the instantaneous velocity in the $i$th Cartesian direction, $\left \langle {\cdot }\right \rangle$ denotes streamwise and ensemble averaging, and $u_i'$ is the fluctuation. The width and height of the wake are denoted as

(3.2a–e)\begin{equation} L_{H}, \quad L_{V}, \quad L_{Hk}, \quad L_{Vk} \quad {\rm and} \quad R_{i}, \end{equation}

where $L_H$ and $L_V$ are the distances from the centre of the wake to where the velocity deficit is half of its centreline value in the horizontal and vertical directions, $L_{Hk}$ and $L_{Vk}$ are the distances from the centre of the wake to where the TKE is half of its centreline value, and $R_{i}$ are given by

(3.3a,b)\begin{equation} R_i^2= A\frac{\displaystyle\int_C (x_i-x_i^c)^2{\langle{u_1}\rangle}^2 \,{\rm d}C}{{\displaystyle\int_C {\langle{u_1}\rangle}^2 \,{\rm d}C}}, \quad x_i^c = \frac{\displaystyle\int_C x_i{\langle{u_1}\rangle}^2 \,{\rm d}C}{{\displaystyle\int_C {\langle{u_1}\rangle}^2 \,{\rm d}C}}, \end{equation}

where $i=2, 3$, $x_i^C$ is the centre location of the wake, integration is in the entire vertical transverse plane, and $A=2$ is a constant (de Stadler et al. 2010).

3.2. Benefits of ensemble averaging

We compare the results from one and many realizations to illustrate the benefit of ensemble averaging. For illustration purposes, the discussion here will focus on case R20F02.

Figure 6 shows the contours of mean velocity deficit at $Nt=4$ (figures 6a,e), $Nt=30$ (figures 6b,f), $Nt=110$ (figures 6c,g) and $Nt=400$ (figures 6d,h) in the vertical transverse plane in case R20F02. The wake is in the NW regime at $Nt=4$, in the NEQ regime at $Nt=30$, and in the Q2D regime at $Nt=110$ and $400$. The results obtained from one hundred samples are shown in figures 6(a–d), and the results obtained from one realization are shown in figures 6(e–h). We see that the results in figures 6(a,e) are similar, and thus the results from one sample and one hundred samples are similar in the NW regime. The difference between the two becomes evident in the late wake. The results from one hundred samples exhibit symmetry with respect to $y=0$ and have the maximum velocity deficit located at the centre of the wake. In contrast, the results from one realization do not show such symmetry and have the maximum velocity deficit shifted away from the centre; see figures 6(f–h). This lack of symmetry and deviation from the expected behaviour suggest that results obtained from one realization are not statistically converged. A lack of statistical convergence partly explains the scattering of the centreline velocity deficit data in table 1.

Figure 6. Contours of the velocity deficit in the vertical transverse plane in case R20F02. (a–d) Results obtained from one hundred samples, and (e–h) results obtained from one sample, for (a,e) $Nt=4$, (b,f) $Nt=30$, (c,g) $Nt=110$, (d,h) $Nt=400$. For presentation purposes, we show only part of the domain. The white dashed lines indicate the core of the wakes, and they extend $2L_{Hk}$ and $2L_{Vk}$ in the $y$ and $z$ directions. The $+$ symbols indicate the wake centre defined in (3.3). The red circle symbols indicate the geometric centre of the wake, which is at $y=0$, $z=0$. Black lines indicate the contour lines of $1/4U_0$. The wake forms non-self-similar, diamond-shaped profiles in the late wake ($Nt=110, 400$).

Figure 7 shows the contours of the TKE in case R20F02 at $Nt=4$ (figures 7a,e,i), $Nt=30$ (figures 7b,f,j), $Nt=110$ (figures 7c,g,k) and $Nt=400$ (figures 7d,h,l). The results from one hundred samples are shown in figures 7(a–d), and the results from two independent realizations are shown in figures 7(e–h) and 7(i–l), respectively. The TKE is a second-order statistic, therefore obtaining statistically converged TKE is more challenging than obtaining statistically converged mean velocity deficit. While there was no apparent difference in the mean velocity deficit obtained from one and one hundred realizations in the NW regime, here we see a difference between one and one hundred ensembles as early as the NW regime. Insufficient averaging in the one-sample case leads to irregular variations across the vertical transverse plane, as seen in figures 7(e–h). In contrast, the results from one hundred samples, shown in figures 7(a–d), are regular and symmetric with respect to $y=0$. Specifically, the one-hundred-ensemble results display a single plateau at the centre of the wake, while the results from one realization show two-peak behaviour, particularly noticeable in figure 7(l). This discrepancy suggests that the two-peak behaviour reported in some realizations and in Brucker & Sarkar (2010) is a result of insufficient averaging.

Figure 7. Contours of the TKE in the vertical transverse plane at (a,e,i) $Nt=4$, (b,f,j) $Nt=30$, (c,g,k) $Nt=110$, and (d,h,l) $Nt=400$, in case R20F02. The dashed lines are contour lines. (a–d) One hundred ensembles. (e–h) Results from one realization. (i–l) Results from another realization.

Although it is not the focus of this study, higher-order statistics such as the terms in the TKE transport equation benefit more from ensemble averaging than low-order statistics such as the velocity deficit. Furthermore, the contours of the production term in the TKE transport equation are shown in figure 8. Comparing the results from one realization (figures 8e–h) with the results from one hundred realizations (figures 8a–d), it can be observed that ensemble averaging among one hundred realizations reveals structures of the production term that are not discernible from the result of a single realization. Specifically, the production term peaks towards the transverse boundary of the wake, and attains its minimum at the centre of the wake at $Nt=30$, 110 and 400, as shown in figures 8(b–d).

Figure 8. Same as figure 6 but for the production term in the TKE transport equation: (a,e) $Nt=4$, (b,f) $Nt=30$, (c,g) $Nt=110$, (d,h) $Nt=400$.

Figures 9(a–c) show the time evolution of the centreline velocity deficit, the wake width $R_2$, and the wake height $R_3$, in case R10F02. The results from one realization may be anywhere between the blue and red lines in the figures, incurring significant uncertainties for the centreline velocity deficit and the width of the wake. Specifically, the uncertainty in the centreline velocity deficit and the wake width is as large as 100 % in the late wake. The results in figure 9 also help to explain the scattering of the data in tables 1 and 3.

Figure 9. Evolution of (a) the centreline velocity deficit, (b) the wake width and (c) the wake height, in case R10F02. Here, SSB10 denotes the results in de Stadler et al. (2010), MAX and MIN are the maximum and minimum among the one hundred realizations, and AVG denotes results obtained from one hundred samples. Error bars in (a) and (b) characterize ${+}100\,\%$ error about the min value, and the error bar in (c) characterizes ${+}15\,\%$ error with respect to the min value; $R_2$ and $R_3$ are defined in (3.3).

The findings in this subsection prompt re-evaluation of how computational resources are utilized in DNS of stratified wake flows. Currently, the prevailing approach in the literature is to allocate all available resources to a single calculation at the highest practically attainable Reynolds number (Diamessis et al. 2011; de Bruyn Kops & Riley 2019; Zhou & Diamessis 2019). The rationale is that data at high Reynolds numbers will help to clarify the Reynolds number scalings of flow statistics. Nonetheless, we can gain meaningful insights only if the convergence error in the results is smaller than the change induced by varying the Reynolds number. This might not be the case when there is only one realization. Therefore, it might be worthwhile to allocate resources to increase the number of statistical samples. Doing so will likely lead to more reliable data for further analysis. In this section, we will revisit the scaling behaviours of statistics including the centreline velocity deficit, the height and width of the wake, and terms in the TKE budget equation, with the goal of extracting more accurate scaling estimates.

3.3. The convergence error and the number of realizations

The convergence error depends on the number of independent samples, which further depends on the domain size and the number of independent realizations. We have kept the domain size the same between different realizations, and here, we focus on the effect of the number of independent realizations (or statistical samples). The discussion here will focus on case R20F02, which, as we will see, is the most challenging case. The effects of the Froude number and the Reynolds number will be discussed towards the end of the subsection.

Transverse symmetry of the flow dictates that the transverse location of the wake centre $\langle x_2^c\rangle$ must be 0, therefore $|x_2^c|/L_{Hk}$ provides a measure of the convergence error in $\langle x_2^c\rangle$. Figures 10(a–d) present $|x_2^c|/L_{Hk}$ for all one hundred realizations for case R20F02 at $Nt=4$, $30$, $110$ and $400$, respectively. The dashed lines and dash-dotted lines in the figure indicate the 5 % and 50 % margins, representing the probabilities that data from one realization would incur errors above these lines. We see from figure 10 that there is a 5 % probability that data from one realization incur an error larger than 0.028, 0.053, 0.098 and 0.17 at $Nt=4$, 30, 110 and 400, respectively; and there is a 50 % probability that data from one realization incur an error larger than 0.0088, 0.019, 0.042 and 0.045 at these time instants, respectively. We see that data from one realization incurs large convergence errors. Furthermore, the convergence error increases as the wake evolves. The solid lines represent the ensemble average of one hundred samples. The error is between 0.002 and 0.003 for the four time instants investigated here. Although there has not been much discussion on an acceptable level of convergence error for DNS of temporally evolving wakes in stratified environments, some discussion is available on the acceptable level of error for DNS of channel flow (Oliver et al. 2014; Yang et al. 2021; Chen et al. 2023). There, a 1 % error in the mean velocity profile with respect to the bulk velocity is considered acceptable. Applying the same standard to the current study and limiting the discussion to the convergence error only, we can conclude that one hundred realizations lead to statistically converged $x_2^c$ statistics.

Figure 10. Plots of $|x_2^c|/L_{Hk}$ for all one hundred realizations for case R20F02 at (a) $Nt=4$, (b) $Nt=30$, (c) $Nt=110$, (d) $Nt=400$. The dashed line indicates the 5 % margin, the dash-dotted line indicates the 50 % margin, and the solid line indicates the ensemble average.

We discuss the number of realizations needed to obtain converged $\langle x_2^c\rangle$. First and foremost, any finite number of statistical samples gives only an estimate of the true $\langle x_2^c\rangle$. The estimator itself is a random variable with mean being the true value of $\langle x_2^c\rangle$ (0 in this case), and standard deviation inversely proportional to the square root of the number of samples, i.e. $\sigma \sim 1/\sqrt {n}$. If there is a $p\%$ probability that the error is above $e$ when there are $N$ samples, then one needs $N (e/e')^2$ samples to claim that there is a $p\%$ probability that the error is above $e'$. From the data presented in figure 10, we can estimate that there would be a 50 % probability that the error is above 0.01 if we take averages among 1, 4, 17 and 20 samples at $Nt=4$, 30, 110 and 400, respectively; and there would be a 5 % probability that the error is above 0.01 if we take averages among 8, 28, 95 and 289 samples at these four time instants. We see that more samples are needed to get statistically converged data in the late wake than in the early wake. This is likely because of the large-scale structures in the late wake that lead to a reduction in the number of effective statistical samples in the domain.

Before we proceed further, we examine the unstratified wake results shown in figure 11. These results pertain to case R20F$\infty$ at the time instants $t=8$ and $800$, which correspond to $Nt=4$ and $400$ in case R20F02. The results exhibit similarities to those in figure 10. With 40 samples, the errors in the ensemble average are 0.21 % and 0.4 % at times $t=8$ and 800. There are 5 % and 50 % probabilities that data from one realization incur errors above 0.05 and 0.021 at $t=800$, and above 0.024 and 0.0096 at $t=8$. Comparing these numbers to those for R20F02, we can conclude that it is more straightforward to get converged statistics for an unstratified wake (for a given $t$).

Figure 11. The centre location $x_2^c/L_{Hk}$ for all realizations for case R20F$\infty$ at (a) $t=8$ and (b) $t=800$. The dashed line indicates the 5 % margin, the dash-dotted line indicates the 50 % margin, and the solid line indicates the ensemble average.

In addition to $\langle x_2^c\rangle /L_{Hk}$, err defined below measures the spanwise asymmetry of the TKE, which also serves as a measure of the statistical convergence:

(3.4) \begin{equation} {\textit{err}}=\frac{\displaystyle\int_{{-}2L_{Vk}}^{2L_{VK}}\displaystyle\int_{{-}2L_{Hk}}^{0}\left\langle k\right\rangle {\rm d}y\,{\rm d}z - \displaystyle\int_{{-}2L_{Vk}}^{2L_{VK}}\displaystyle\int_{0}^{2L_{Hk}}\left\langle k\right\rangle {{\rm d}y}}{{\displaystyle\int_{{-}2L_{Vk}}^{2L_{VK}}}\displaystyle\int_{{-}2L_{Hk}}^{2L_Hk}\left\langle k\right\rangle{\rm d}y\,{\rm d}z}\times 100\,\%. \end{equation}

Figure 12 shows err as a function of the number of samples for case R20F02 at $Nt=4$, 30, 110 and 400. For any given number of samples, e.g. $n$, we randomly draw $n$ samples from all available realizations five times. As a result, for any given $n$, there are five data points. The trends are similar at the four time instants. The error decreases rapidly as the number of statistical samples increases from 1 to 20, and does not decrease significantly beyond 40–60 statistical samples. For the four time instants investigated here, the error is approximately 1 %–2 % when there are approximately 40–60 samples. Again, there is not a lot of discussion about the acceptable level of error in the context of DNS of stratified wake flows, but there is discussion in the context of DNS of channel flow. There, an approximately 2 %–3 % error in TKE is considered acceptable. We see that approximately 40–60 samples bring us to that regime, with the late wake requiring more realizations.

Figure 12. Plots of err in case R20F02 at (a) $Nt=4$, (b) $Nt=30$, (c) $Nt=110$, (d) $Nt=400$. For any given number of samples $n$, err is computed five times by randomly drawing $n$ samples from the available realizations. The solid line is at $err=2\,\%$.

Finally, figure 13 shows the results for cases R20F10, R20F50, R50F50 and R20F$\infty$ at time $t=800$ (corresponding to $Nt=400$ for case R20F02). We observe the following. First, the ensemble average in these cases reduces the convergence error to 2 % at $t=800$. Second, stratification has a notable effect on the result. Comparing figures 12(d) and 13(a–d), we see that for a given number of statistical samples, case R20F02 incurs a significantly larger convergence error than the other stratified cases (R20F10, R20F50 and R50F50), and the other stratified cases incur larger convergence errors than the unstratified case. Third, comparing figures 13(b,c), we see that there is no apparent Reynolds number effect from varying the Reynolds number by a factor 2.5.

Figure 13. Plots of err at $t=800$ in cases (a) R20F10, (b) R20F50, (c) R50F50, (d) R20F$\infty$. For any given number of samples $n$, err is computed five times by randomly drawing $n$ samples from the available realizations. The solid line is at $err=2\,\%$.

3.4. Centreline velocity deficit

Figure 14(a) presents the time evolution of the centreline velocity deficit $U_0$, the horizontal and vertical velocity fluctuations $u_h'$ and $w'$, and the square root of the TKE for cases R20F50 and R50F50. We see no apparent Reynolds number effect at the given range of the Reynolds numbers. In both cases, the centreline velocity decays as $Nt^{-2/3}$ in the NW regime up to $Nt\approx 1$; the decay rate decreases between $Nt=1$ and $Nt=50$ in the NEQ regime, and increases beyond $Nt=50$ in the Q2D regime. Contrary to the existing literature, we see no clear evidence of a $-0.25$ scaling or a $-0.76$ scaling. In fact, it is hard to find any power-law scaling. The decay rate changes continuously without staying constant for an extended period of time. This becomes more clear from figure 14(b), where we show the premultiplied centreline velocity deficit, i.e. $t^{0.25} U_0$. If there is a region within which the centreline velocity deficit follows $U_0\sim (Nt)^{-0.25}$, then we should observe a plateau in figure 14(b). However, no such a region can be identified. The $t^{0.76}$ compensated centreline velocity deficit is shown in the Appendix, and no plateau can be observed there either.

Figure 14. (a) Evolution of centreline velocity deficit $U_0$, horizontal and vertical velocity fluctuations $u_h'$ and $w'$, and root mean square of the TKE in cases R20F50 and R50F50. The solid lines represent the R50F50 case, and the dashed lines represent the R20F50 case. (b) Premultiplied centreline velocity deficit.

Figure 15(a) shows the results for case R20F10, and figure 15(b) shows the results for cases R10F02 and R20F02. Similar to previous cases, no significant Reynolds number effect is evident in these figures. Due to the strong buoyancy effect, a distinct NW regime is barely observable in the two F02 cases. If one defines the beginning of the NEQ regime to be the time instant at which the $U_0$ decay rate begins to decrease, then it is clear from figures 14(a) and 15(a,b) that the onset of the NEQ regime depends on the Froude number. Figures 15(c,d) show the premultiplied centreline velocity deficit. Once again, we do not see a plateau. We also examined the $U_0t^{2/3}$ compensated scalings in the Appendix, and no plateaus are identified there either.

Figure 15. (a,b) Evolution of centreline velocity deficit $U_0$, horizontal and vertical velocity fluctuations $u_h'$ and $w'$, and root mean square of the TKE in cases (a) R20F10, (b) R20F02 and R10F02. We use dashed lines for R20F02 and solid lines for R10F02. (c,d) Premultiplied centreline velocity deficit.

Finally, figure 16 shows the power-law exponent $b\equiv {\rm d}\log (U_0)/{\rm d}\log (t)$, or the decay rate of the centreline velocity deficit, as a function of the non-dimensional time $Nt$ for all four cases. The plots are shown in both log and linear scales. The estimates in Meunier et al. (2006) are included for comparison. We notice a discontinuity in the estimate at $Nt\approx 2$, prior to which is the NW regime and $b=-2/3$ (Tennekes & Lumley 1972; Meunier et al. 2006). We make the following observations. First, the initial period of the wake evolution is affected by the initialization. The gravitational ramping lasts until $t=8$, which corresponds to $Nt=4$, $0.8$ and $0.16$ for $Fr=2$, $10$ and $50$, respectively. In a plot whose $x$ axis is $Nt$, initialization will appear to have a more significant impact on wakes with lower Froude numbers. Second, the computed power-law exponents show a reasonable level of collapse. The value of the power-law exponent is approximately $-2/3$ for $Nt\lesssim 1$, decreases to approximately $-0.1$ between $Nt\approx 5$ and $Nt\approx 10$, then gradually increases. Third, the estimate in Meunier et al. (2006) provides a reasonable approximation of the data. According to the estimate, the decay rate attains its minimum 0.1 at $Nt=2$. Although the data suggest that the decay rate attains its minimum between $Nt\approx 5$ and $Nt\approx 10$, the value is indeed approximately 0.1. The estimate also captures the gradual increase of the decay rate in the stratified turbulence regime – although there are discrepancies, particularly when the flow is strongly stratified. Finally, the anticipated Q2D regime asymptotics, where $b$ approaches $-0.76$, are not observed in our data.

Figure 16. Power-law exponent $b={\rm d}\log (U_0)/{\rm d}\log (t)$ as a function of $Nt$. Solid lines represent the theoretical predictions in Meunier et al. (2006) and the symbols are our DNS results: (a) linear-log scale; (b) linear-linear scale.

3.5. Velocity fluctuations

The evolution of the horizontal and vertical velocity fluctuations $u_h'$, $w'$ and the square root of the TKE $K^{1/2}$ in R20F50 and R50F50 is shown in figure 14(a). From the plot, we observe that the results are not affected significantly by the Reynolds number, at least within the range of Reynolds numbers investigated in this study. Both the vertical and horizontal velocity fluctuations follow approximately a power-law decay with an exponent close to $-1.08$ in the NW regime between $Nt=0.8$ and $Nt=4$. This is consistent with the ${-1}$ scaling reported in Chongsiripinyo & Sarkar (2020) and is in agreement with the $t^{-1.08}$ scaling reported in Brucker & Sarkar (2010). The vertical velocity fluctuation continues to decay as $Nt^{-1.08}$ in the NEQ and Q2D regimes until the end of the simulation. On the other hand, the decay rate of the horizontal velocity fluctuations decreases as the flow enters the NEQ regime, and remains almost constant between $Nt=20$ and $Nt=80$. Due to this disparate decay rate between $u_h'$ and $w'$ in the NEQ and Q2D regimes, $u_h'$ becomes increasingly larger than $w'$ as the wake evolves. Consequently, the TKE is $K\approx u_h'/\sqrt {2}$.

Figure 15(a) shows the time evolution of $u_h'$, $w'$ and $K'$ in case R20F10, and figure 15(b) shows the corresponding results for cases R10F02 and R20F02. In the two F02 cases, a clear NW regime can barely be observed as the flow quickly enters the NEQ regime after the initialization period. In the NEQ and Q2D regimes, the decay rates of the horizontal velocity fluctuations and the square root of the TKE continue to increase. It is worth noting that the undulations in $w'$ observed in the two F02 cases are the result of rapid exchange between TKE and available potential energy. Interestingly, these undulations are limited to $w'$ in our DNS, while others have reported undulations in $U_0$ as well (Pal et al. 2017; Chongsiripinyo & Sarkar 2020). The vertical velocity fluctuation $w'$ follows a power law at these cases shortly after the calculation begins and all the way into the viscous regime. A least squares fit of the data gives power-law exponents $-1.19$ and $-1.29$ at the F10 and F02 cases, which are slightly different from the exponent in the F50 case.

Furthermore, by comparing the decay of $w'$ at different Froude numbers, it becomes evident that the exponent increases as $Fr$ decreases. An empirical fit of decay of $w'$ as a function of $Fr$ gives $w'=x^{-1.33+0.065\ln (Fr)}$. Nevertheless, further study is needed to clarify the exact effects of the Froude number and consolidate this empirical scaling estimate.

3.6. Length scale statistics

Figure 17 shows the time evolution of the wake width and height, denoted as $L_H$, $L_V$, $L_{Hk}$ and $L_{Vk}$. As before, $L_H$ and $L_V$ represent the distance from the wake centre to where the velocity deficit is half of its centreline value, while $L_{Hk}$ and $L_{Vk}$ measure the distance from the wake centre to where the TKE is half of its centreline value. First, we observe that the Reynolds number has a negligible effect on the data. From figure 17(a), we can see that the wake width $L_H$ grows as $(Nt)^{1/3}$ in the NW regime; its growth rate decreases at $Nt\approx 5$ and increases shortly after that. The $1/3$ power law provides a good working approximation of $L_H$ throughout the NW, NEQ and Q2D regimes. From figure 17(b), we see that the wake height $L_V$, also grows as $(Nt)^{1/3}$ in the NW regime; it then decreases in the NEQ regime, and reaches its minimum at $Nt\approx 50$ before it begins to increase again. The decrease of $L_V$ in the NEQ regime was also found in Davidson (2010) and Chongsiripinyo & Sarkar (2020), where it was attributed to buoyancy effects. The estimates in Meunier et al. (2006) are included in figures 17(a,b) for comparison. Despite some quantitative differences, the estimates from Meunier et al. (2006) capture the overall trend quite well. For instance, the estimate in Meunier et al. (2006) captures accurately the decrease and subsequent increase in the $L_H$ decay rate in the NEQ regime. Additionally, it also captures the suppression of the vertical length scale in the stratified flow regime. Finally, figures 17(c,d) depict $L_{Hk}$ and $L_{Vk}$. The two length scales behave similarly to $L_H$ and $L_V$, with $L_{Hk}$ growing faster than $L_H$, and $L_{Vk}$ decreasing more rapidly in the NEQ regime following a $-0.4$ power law (Chongsiripinyo & Sarkar 2020).

Figure 17. Wake width and height as functions of the non-dimensional time $Nt$: (a) $L_H$, (b) $L_V$, (c) $L_{Hk}$, (d) $L_{Vk}$. The solid lines represent estimates in Meunier et al. (2006). Again, normalization is by the diameter of the wake-generating body.

Before we conclude this subsection, we comment on the interplay between $U_0$, $L_H$ and $L_V$. If the wake profile is such that $U=U_0(x)\,f(\kern0.7pt y/L_H,z/L_V)$, where $f$ is generic function, then the momentum conservation would imply that $U_0L_HL_V={\rm const.}$ (Meunier et al. 2006). Figure 18(a) illustrates the evolution of $U_0L_VL_H$, which shows that this statistic is approximately constant after an initial transient period. We denote this constant as $C_m$. In figures 18(b,c), we compare $U_0$ with $C_m/(L_HL_V)$ and $C_m/[L_VC_H(Nt)^{1/3}]$, where $L_H\approx C_H(Nt)^{1/3}$. Here, $C_H$ is obtained from fitting the wake width data. We observe that the three quantities agree well. With the decay of $L_H$ contributing to a $-1/3$ power-law decay throughout, the results in figures 18(b,c) suggest that the decrease in the $U_0$ decay rate in the NEQ regime is mainly due to the decreasing $L_V$, while the increase in the $U_0$ decay rate in the Q2D regime is a result of the increasing $L_V$.

Figure 18. Plots of (a) $U_0L_vL_H$, (b) $C_m/(L_V L_H)$ and (c) $C_m/[L_VC_H(Nt)^{-1/3}]$, as functions of time. The plots (b,c) are shifted vertically by an arbitrary distance so that they collapse at $Nt=10$. Solid lines represent the estimated $U_0$, whereas the dashed line denotes $U_0$ in the DNS.

3.7. TKE budget

We devote this section to the TKE budget. The TKE transport equation reads

(3.5)\begin{equation} \frac{\partial k}{\partial t}={-}\mathcal{C}+\mathcal{P}-\mathcal{D}+\mathcal{B}-\mathcal{T}. \end{equation}

Here, $\mathcal{C}$ is the convection term

(3.6)\begin{equation} \mathcal{C} \equiv \langle{u_i}\rangle\,\frac{\partial{k}}{x_i}, \end{equation}

$\mathcal{P}$ is the production term

(3.7)\begin{equation} \mathcal{P} \equiv{-}\langle{u_i'u_j'}\rangle\frac{\partial{u_i}}{\partial{x_j}}, \end{equation}

$\mathcal{D}$ is the dissipation term

(3.8)\begin{equation} \mathcal{D} \equiv \frac{2}{Re_0}\,\langle{s_{ij}'s_{ij}'}\rangle, \end{equation}

where

(3.9)\begin{equation} s_{ij}'=\frac{1}{2}\left(\frac{\partial{u_i}'}{\partial{x_j}}+\frac{\partial{u_j}'}{\partial{x_i}}\right), \end{equation}

$\mathcal{B}$ is the buoyancy flux term

(3.10)\begin{equation} \mathcal{B} \equiv{-}\frac{1}{Fr^2}\,\langle{u_3'T'}\rangle, \end{equation}

and $\mathcal{T}$ is the transport/diffusion term

(3.11)\begin{equation} \mathcal{T} \equiv \frac{\partial{\langle{u_i'p'}}\rangle}{\partial{x_i}}+\frac{1}{2}\,\frac{\partial{\langle{u_i'u_j'u_j'}}\rangle}{\partial{x_i}}-\frac{2}{Re_0}\,\frac{\partial\langle u_j's_{ij}'\rangle}{\partial{x_i}}. \end{equation}

Pal & Sarkar (2015), Redford et al. (2015) and Chongsiripinyo & Sarkar (2020), among others, have reported the time evolution of the transverse and vertical integrated (integrated in both directions) terms in the TKE equation. Following these authors, we first report the transverse and vertical integrated terms. The integration is over $-2L_{Hk}< y<2L_{Hk}$ and $-2L_{Vk}< y<2L_{Vk}$, representing the core of the wake. Also, since the flow is not affected significantly by the Reynolds number in the investigated range, the analysis here will focus on a high-Froude-number case R50F50 and a low-Froude-number case R20F02.

Figures 19(a–d) show the evolution of the budget terms for case R50F50; the values of the budget terms vary vastly from one time to another, so we have plotted separately $0.2< Nt<1$, $1< Nt<3$, $3< Nt<15$ and $15< Nt<110$, respectively. The production term and the dissipation terms are the dominant terms for $Nt\lesssim 0.5$. Dissipation is by far the most dominant term between $Nt\approx 1$ and $Nt\approx 10$. From $Nt\approx 15$ until the end of the simulation, the production, transport and dissipation terms are the most significant terms. We note that buoyancy is never significant. It is also worth noting that our results align closely with the budget reported in Brucker & Sarkar (2010), suggesting that even one statistical sample is sufficient for transverse and vertically integrated budget terms. Figure 20 shows the transverse and vertically integrated budget terms for case R20F02. Unlike in case R50F50, the buoyancy term is a dominant term from the early stage until $Nt\approx 30$, which is consistent with the findings in Chongsiripinyo & Sarkar (2020). Additionally, the buoyancy term exhibits rapid sign changes during the early stage, alternating between acting as a sink and a source, which manifests the exchange between the available potential energy and TKE (Redford et al. 2015; Chongsiripinyo & Sarkar 2020). This behaviour explains the observed undulations in $K^{1/2}$ and $w'$ as depicted in figure 15(b). In the late wake, the dissipation and production terms emerge as dominant, which is similar to R50F50. Finally, the transport term is comparable to the production term between $Nt\approx 15$ and $Nt\approx 35$, but becomes negligibly small beyond $Nt\approx 100$.

Figure 19. The volume integrated TKE budget terms in case R50F50, $X=\int _V \mathcal {X}\,{\rm d}V$, where $\mathcal {X}$ represents the terms in the budget equation (3.5). Note that we did not divide the volume, following Brucker & Sarkar (2010). Plots of the budget terms as a function of $Nt$ (a) from 0.2 to 1, (b) from 1 to 3, (c) from 3 to 15, (d) from 15 to 110. Here, P, C, B, T and D represent the production, convective, buoyancy, transport and dissipation terms.

Figure 20. Same as figure 19 but for R20F02. Plots of the budget terms as a function of $Nt$ (a) from 4 to 15, (b) from 15 to 50, (c) from 50 to 180, (d) from 180 to 400.

By integrating in both the transverse and vertical directions, one obtains statistics that are less susceptible to convergence errors, but information pertaining to the transverse and vertical distributions is lost. In this study, we have obtained many statistical samples. They allow us to integrate in one direction yet still get converged data.

Figure 21 shows the transverse integrated budget terms as well as the vertically integrated budget terms in R50F50 at three time instants, $Nt=2$, 30 and 110, corresponding to the NW, NEQ and Q2D regimes, respectively. The transverse integrated terms are shown in figures 21(a,c,e) as functions of the vertical coordinate $z$. The vertically integrated terms are shown in figures 21(b,d,f) as functions of the transverse coordinate $y$. We first examine the terms in the NW regime. The dissipation term is the dominant term, and it has similar distributions in the vertical and the transverse directions – in the sense that one could not tell if the quantity is vertically integrated or transverse integrated from the plots themselves. The other terms are small. Next, we examine the terms in the NEQ regime. The production, dissipation and transport terms are the dominant terms. The wake is anisotropic, with significant differences between the transverse and vertically integrated terms. The transverse integrated production and transport terms have one peak at the wake centre, whereas the vertically integrated production and transport terms have two peaks at $y/L_{Hk}\approx \pm 0.6$. Additionally, the transverse integrated dissipation term shows double peaks at $z/L_{Vk}\approx \pm 0.6$, while the vertically integrated dissipation term has a single peak at the wake centre. The destruction terms, i.e. the dissipation and transport terms, are largely balanced by the production term. This explains the reduced decay rate of $K^{1/2}$ in the NEQ regime. Finally, we examine the terms in the Q2D regime: they are not qualitatively different from their counterparts in the NEQ regime. The dominant terms remain production, dissipation and transport. However, the balance between these terms changes. The production term decreases faster than the destruction terms, leading to a faster decay of $K^{1/2}$ in the Q2D regime. An interesting observation is that the vertically integrated production term becomes negative at the wake centre, indicating that the mean flow extracts energy from turbulence in that region.

Figure 21. The terms in the TKE budget equation in case R50F50: (a,b) $Nt=2$, (c,d) $Nt=30$, and (e,f) $Nt=110$. (a,c,e) Integrated in the transverse direction, $X=\int _y \mathcal {X}\,{{\rm d} y}$. (b,d,f) Integrated in the vertical direction, $X=\int _z \mathcal {X}\,{\rm d}z$. Here, P, C, B, T and D denote the production, convective, buoyancy, transport and dissipation terms. Considering the symmetry with respect to the centreline, we show only results from the wake centre.

Figure 22 presents the transverse integrated budget terms and the vertically integrated terms in case R20F02. The buoyancy term is the dominant term in the NW regime. From figure 20, we know that the term flips sign in the early wake. While the term is positive at $Nt=4$, it is negative at $Nt=5.5$ (not shown). The production, dissipation and transport terms are the dominant terms in the lake wake, similar to R50F50. However, unlike R50F50, the destruction terms in R20F02 are noticeably larger than the generation terms. This leads to a fast decaying $K^{1/2}$, as is evident in figure 15(b). It is interesting that the transport terms take energy from the wake centre and put it near the edge of the wake in the late wake, as shown in figures 21(d,f) and 22(d,f).

Figure 22. The terms in the TKE budget equation in case R20F02: (a,b) $Nt=4$, (c,d) $Nt=30$, and (e,f) $Nt=110$. (a,c,e) Integrated in the transverse direction. (b,d,f) Integrated in the vertical direction. Legend is the same as in figure 21.

4. Concluding remarks

We report DNS results of temporally evolving wake at Reynolds numbers from 10 000 to 50 000, and Froude numbers from 2 to 50. Instead of computing statistics from one realization, we take ensemble averages among 80–100 realizations. To assess the statistical convergence of our data, two statistical objects are considered: the transverse location of the wake centre, denoted as $x_2^c$, and the transverse asymmetry in the TKE, denoted as err. Both should be 0 when there are sufficiently many statistical samples. The analysis shows high convergence error when data are taken from one realization. The convergence error quickly decreases as the number of realizations/samples increases. The benefit, however, diminishes beyond 20 and 60 realizations before and after entering Q2D regimes, respectively.

The data allow us to assess the previously established scaling laws and establish new scaling estimates. First, the velocity deficit, denoted as $U_0$, does not seem to follow any power-law scaling in the late wake. The decay rate, defined as $-{\rm d}\log (U_0)/{\rm d}\log (t)$, attains its minimum 0.1 at a time between $Nt=5$ and $Nt=10$, and increases continuously as the wake evolves. No power-law scalings are found in $L_H$ or $L_V$ either, but $L_H\sim (Nt)^{1/3}$ proves to be a good working approximation of the wake width. Furthermore, while the wake is not self-similar, $U_0L_HL_V$ is approximately a constant after the initial transient period. It follows that $U_0\sim (Nt)^{-1/3}/L_V$. Hence we may conclude that the decrease in the $U_0$ decay rate in the NEQ regime is a result of the decrease in the wake's vertical size, and that the increase in the $U_0$ decay rate in the Q2D regime is a consequence of the increase in the wake's vertical size. Finally, the power-law decay of the vertical velocity fluctuation is verified, but there is no power-law decay in the horizontal velocity fluctuation.

In addition to low-order statistics, we report the transverse integrated budget terms, the vertically integrated budget terms, and the budget terms integrated in both cross-sectional directions. The behaviours of these terms as functions of the vertical and transverse coordinates have not received much attention due to the high number of ensembles needed to get converged statistics. Our data suggest that there are two peaks in the vertically integrated production and transport terms at $y=\pm 0.6 L_{Hk}$, and one peak at the wake centre when the terms are transverse integrated. It is interesting that the transverse integrated production term is negative at the wake centre, suggesting energy transfer from turbulence to the mean flow.

Finally, we add a caveat. Although the results do not seem to be affected significantly by the Reynolds number, the Reynolds numbers investigated here (i.e. 10 000 to 50 000) are quite limited. The Reynolds number of a wake behind an underwater vehicle can be as high as $O(10^9)$ (Li, Yang & Kunz 2022). Work by Meunier et al. (2006) and Brucker & Sarkar (2010) suggests that the NEQ regime extends as the Reynolds number increases. A lack of power-law decay of the centreline velocity deficit in the NEQ regime can be a result of the limited Reynolds number, which needs to be investigated further.