1. Introduction
Wakes in a stratified environment are common in both the ocean and the atmosphere. They are observed behind aircraft, wind turbines and underwater vehicles. Parametrising a stratified wake involves considering the inertia force, buoyancy force and viscous force, leading to the emergence of two non-dimensional numbers: the Froude number and the Reynolds number. The bulk Froude number, denoted as
$Fr$
, characterises the ratio of inertia to buoyancy and is expressed as
$Fr = U_b / (ND)$
. Here,
$U_b$
is the bulk velocity,
$D$
is the size of the wake-generating body, and
$N$
is the Brunt–Väisälä frequency defined as
$N\equiv \sqrt {(-g/\rho _0)\, {\rm d}\bar {\rho }(z)/{\rm d}z}$
, where
$\rho _0$
is the reference density, and
${\rm d}\bar {\rho }(z)/{\rm d}z$
is the density gradient along the vertical direction. The bulk Reynolds number, denoted as
$Re$
, delineates the ratio of inertia to viscous force and is defined as
$Re = U_b D / \nu$
, where
$\nu$
is the kinematic viscosity. The bulk Froude number and bulk Reynolds number define the initial state of a stratified wake. Additionally, one can define a local Froude number
$Fr_v = U_0 / (N L_V)$
and a local Reynolds number
$Re_v = U_0 L_V / \nu$
, where
$U_0$
is the centreline velocity deficit, and
$L_V$
is the height of the wake, defined as the distance from the centre of the wake to where the velocity deficit is half of its centreline value in the vertical directions. The centreline velocity deficit
$U_0$
decays as the wake evolves, while the height of the wake
$L_V$
remains roughly constant in the late wake. Consequently, both the local Froude number and the local Reynolds number decrease as the wake evolves.
Figure 1. Visualisation of the instantaneous vertical and transverse vorticities on the horizontal and vertical planes (top and bottom rows), with red indicating positive vorticity, and blue indicating negative vorticity. The three time instances roughly cover the NW, NEQ and Q2D regimes, respectively. The wake is a temporally evolving one.
Substantial advancements in understanding the flow physics of stratified wakes have been achieved since early work contributed to and reviewed in Lin & Pao (Reference Lin and Pao1979). In contrast to their unstratified counterparts, stratified wakes undergo multi-stage evolution. Spedding (Reference Spedding1997) categorises this evolution into three distinct flow regimes based on variations in the decay rates of the centreline velocity deficit: the near-wake (NW) regime, the non-equilibrium (NEQ) regime, and the quasi-two-dimensional (Q2D) regime. Buoyancy effects are negligible in the NW regime, but become increasingly significant in the NEQ and Q2D regimes. The transition points between these regimes are still debated (Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020; Li, Yang & Kunz Reference Li, Yang and Kunz2024). For instance, Spedding (Reference Spedding1997) reported the end of the NW regime at approximately
$Nt=1.7\pm 0.3$
, with the Q2D regime commencing at approximately
$Nt=50\pm 15$
. Conversely, Brucker & Sarkar (Reference Brucker and Sarkar2010) identified the termination of the NW regime at
$Nt=5$
and the onset of the Q2D regime at
$Nt=100$
.
Much of the prior work has focused on the scaling behaviours of some important low-order flow statistics, such as the centreline velocity deficit
$U_0$
, and the width and height of the wake. In the following, we review some of the prior work. In the NW regime, a stratified wake remains unaffected by buoyancy and behaves similarly to an unstratified wake, leading to
$U_0 \sim x^{-2/3}$
,
$L_v \approx L_h \sim x^{1/3}$
(Tennekes & Lumley Reference Tennekes and Lumley1972; Townsend Reference Townsend1976). Here,
$x$
measures the downstream distance, and
$L_v$
and
$L_h$
are the vertical and horizontal dimensions of the wake. Additionally, an
$x^{-1}$
scaling in
$U_0$
was observed in some studies and attributed to the anisotropic dissipation of turbulent kinetic energy (TKE) (Nedić et al. Reference Nedić, Vassilicos and Ganapathisubramani2013; Dairay, Obligado & Vassilicos Reference Dairay, Obligado and Vassilicos2015; Pal et al. Reference Pal, Sarkar, Posa and Balaras2017; Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020; Ortiz-Tarin, Nidhan & Sarkar Reference Ortiz-Tarin, Nidhan and Sarkar2021). In the NEQ regime, Spedding (Reference Spedding1997) reported
$U_0 \sim x^{-0.25}$
, and in the Q2D regime,
$U_0 \sim x^{-0.75}$
. These scalings were corroborated by several numerical studies (Diamessis, Domaradzki & Hesthaven Reference Diamessis, Domaradzki and Hesthaven2005; Brucker & Sarkar Reference Brucker and Sarkar2010; Diamessis, Spedding & Domaradzki Reference Diamessis, Spedding and Domaradzki2011), with slight variations that could be attributed to a lack of statistical convergence (Li et al. Reference Li, Yang and Kunz2024).
In addition to the centreline velocity deficit, the evolution of Reynolds stresses and their budgets are also important statistics. Brucker & Sarkar (Reference Brucker and Sarkar2010), de Stadler, Sarkar & Brucker (Reference de Stadler, Sarkar and Brucker2010) and Redford, Lund & Coleman (Reference Redford, Lund and Coleman2015), among others, conducted direct numerical simulations (DNS) of temporally evolving stratified wake flows in a streamwise periodic domain, and reported the planar-integrated TKE budgets as a function of time. Brucker & Sarkar (Reference Brucker and Sarkar2010) observed reduced vertical shear production
$P = -\langle u_x' u_z' \rangle\, \partial _z \langle u_x \rangle$
in the NEQ regime, and argued that it is responsible for the slow decay of
$U_0$
. Redford et al. (Reference Redford, Lund and Coleman2015) noted that the area-integrated TKE dissipation decays at the same rate in the NW and Q2D regimes, confirming that the Q2D regime is not two-dimensional turbulence (Godeferd & Cambon Reference Godeferd and Cambon1994). Historically, such studies were conducted with single DNS realisations of a given temporally evolving wake configuration. However, statistics derived from single DNS realisations often lack sufficient statistical convergence (Li et al. Reference Li, Yang and Kunz2024), making it difficult to examine the spatial structure of Reynolds stresses and their budgets. Consequently, the spatial structures of the Reynolds stress and the associated budget terms in their transport equations remain largely unexplored. Compared to planar-integrated terms, the spatial distribution of a term provides valuable additional information. The planar integration of the divergence of a flux equals the flux at the domain boundary, which is typically small. Thus planar integration alone does not reveal whether the term is uniformly small across the domain, or if it fluctuates locally and integrates to a small value.
Although the Reynolds stress budget in stratified wakes remains largely underexplored, it has been investigated extensively in other canonical flows, such as channel flows. Mansour, Kim & Moin (Reference Mansour, Kim and Moin1988) performed Reynolds stress budget analysis for turbulent channel flow at friction Reynolds number
$ Re_{\tau }=180$
; Hoyas & Jiménez (Reference Hoyas and Jiménez2008) examined the influence of Reynolds number on the Reynolds stress budget in channel flow at a higher Reynolds number, up to
$ Re_\tau =2000$
; Yuan & Piomelli (Reference Yuan and Piomelli2014) studied the Reynolds stress budget for open-channel flow over sandgrain roughness; and Lee & Moser (Reference Lee and Moser2019) conducted spectral analyses of TKE budget terms, exploring interactions between large and small scales. Analysing the Reynolds stress budget provides valuable insights into the mechanisms driving the production, dissipation and redistribution of turbulent stresses. This understanding is also crucial for advancing turbulence modelling, which aims to model the Reynolds stress accurately.
Turbulence modelling, especially within the framework of Reynolds-averaged Navier–Stokes (RANS), has a rich history (Pope Reference Pope2000; Kalitzin et al. Reference Kalitzin, Medic, Iaccarino and Durbin2005; Durbin Reference Durbin2018). Eddy-viscosity-based RANS models are grounded in the Boussinesq hypothesis,
$ -\langle u'_iu'_j\rangle = 2\nu _tS_{ij}- {2}/{3}k\delta _{ij}$
, where
$ \nu _t$
is the eddy viscosity,
$ S_{ij}$
is the strain rate tensor,
$ k$
represents TKE, and
$ \delta _{ij}$
is the Kronecker delta. Early eddy-viscosity models were algebraic (Smith & Cebeci Reference Smith and Cebeci1967; Baldwin & Lomax Reference Baldwin and Lomax1978), limiting their applicability. Subsequent developments in the 1980s and 1990s led to the widely adopted transport models used today (Launder & Sharma Reference Launder and Sharma1974; Chien Reference Chien1982; Wilcox Reference Wilcox1988; Spalart & Allmaras Reference Spalart and Allmaras1992; Menter Reference Menter1994). Beyond eddy-viscosity models, Reynolds stress models solve the transport equations for Reynolds stress components directly (Launder, Reece & Rodi Reference Launder, Reece and Rodi1975; Speziale, Sarkar & Gatski Reference Speziale, Sarkar and Gatski1991). A detailed overview of the terms in these transport equations is provided in § 3. Despite their widespread use in engineering applications at high Reynolds numbers, RANS simulations often underperform in predicting stratified wakes (Wall & Paterson Reference Wall and Paterson2020; Jain et al. Reference Jain, Huang, Li, Yang and Kunz2023), largely due to the anisotropic nature of the late wake and the failure in capturing such anisotropy. One objective of this study is to identify the sources of error in existing turbulence models and identify potential avenues for improvements.
In this study, we present budget analyses of the Reynolds stresses for stratified wakes. The DNS are temporally evolving wakes. Ensemble averages are conducted to mitigate statistical errors that are otherwise present in data obtained from single DNS realisations. We report Reynolds stresses and their budgets in the transverse-vertical plane for the NW, NEQ and Q2D regimes. Additionally, we assess the Boussinesq hypothesis and the model in Speziale et al. (Reference Speziale, Sarkar and Gatski1991). The primary objective is to enhance our understanding of flow physics and inform future modelling efforts.
The remainder of this paper is organised as follows. We present the details of the DNS data in § 2. The Reynolds stress transport equations and various turbulence closures are shown in §§ 3 and 4, respectively. Results of the budget analysis are presented in § 5. Finally, we conclude in § 6.
2. Direct numerical simulation details
We consider a temporally evolving wake in a linearly stratified environment. A temporally evolving frame is related to a spatially evolving one through
$X = U_b t$
, where
$X$
measures the distance from the wake-generating body in the streamwise direction,
$U_b$
is the bulk velocity, and
$t$
is the elapsed time of the temporally evolving wake. The wake is statistically homogeneous in the
$x$
direction, and decays over time. The gravitational force acts in the
$-z$
direction. The domain extends
$60D$
in the streamwise direction (
$x$
). The transverse and vertical sizes are such that the boundaries are sufficiently far from the wake, following the practice in Diamessis et al. (Reference Diamessis, Spedding and Domaradzki2011). The grid resolution ensures
$\Delta x / \eta \lt 4$
and
$\Delta y / \eta , \Delta z / \eta \lt 2$
, where
$\eta$
is the Kolmogorov length scale, following Moin & Mahesh (Reference Moin and Mahesh1998) and Brucker & Sarkar (Reference Brucker and Sarkar2010). Sponge layers with thickness
$1D$
are added at the transverse and vertical boundaries of the computational domain to absorb internal gravity waves, following Brucker & Sarkar (Reference Brucker and Sarkar2010) and de Stadler et al. (Reference de Stadler, Sarkar and Brucker2010). In the sponge layer, the velocity and temperature/density are damped using Rayleigh damping functions,
\begin{equation} -\phi (x_i) \left [U_{i}(x_{i},t)-U_{i,\infty }\right ], \quad -\phi (x_i) \left [ T(x_{i},t)-T_{\infty }\right ], \quad i=2,3, \end{equation}
which are added to the right-hand sides of the momentum and temperature transport equations, respectively. Here,
$U_{i,\infty }$
is the free stream velocity deficit,
$T_{\infty }$
is the background temperature, and
$\phi (x_i)$
is the damping ratio, which is a constant.
Figure 2. A sketch of the computational set-up of a temporally evolving stratified wake flow.
The initial condition comprises the mean wake profile and turbulent fluctuations. The mean streamwise velocity profile is
\begin{equation} \left \lt u \right \gt (r) = U_{0} \exp \left ( -\frac {1}{2}\left ( \frac {r}{r_0} \right )^2 \right ), \end{equation}
where
$U_{0}$
is the initial centreline velocity deficit,
$r$
is the radial distance from the wake centre, and
$r_0$
is the initial wake radius. Here, we set
$U_0 = 0.11U_b$
and
$r_0 = {1}/{2}D$
, following Dommermuth et al. (Reference Dommermuth, Rottman, Innis and Novikov2002) and Brucker & Sarkar (Reference Brucker and Sarkar2010). Note that we solve for the velocity deficit and therefore
$u=0$
in the freestream. The turbulent fluctuations are initialised as follows. First, a random field is generated from a uniform distribution between 0 and 1. Second, the field is made divergence-free and digitally filtered such that it follows the energy spectrum
\begin{equation} E(k) \sim \left ( \frac {k}{k_0} \right )^4 \exp\left({-2\left ( \frac {k}{k_0} \right )^2}\right), \end{equation}
where
$k_0 = 4$
(de Stadler et al. Reference de Stadler, Sarkar and Brucker2010). Third, the fluctuations are damped exponentially away from the centre according to the damping function
\begin{equation} g(r) = I\left ( 1 + \frac {r}{r_0} \right )^2 \exp \left ( -\frac {1}{2} \left ( \frac {r}{r_0} \right )^2 \right ), \end{equation}
where
$I$
is the initial turbulence intensity, set to
$0.055$
. During the initialisation period, both the mean flow in (2.2) and the turbulence intensity in (2.4) are maintained. The turbulence spectrum in (2.3) is, however, not maintained, and is allowed to evolve until it equilibrates with the specified mean flow and turbulence intensity.
Further details of the initialisation and validation of the methodology can be found in Brucker & Sarkar (Reference Brucker and Sarkar2010) and Li et al. (Reference Li, Yang and Kunz2024), and are not detailed here for brevity. Applying the above initialisation 100 times at a given flow condition leads to 100 independent realisations, and we take ensemble averages among these realisations. The purpose of performing 100 independent realisations is to obtain sufficient data for statistical convergence. While extending the computational domain by a factor of 100 could, in theory, yield the same amount of data, this approach is far less computationally efficient. Parallelisation is more effective for multiple small-scale simulations, as it reduces inter-processor communication. Additionally, a longer computational domain results in long-wavelength modes, significantly slowing the convergence rate of the pressure Poisson equation. This trade-off is the primary reason why we opted for multiple independent realisations rather than a single large-scale calculation. The rationales behind our choice of the number of ensembles are detailed in Li et al. (Reference Li, Yang and Kunz2024), where a statistical convergence study is presented. Specifically, by computing the spanwise asymmetry in TKE as a function of the number of ensembles, it is concluded that for the domain sizes used in this study, approximately 60–80 ensembles are sufficient to reduce the asymmetry error to less than 2 % in the mean flow and the Reynolds stresses. Here, 80–100 ensembles are used.
Table 1 shows the details of our DNS, including the initial Reynolds and Froude numbers, the grid number, the domain size, the total integration time (in terms of
$t=t^*U_b/D$
, where
$t^*$
is the dimensional time), the buoyancy scaled time
$Nt$
, and the number of independent realisations. The nomenclature of the cases is R[Re]F[Fr], where Re is the initial Reynolds number divided by 1000, and Fr is the initial Froude number.
Table 1. The DNS details. The superscripts 1, 2, 3 denote the different stages of a case. For a given stage that runs from
$t_s$
to
$t_e$
, the time
$t$
listed in the table corresponds to
$t_e$
rather than
$t_e-t_s$
. That is,
$t$
is the non-dimensional simulation time at the end of a given stage. Here,
$Nt$
is the non-dimensional buoyancy scaled total runtime,
$L_x$
,
$L_y$
and
$L_z$
are non-dimensionalised using
$D$
, and ‘No. ens.’ is the number of independent realisations.
The size of the Kolmogorov length scale increases as the wake decays. Consequently, the initial grid becomes overly fine as the wake evolves. Following the practices in Brucker & Sarkar (Reference Brucker and Sarkar2010), de Stadler et al. (Reference de Stadler, Sarkar and Brucker2010) and Redford et al. (Reference Redford, Lund and Coleman2015), we periodically coarsen the DNS grid to improve computational efficiency. This leads to multiple stages for one case. Specifically, we coarsen the grid once at
$t = 110$
for case R20F02, twice at
$t = 110$
and
$t = 570$
for case R20F50, and twice at
$t = 68$
and
$t = 265$
for case R50F50, leading to two and three stages for F02 and F50 cases, respectively. During regridding, the domain size is also increased in the transverse and vertical directions to accommodate the wake growth in the corresponding directions. In the region outside the old domain, the velocity and temperature fields are set equal to
$\psi _{i}=\psi _{{old}}\exp[-(r_i - r_{{old}})^2]$
, where
$\psi _{{old}}$
is the value of the flow variables, velocity and temperature, at the boundary of the old domain excluding sponge layers, and
$\psi _i$
is the value of the flow variables in the new domain. This allows the flow field to decay smoothly to zero in the new domain. Note that because we solve for the velocity deficit and the temperature fluctuations, their values are 0 in the freestream.
An in-house adaptation of the open-source code AFiD (Van Der Poel et al. Reference Van Der Poel, Ostilla-Mónico, Donners and Verzicco2015) is employed for the DNS. Non-dimensionalisation is achieved using the bulk velocity
$U_b$
, the freestream temperature
$T_0$
, the freestream density
$\rho _0$
, and the diameter of the wake generator
$D$
. The non-dimensional flow variables are
\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}
where the superscript
$^*$
denotes the dimensional variables. All variables reported hereafter are non-dimensional unless otherwise noted. The code solves the following non-dimensional Navier–Stokes equation with the Boussinesq approximation:
\begin{equation} \frac {\partial {u_{i}}}{\partial {t}} + \frac {\partial {\left ( u_{j} u_{i} \right )}}{\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{equation}
Here,
$\rho '$
is the instantaneous density fluctuation defined in the decomposition
\begin{equation} \rho =\rho _0+\bar {\rho }(x_3)+\rho '(x_i,t), \end{equation}
where the reference density is
$\rho _0=1$
according to the present non-dimensionalisation, and
$\bar {\rho }(x_3)$
is the average of the deviation of the density from the reference density at a given vertical location
$x_3$
, not a function of time. In our DNS, it represents a linear density stratification in the background flow field. A similar decomposition can be defined for the temperature
\begin{equation} T=T_0 + \bar {T}(x_3)+T'(x_i,t), \end{equation}
where the reference temperature is
$T_0=1$
per the present non-dimensionalisation,
$\bar {T}(x_3)$
is the average of the temperature deviation from the reference temperature at a given vertical location
$x_3$
, and
$T'$
is the temperature fluctuation. The density fluctuation is related to the temperature fluctuation according to
\begin{equation} \bar {\rho }(x_3)=-\beta\, \bar {T}(x_3), \quad\rho '= -\beta T',\quad\beta = -\left ( \frac {\partial {\rho }}{\partial {T}} \right )_P, \end{equation}
where
$\beta$
is the non-dimensional thermal expansion coefficient and equals 1. The temperature fluctuation itself is solved according to
\begin{equation} \frac {\partial {T'}}{\partial {t}} + u_{i} \frac {\partial ({\bar {T}(x_3)+T'})}{\partial {x_{i}}} = \frac {1}{Re\,Pr} \frac {\partial ^2({\bar {T}(x_3)+T'})}{\partial {x_{i}} \partial {x_{i}}}. \end{equation}
Here,
$Pr = \nu / \alpha$
is the Prandtl number, which is set to unity in the current simulations, and
$\alpha$
is the thermal diffusivity. The code employs second-order central discretisation with a staggered grid and a low-storage third-order Runge–Kutta method for time stepping. Further details of the code can be found in Ostilla-Mónico et al. (Reference Ostilla-Mónico, Yang, Van Der Poel, Lohse and Verzicco2015) and Kim & Moin (Reference Kim and Moin1985), and the references cited therein. The code has been used previously for stratified flows; see e.g. Jain et al. (Reference Jain, Pham, Huang, Sarkar, Yang and Kunz2022b
),Yang et al. (Reference Yang, Van Der Poel, Ostilla-Mónico, Sun, Verzicco, Grossmann and Lohse2015, Reference Yang, Chen, Verzicco and Lohse2020), Du, Zhang & Yang (Reference Du, Zhang and Yang2021), Bin et al. (Reference Bin, Yang, Yang, Ni and Shi2022) and Li & Yang (Reference Li and Yang2021), among others.
3. The RANS and Reynolds stress transport equations
This study focuses on Reynolds stress budget analysis. For completeness, we present Reynolds stress transport equations in this section. We define the Reynolds decomposition
\begin{equation} {u_i}={\left \lt {u_i}\right \gt } + u_i', \end{equation}
where
$u_i$
is the instantaneous velocity in the
$i$
th Cartesian direction,
$\langle \cdot \rangle$
denotes streamwise and ensemble averaging, and
$u_i'$
is the velocity fluctuation. The Reynolds stress tensor
$R_{ij}$
is
$\langle u_i'u_j'\rangle$
, which appears in the following Reynolds averaged NS equation. The streamwise-and-ensemble-averaged streamwise momentum equation reads
\begin{equation} \frac {\partial \left \lt u \right \gt }{\partial t} + \left \lt v \right \gt \frac {\partial \left \lt u \right \gt }{\partial y} + \left \lt w \right \gt \frac {\partial \left \lt u \right \gt }{\partial z} = \nu \left ( \frac {\partial ^2 \left \lt u \right \gt }{\partial y^2} + \frac {\partial ^2 \left \lt u \right \gt }{\partial z^2} \right ) - \left ( \frac {\partial \left \lt u'v' \right \gt }{\partial y} + \frac {\partial \left \lt u'w' \right \gt }{\partial z} \right ). \end{equation}
Note that since the streamwise direction is homogeneous, there is no mean gradient in the streamwise direction, i.e.
$\partial /\partial x=0$
.
The transport equation of the Reynolds stress
$R_{ij}=\langle u'_iu'_j \rangle$
reads
\begin{equation} \frac {\partial {R_{ij}}}{\partial {t}}+U_{k}\frac {\partial {R_{ij}}}{\partial {x_k}} = P_{ij} +(D_{T,ij}+D_{v,ij}+D_{p,ij})+B_{ij}+\phi _{ij}-\epsilon _{ij}, \end{equation}
where
$P_{ij}$
,
$D_{T,ij}$
,
$D_{v,ij}$
,
$D_{p,ij}$
,
$B_{ij}$
,
$\phi _{ij}$
and
$\epsilon _{ij}$
are the production, turbulent diffusion, viscous diffusion, pressure diffusion, buoyancy flux, pressure–strain correlation and turbulence dissipation. Specifically, the production term is
\begin{equation} P_{ij}=-R_{ik}\frac {\partial {U_j}}{\partial {x_k}}-R_{jk}\frac {\partial {U_i}}{\partial {x_k}}, \end{equation}
the diffusion terms are
\begin{equation} D_{p,ij}=-\frac {\partial }{\partial {x_k}}\left ( \big \lt u_i'p'\big \gt \delta _{jk}+\big \lt u_j'p'\big \gt \delta _{ik} \right ), \quad D_{v,ij}=\nu \frac {\partial ^2 {R_{ij}}}{\partial {x^2_k}}, \quad D_{T,ij}=-\frac {\partial {R_{ijk}}}{\partial {x_k}}, \end{equation}
the buoyancy term is
\begin{equation} B_{ij}= \frac {1}{Fr^2}\left ( \delta _{i3}\big \lt u_j'T'\big \gt +\delta _{j3}\big \lt u_i'T' \big \gt \right ), \end{equation}
the pressure–strain correlation term is
\begin{equation} \phi _{ij}=\left \lt {p'\left (\frac {\partial {u'_i}}{\partial {x_j}}+\frac {\partial {u'_j}}{\partial {x_i}} \right )}\right \gt , \end{equation}
and the turbulent dissipation term is
\begin{equation} \epsilon _{ij}=2\nu \left \lt {\frac {\partial {u'_i}}{\partial {x_k}}\frac {\partial {u'_j}}{\partial {x_k}}} \right \gt . \end{equation}
The terms on the right-hand side of the Reynolds stress transport equations are categorised as source/sink terms, spatial redistribution terms, or inter-component redistribution terms. The production, dissipation and buoyancy terms are source/sink terms. Their spatial integration is not necessarily zero, therefore these terms lead to net generation or net destruction of the Reynolds stresses. The viscous diffusion, pressure diffusion and turbulent diffusion terms are spatial redistribution terms that, according to Green’s theorem, do not contribute to the net generation or destruction of Reynolds stresses. The pressure–strain correlation term is an inter-component redistribution term for the three normal stresses, as it transfers energy from one normal component to another with
$\sum _{k=1}^3\phi _{kk}=0$
. For Reynolds shear stresses, however, the pressure–strain term serves as a source/sink, leading to net generation or destruction of the stresses.
4. Closure models
Comparing data to the existing closures allows us to examine the current understanding of the stratified wake flow physics. We consider two types of models: eddy-viscosity type models and Reynolds stress type models. Eddy-viscosity type models solve eddy viscosity and approximate Reynolds stress with the Boussinesq hypothesis. On the other hand, Reynolds stress type models solve the Reynolds stress transport equations, and the diffusion, dissipation and pressure–strain correlation terms must be closed. A commonly used diffusion model is the Daly–Harlow model (Daly & Harlow Reference Daly and Harlow1970). This model makes use of the Reynolds stress tensor to define an anisotropic diffusion coefficient, which is subsequently used to close the turbulent diffusion term:
\begin{equation} D_{T,ij} = \frac {\partial }{\partial x_k} \left [ \left ( \nu \delta _{kl} + C_s \frac {k \left \lt u_k' u_l' \right \gt }{\epsilon } \right ) \frac {\partial \big \lt u'_i u'_j \big \gt }{\partial x_l} \right ], \end{equation}
where
$C_s$
is a constant (Launder Reference Launder and Lumley1990). A commonly used model for the dissipation tensor is the isotropic dissipation model
\begin{equation} \epsilon _{ij} = \tfrac {2}{3} \epsilon \delta _{ij}, \end{equation}
where
$\epsilon$
is the modelled dissipation rate (Launder et al. Reference Launder, Reece and Rodi1975) and is solved separately. This model assumes isotropy, which, as we will see, fails in stratified wake flows. Anisotropic dissipation models have been explored recently (Wall & Paterson Reference Wall and Paterson2020; Jain et al. Reference Jain, Huang, Li, Yang and Kunz2023), but are not examined here. Finally, the pressure–strain correlation term has received much attention. Speziale et al. (Reference Speziale, Sarkar and Gatski1991) modelled the pressure–strain correlation as a function of the strain-rate tensor. Although not intended to model flow with buoyancy effects, the Speziale, Sarkar and Gatski (SSG) model remains one of the most widely applicable models for the pressure–strain correlation term. Conducting a Taylor expansion and keeping the first-order terms only, they obtained
\begin{align} \phi _{ij} =& - C_1 \epsilon + C_1^* \mathcal {P}_{kk} b_{ij} + C_2 \epsilon \left ( b_{ik} b_{kj} - \tfrac {1}{3} b_{mn} b_{mn} \delta _{ij} \right ) + \left ( C_3 - C_3^* \sqrt {b_{kl} b_{kl}} \right ) k S_{ij} \nonumber\\ &{}+ C_4 k \left ( b_{ik} S_{jk} + b_{jk} S_{ik} - \tfrac {2}{3} b_{mn} S_{mn} \delta _{ij} \right ) + C_5 k \left ( b_{ik} W_{jk} + b_{jk} W_{ik} \right ), \end{align}
where
$b_{ij}$
is the Reynolds stress anisotropy tensor defined as
$b_{ij}\equiv \langle u'_iu'_j \rangle /2k-1/3\delta _{ij}$
, and
$W_{ij}$
is the rotation tensor
$W_{ij} \equiv 1/2(\partial u_i/\partial x_j - \partial u_j/ \partial x_i)$
. The constants are
\begin{align} C_1 = 3.4, \!\quad\! C_1^* = 1.8, \!\quad\! C_2 = 4.2, \!\quad\! C_3 = \tfrac {4}{5}, \!\quad\! C_3^* = 1.3, \!\quad\! C_4 = 1.25, \!\quad\! C_5 = 0.4, \end{align}
whose calibration involves return to isotropy (Sarkar & Speziale Reference Sarkar and Speziale1990), rapid distortion theory (Crow Reference Crow1968), homogeneous shear flow (Tavoularis & Corrsin Reference Tavoularis and Corrsin1981) and rotating shear flow (Bertoglio Reference Bertoglio1982).
5. Results
We present ensemble-averaged Reynolds stresses and the terms in their transport equations, with a focus on the core of the wake:
$ y/D \leqslant 2L_{Hk}$
and
$ z/D \leqslant 2L_{Vk}$
, where
$ L_{Hk}$
and
$ L_{Vk}$
represent the wake’s width and height, defined as the distances from the geometric centre to the locations where the TKE is half of its centreline value. In principle, we can examine all budget terms for the six Reynolds stress components, resulting in
$ 7 \times 6 = 42$
terms (although some are zero or negligibly small). Each budget term, such as production, is a function of
$ t$
,
$ y$
,
$ z$
, as well as the Reynolds and Froude numbers. A comprehensive presentation of all results offers insights into the spatial organisation of the budget terms, their temporal evolution, inter-term balances, and their dependence on Reynolds and Froude numbers. However, presenting all results would lead to an overwhelming number of figures and plots. To maintain a concise narrative, this study focuses on one representative case, R50F50, while referencing results for other flow conditions as appropriate. We illustrate the spatial structure of the budget terms and their temporal evolution, using separate plots for different Reynolds stress components. This allows readers to compare the behaviour of the same term, such as production, across different budget equations by referencing various figures in this section. Additionally, readers can assess the dependence of these budget terms on Reynolds and Froude numbers by comparing results in this section with those in Appendix A. This approach is justified by the observation that most budget terms exhibit limited sensitivity to variations in Reynolds and Froude numbers – at least within the Reynolds number range investigated in this study. Our presentation also aligns with prior studies of Reynolds stress budgets in channel flows, where investigations initially focus on a single flow condition (Mansour et al. Reference Mansour, Kim and Moin1988) before exploring Reynolds number dependence (Hoyas & Jiménez Reference Hoyas and Jiménez2008; Lee & Moser Reference Lee and Moser2019).
5.1. Reynolds stresses and the Boussinesq hypothesis
Figure 3 presents the Reynolds stress profiles of case R50F50 at
$Nt=1$
, 4, 30, 110, roughly covering the NW, NEQ and Q2D regimes. Although the exact range of each flow regime is still an open topic, the structures of the flow are expected to be similar within a given flow regime. By presenting results at only four time instances, the discussion here and in subsequent subsections emphasises the spatial structure rather than the temporal evolution. Topics such as the decay rate of
$ \langle w'w' \rangle$
, which have been addressed extensively in previous studies (Zhou & Diamessis Reference Zhou and Diamessis2019; Brucker & Sarkar Reference Brucker and Sarkar2010; Li et al. Reference Li, Yang and Kunz2024), are not covered. We apply the same colour map range for data at each time instant, but adjust the colour map between different time instants to keep information regarding wake decay. The top three rows of figure 3 show the evolution of the three normal stresses. Anisotropy develops as the wake evolves, as indicated by the increasing disparity between the three normal stresses. At
$Nt = 110$
, the last reported time in figure 3, the vertical component of the normal stress,
$\langle w'w' \rangle$
, is negligible compared to the longitudinal and transverse components,
$\langle u'u' \rangle$
and
$\langle v'v' \rangle$
, due to the suppression of vertical motion by stratification in the late wake. Between the longitudinal and transverse components,
$\langle u'u' \rangle$
dominates
$\langle v'v' \rangle$
at
$Nt = 4$
, whereas
$\langle v'v'\rangle$
dominates
$\langle u'u'\rangle$
at
$Nt = 110$
. In addition to their magnitudes, the three normal stresses also differ in their structures. In the late wake, we observe two peaks in the longitudinal component
$\langle u'u' \rangle$
, and only one peak in the transverse component
$\langle v'v' \rangle$
. According to figure 3, the two peaks in
$\langle u'u' \rangle$
already appear at
$Nt = 30$
. Prior to
$Nt=30$
, a local minimum in
$\langle u'u' \rangle$
emerges at the geometric centre of the wake at
$Nt = 4$
. We will revisit this peculiar behaviour of
$\langle u'u' \rangle$
in § 5.2, and will demonstrate that the two peaks in
$\langle u'u' \rangle$
are due to the production term.
Figure 3. Contours of all six Reynolds stresses in R50F50 at
$Nt=1$
, 4, 30 and 110. The plot does not reflect the true aspect ratio of the wake, which is indicated by the rectangles in the first row. Different rows are for different Reynolds stresses, and different columns for different time instants. Top to bottom:
$ \langle u'u' \rangle$
,
$ \langle v'v' \rangle$
,
$ \langle w'w' \rangle$
,
$ \langle u'v' \rangle$
,
$ \langle u'w' \rangle$
and
$ \langle v'w' \rangle$
. Left to right:
$Nt=1$
, 4, 30 and 110.
The bottom three rows of figure 3 illustrate the evolution of the three Reynolds shear stresses. Unlike the normal stresses that are positive definite, Reynolds shear stress can be positive or negative. From figure 3, we observe that the
$\langle u'v' \rangle$
stress is positive for
$y \gt 0$
and negative for
$y \lt 0$
at all four time instants. In other words, the
$-\langle u'v' \rangle$
term, which appears on the right-hand side of the Reynolds-averaged streamwise momentum equation, is negative for
$y \gt 0$
and positive for
$y \lt 0$
. This implies positive correlations between positive (negative)
$u'$
and positive (negative)
$v'$
when
$y \gt 0$
, and positive (negative)
$u'$
and negative (positive)
$v'$
when
$y \lt 0$
, suggesting a net momentum deficit flux from the geometric centre of the wake to the
$+y$
and
$-y$
directions, which subsequently leads to wake expansion in the transverse direction. Note that we solve for the velocity deficit, therefore a positive
$U_0$
or
$u$
indicates a velocity deficit. We may conduct a similar analysis for the other two Reynolds shear stress components. For instance, the
$-\langle u'w' \rangle$
component is negative for
$z \gt 0$
and positive for
$z \lt 0$
at
$Nt = 1$
. The observed correlation between positive (negative)
$u'$
and positive (negative)
$w'$
for
$z \gt 0$
, and positive (negative)
$u'$
and negative (positive)
$w'$
for
$z \lt 0$
, leads to a net momentum deficit from the geometric centre of the wake to the
$+z$
and
$-z$
directions in the
$z \gt 0$
region and the
$z \lt 0$
region, respectively, leading to wake expansion in the vertical direction at
$Nt = 1$
in the early wake.
These correlations lead to positive correlations between
$v'$
and
$w'$
in the first and third quadrants, and negative correlations between them in the second and fourth quadrants. However, as buoyancy becomes increasingly important in the late wake, vertical turbulent motions are suppressed. Consequently, the two Reynolds shear stress components that involve
$w'$
, i.e.
$\langle u'w' \rangle$
and
$\langle v'w' \rangle$
, become negligibly small compared to the
$\langle u'v' \rangle$
component in the late wake. Despite their small magnitudes, both
$\langle u'w' \rangle$
and
$\langle v'w' \rangle$
have complex internal structures. For instance, we see horizontal stripes with alternating signs in the
$\langle u'w' \rangle$
component. This is more evident in figure 4, where we show contours of
$\langle u'w' \rangle$
at
$Nt = 1.5$
, 2, 2.5 and 3. As positive and negative
$\langle u'w' \rangle$
in the
$z \gt 0$
and
$z \lt 0$
regions lead to vertical wake expansion, the negative and positive stripes in the
$z \gt 0$
and
$z \lt 0$
regions at the time instant
$Nt = 3$
necessarily suggest a narrowing process of the wake in the vertical direction. This is consistent with previously reported results on the vertical size of the wake, where the height of the wake slightly decreases beyond the near-wake regime (Brucker & Sarkar Reference Brucker and Sarkar2010; Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011; Li et al. Reference Li, Yang and Kunz2024).
Figure 4. Plots of
$ \langle u'w' \rangle$
at several early time instances: from left to right,
$Nt=1.5$
, 2, 2.5 and 3.
These data allow us to assess the Boussinesq hypothesis (the eddy-viscosity hypothesis). We compute the eddy viscosity as
\begin{equation} \nu _{t,ij} =- \frac {R_{ij}^d}{2S_{ij}}, \end{equation}
where we separately evaluate eddy viscosity for each Reynolds stress in the
$R_{ij}^d$
tensor. Here,
$R_{ij}^d$
is the deviatoric part of the Reynolds stress, and no summation of the repeated indices is implied.
Figure 5 shows the
$\nu _{t,12}$
and
$\nu _{t,13}$
components, which are relevant to the Reynolds-averaged streamwise momentum equation, at
$Nt = 1$
, 4, 30 and 110. Again, we have adjusted the ranges of the colour bars and the axes to accommodate the changing states of the flow. For comparison purposes, the colour bar ranges are kept the same between the two eddy viscosity components at the same time instant. Notice that for any given eddy-viscosity component, the plots exhibit minimal variation within the wake. This is because the variation in eddy viscosity between different components far exceeds the variation within a single component across the wake.
Figure 5. Eddy viscosity
$\nu _{t,12}$
and
$\nu _{t,13}$
at
$Nt=1$
, 4, 30, 110.
Before discussing the results in figure 5, it is important to note that the magnitude of eddy viscosity does not necessarily signify its significance in turbulence modelling. For example, in a turbulent boundary layer, the 11 component of the eddy viscosity,
$\nu _{t,11}$
, is large. However,
$\nu _{t,11}$
has little to no significance for modelling boundary layer flows. Now consider
$\nu _{t,13}$
in the stratified wake flow. The
$\langle u'w' \rangle$
stress is negligibly small in the late wake. Yet modelling
$\nu _{t,13}$
is critical for modelling stratified wakes.
Bearing the above in mind, we now examine the results in figure 5. First and foremost,
$\nu _{t,12}$
and
$\nu _{t,13}$
exhibit distinctly different behaviours, with negative eddy viscosities found in
$\nu _{t,13}$
at
$Nt = 4$
and
$Nt = 30$
. Hence the Boussinesq hypothesis, which assumes
$\nu _{t,ij} = \nu _{t,i'j'}$
, is a poor approximation of reality, necessitating more sophisticated models for stratified flows, such as Reynolds stress models. It is also interesting to observe that the eddy viscosity
$\nu _{t,12}$
remains largely constant without any complex structures. This simplification allows us to approximate the Reynolds-averaged
$\langle u \rangle$
equation in the late wake as
\begin{equation} \frac {\partial \left \lt u \right \gt }{\partial t} = \nu \left ( \frac {\partial ^2 \left \lt u \right \gt }{\partial y^2} + \frac {\partial ^2 \left \lt u \right \gt }{\partial z^2} \right ) + \nu _{t,12} \frac {\partial ^2 \left \lt u \right \gt }{\partial y^2}. \end{equation}
Here, we have neglected the vertical Reynolds stress
$\partial \langle u'w' \rangle / \partial z$
since
$\langle v \rangle$
,
$\langle w \rangle$
and
$\langle u'w' \rangle$
are small in the late wake. Equation (5.2) is, in fact, the starting point of the phenomenological model in Meunier, Diamessis & Spedding (Reference Meunier, Diamessis and Spedding2006).
5.2. Budgets of the normal Reynolds stresses
We present the budget of the normal Reynolds stresses. The discussion here covers R50F50, R20F50 and R20F02, and we explore the effects of both the Reynolds and Froude numbers. To ensure a fair comparison, we employ the same colour bar ranges for all budget terms and all three Reynolds stresses at a given time instant.
The budget terms in the transport equation of
$\langle u'u' \rangle$
for R50F50 are depicted in figure 6 at
$Nt = 1$
, 4, 30 and 110. Note that the buoyancy term does not appear in the
$\langle u'u'\rangle$
transport equation. The pressure diffusion term
$D_{p,11} = - \partial ( \langle u_1' p' \rangle \delta _{1k} + \langle u_1' p' \rangle \delta _{1k} )/\partial x_k$
is zero due to the absence of a mean gradient along the streamwise direction. From figure 6, we see that the production, dissipation, turbulent diffusion and pressure–strain correlation terms are significant. The viscous diffusion term is not significant until very late at
$Nt = 110$
. In prior budget analyses, such as plane channel flows, the residual is computed by summing the budget terms, with deviations from zero indicating convergence error. In the present study, the unsteady nature of the flow means that the sum of the budget terms yields the unsteady term (up to discretisation errors). Consequently, we must resort to other measures to assess statistical convergence. Following Li et al. (Reference Li, Yang and Kunz2024), we examine the transverse asymmetry,
\begin{equation} {\textit {err}}=\frac {\left | \int _{-2L_{Vk}}^{2L_{Vk}}\int _{-2L_{Hk}}^{0}\left \lt |\theta |\right \gt {\rm d}y{\rm d}z - \int _{-2L_{Vk}}^{2L_{Vk}}\int _{0}^{2L_{Hk}}\left \lt |\theta |\right \gt {\rm d}y{\rm d}z\right |}{{\int _{-2L_{Vk}}^{2L_{Vk}}}\int _{-2L_{Hk}}^{2L_Hk}\left \lt |\Theta | \right \gt {\rm d}y\,{\rm d}z}\times 100\,\%, \end{equation}
where
$\theta$
is an arbitrary budget term, and
$\Theta$
is the dominant term in the transport equation of the corresponding Reynolds stress. According to this measure, the error is well below 2 % for all budget terms.
Figure 6. The budget of
$\langle u'u' \rangle$
in R50F50. The black rectangles at the top-right corners of the plots in the first row indicate the true aspect ratio of the wake at the given time step.
We now closely examine the budget terms. The production term is doughnut-shaped at
$Nt = 1$
, and develops into two islands at
$Nt = 4$
. These two islands correspond to the two peaks in
$\langle u'u'\rangle$
in the late wake. The dissipation term is dominant in the early wake at
$Nt = 1$
and 4, but becomes less important in the late wake at
$Nt = 30$
and 110. Comparing figures 3 and 6, we see that the turbulent diffusion term is negative where
$\langle u'u'\rangle$
is large, and positive where
$\langle u'u'\rangle$
is small. The same is true for the viscous diffusion term at
$Nt = 110$
. The two terms differ in that the turbulent diffusion redistributes
$\langle u'u'\rangle$
in the transverse direction, whereas the viscous diffusion term redistributes
$\langle u'u'\rangle$
in the vertical direction. This is consistent with the previous argument that the expansion of the late wake in the vertical direction is a result of viscous diffusion (Meunier et al. Reference Meunier, Diamessis and Spedding2006; Brucker & Sarkar Reference Brucker and Sarkar2010). The pressure–strain correlation term acts primarily as a sink term. Comparing figures 6, 8 and 9, we see that the term takes energy from
$\langle u'u'\rangle$
, and transfers it to
$\langle w'w'\rangle$
in the early wake, and to
$\langle v'v'\rangle$
in the late wake.
Also, figure 7 shows the results from a single realisation. Unlike in figure 6, we see no coherence in the two diffusion terms or in the pressure–strain correlation term in figure 7. This confirms that ensemble averaging is critical for budget analysis.
Figure 7. From left to right: the budget of
$\langle u'u' \rangle$
at
$Nt=1$
, 4, 30 and 110 calculated from a single realisation with streamwise average only.
The budget terms in the transport equation of
$\langle v'v' \rangle$
for R50F50 are illustrated in figure 8. The buoyancy term is zero. Compared to the production in the transport equation of
$\langle u'u'\rangle$
, which is a dominant term, the production term in the transport equation of
$\langle v'v'\rangle$
is at most comparable to the other terms, leading to a faster decay of
$\langle v'v'\rangle$
than
$\langle u'u'\rangle$
in the early wake. Additionally, the term is negative in the centre region of the wake at
$Nt = 4$
, 30 and 110, which explains the negative production of turbulent kinetic energy reported in Li et al. (Reference Li, Yang and Kunz2024). This is interesting, as it suggests energy transfer from the turbulence to the mean flow through
$\langle v'v'\rangle$
. The dissipation term is significant at all four times. It is peculiar that the dissipation term develops two islands above and below
$z = 0$
. The turbulent diffusion term
$D_{T,22}$
in the transport equation of
$\langle v'v'\rangle$
closely resembles the turbulent diffusion term
$D_{T,11}$
in the transport equation of
$\langle u'u'\rangle$
at all four times. It is interesting to note that the turbulent diffusion term
$D_{T,22}$
and the pressure diffusion term
$D_{p,22}$
are almost perfectly anti-correlated, with the former taking energy from the centre of the wake, and the latter adding energy to the centre of the wake. Like the viscous diffusion term
$D_{v,11}$
in the transport equation of
$\langle u'u'\rangle$
, the viscous diffusion term
$D_{v,22}$
in the transport equation of
$\langle v'v'\rangle$
is significant only in the late wake, contributing to wake expansion in the vertical direction.
Figure 8. The budget of
$\langle v'v' \rangle$
in R50F50.
The budget of
$\langle w'w' \rangle$
is depicted in figure 9. The vertical motions are suppressed by buoyancy in the late wake, therefore the budget terms in the transport equation of
$\langle w'w'\rangle$
are significant only at
$Nt = 1$
and 4. Again, we have employed the same colour bar ranges for all the budget terms in the
$\langle u'u'\rangle$
,
$\langle v'v'\rangle$
and
$\langle w'w'\rangle$
equations at a given time. The production term is significant only at
$Nt = 4$
. Comparing the production terms in the transport equations of
$\langle v'v'\rangle$
and
$\langle w'w'\rangle$
at
$Nt = 4$
, it is interesting to notice that the two terms have the same magnitude but opposite signs. The production term measures the energy exchange between the mean flow and the turbulence, therefore energy transfer from
$\langle v'v'\rangle$
to
$\langle w'w'\rangle$
directly through the production term is not possible. Nonetheless, the data suggest that there is energy transfer from
$\langle v'v'\rangle$
to
$\langle w'w'\rangle$
through the mean flow and the production term. It is worth noting that the production terms
$P_{22}$
and
$P_{33}$
are smaller than
$P_{11}$
even in unstratified flows because
$\partial V/\partial x_i$
and
$\partial W/\partial x_i$
are small. The dissipation term is the dominant term at
$Nt = 1$
and 4, leading to a faster decay of
$\langle w'w'\rangle$
than
$\langle u'u'\rangle$
, similar to
$\langle v'v'\rangle$
. The turbulent diffusion term and the pressure diffusion term are still anti-correlated, as seen in the budget of
$\langle v'v'\rangle$
. However, the terms flip signs at
$Nt = 1$
and
$Nt = 4$
: the turbulent diffusion term is negative at the wake centre at
$Nt = 1$
, and positive at
$Nt = 4$
, whereas the pressure diffusion term is positive at the wake centre at
$Nt = 1$
, and negative at
$Nt = 4$
. The buoyancy term appears in the
$\langle w'w'\rangle$
budget but is never significant. Nonetheless, it acts as a sink term and directly suppresses the vertical motions.
Figure 9. The budget of
$\langle w'w' \rangle$
in R50F50.
The R20F50 result is similar to the R50F50 result, which is presented in Appendix A. Here, we present the budget of the normal Reynolds stresses for R20F02 at
$Nt = 4$
, 30, 110 and 400. These times approximately correspond to
$Nt = 1$
, 4, 30 and 110 in R50F50, as far as
$P_{11}$
is concerned. Figures 10, 11 and 12 show the budgets of
$\langle u'u'\rangle$
,
$\langle v'v'\rangle$
and
$\langle w'w'\rangle$
, respectively. The results are largely similar to those in figures 6, 8 and 9 for R50F50. Below, we highlight a few differences. The turbulent diffusion term in the transport equation of
$\langle u'u'\rangle$
is negative at the centre of the wake in R20F02, but is positive in R50F50. The production term in the transport equation of
$\langle v'v'\rangle$
is insignificant at all four times in R20F02, whereas in R50F50, the term is comparable to other terms. Unlike in R50F50, where the buoyancy term is never significant,
$B_{33}$
in the transport equation of
$\langle w'w'\rangle$
is the most significant term in the early wake in R20F02. It is noteworthy that although the buoyancy term plotted here acts as a sink, it can function as both a sink and a source at different time snapshots, as shown in figure 13.
Figure 10. The budget of
$\langle u'u' \rangle$
in R20F02.
Figure 11. The budget of
$\langle v'v' \rangle$
in R20F02.
Figure 12. The budget of
$\langle w'w' \rangle$
in R20F02.
Figure 13. The buoyancy term in R20F02.
Figures 14 and 15 summarise the budget of the three normal Reynolds stresses in the early and late wakes at
$Nt = 1$
and 110, respectively. In the early wake,
$\langle u'u'\rangle$
extracts energy from the mean flow through
$P_{11}$
, dissipates the energy to heat through
$\epsilon _{11}$
, and transfers energy to
$\langle w'w'\rangle$
via
$\phi _{11}$
. In addition, the turbulent diffusion term redistributes energy from the centre of the wake to the periphery;
$\langle v'v'\rangle$
dissipates energy to heat through
$\epsilon _{22}$
. The pressure diffusion term and the turbulent diffusion term redistribute energy, with the former moving energy to the centre of the wake, and the latter taking energy away from the centre of the wake. Also,
$\langle w'w'\rangle$
receives energy from
$\langle u'u'\rangle$
through the pressure–strain correlation, dissipates energy to heat through
$\epsilon _{33}$
, and interchanges with potential energy through
$B_{33}$
, with the pressure diffusion term and the turbulent diffusion term moving energy to and away from the centre of the wake. In the late wake,
$\langle u'u'\rangle$
extracts energy from the mean flow through
$P_{11}$
, dissipates to heat through
$\epsilon _{11}$
, and transfers energy to
$\langle v'v'\rangle$
via
$\phi _{11}$
. The turbulent diffusion term and the viscous diffusion term take energy away from the two peaks in
$\langle u'u'\rangle$
. We have that
$\langle v'v'\rangle$
dissipates to internal energy through
$\epsilon _{22}$
, and transfers energy to the mean flow through
$P_{22}$
. The viscous diffusion and the turbulent diffusion terms take energy away from the centre of the wake, while the pressure diffusion term moves energy to the centre of the wake. The
$\langle w'w'\rangle$
term is negligibly small. There is some energy transfer from
$\langle u'u'\rangle$
, with dissipation taking the energy away.
Figure 14. A summary of the Reynolds normal stress budget in R50F50 at
$Nt=1$
: (a)
$\langle u'u'\rangle$
, (b)
$\langle v'v'\rangle$
, (c)
$\langle w'w'\rangle$
.
Figure 15. Same as figure 14 but for
$Nt=110$
.
5.3. The budget of Reynolds shear stress
The discussion of the Reynolds shear stress will focus on
$\langle u'v' \rangle$
and
$\langle u'w' \rangle$
. These two terms appear in the Reynolds-averaged streamwise momentum equation. The budget of
$\langle u'v' \rangle$
is presented in figure 16 at
$Nt = 1$
, 4, 30 and 110. The significant terms in the transport equation of
$\langle u'v' \rangle$
are the production, turbulent diffusion, pressure diffusion and pressure–strain correlation. The buoyancy term does not appear in the equation. Dissipation and viscous diffusion are present but small. We note the following about the dissipation and viscous diffusion terms. The dissipation term is small both here and in the transport equations of
$\langle u'w'\rangle$
and
$\langle v'w'\rangle$
. Hence while isotropy is not a valid assumption for the dissipation tensor, it is valid, at least to a good approximation, to model
$\epsilon _{ij}$
as diagonal. The viscous term is also small, not only here, but also in the transport equations of
$\langle u'w'\rangle$
and
$\langle v'w'\rangle$
, even in the late wake. This is interesting as viscous diffusion is significant in the transport equations of the normal stresses in the late wake.
We proceed by examining the terms in figure 16. First, we notice that the terms are symmetric with respect to
$z/D = 0$
, and anti-symmetric with respect to
$y/D = 0$
. Comparing figures 3 and 16, we observe that due to the anti-symmetric distribution of
$\langle u'v'\rangle$
with respect to
$y/D = 0$
, the production term is always a source term and never a sink. The turbulent diffusion term
$D_{T,12}$
is mostly in the direction of the mean gradient. This term and the pressure diffusion term have similar spatial distributions. Nonetheless, they are anti-correlated in the early wake, and positively correlated in the late wake, with the transition occurring approximately at
$Nt \approx 4$
. The pressure–strain correlation term acts as a sink term here, balancing the production term.
The significant terms in the transport equation of
$\langle u'w' \rangle$
are plotted in figure 17. Here, the terms are symmetric with respect to
$y/D = 0$
, and anti-symmetric with respect to
$z/D = 0$
. In the early wake, the significant terms are the production, turbulent diffusion, pressure diffusion, pressure–strain correlation and buoyancy. As in the transport equation of
$\langle u'v'\rangle$
, the production term is always a source and never a sink. Similarly, the pressure diffusion and turbulent diffusion are anti-correlated in the early wake, and the pressure–strain correlation acts as a sink. The buoyancy term also acts as a sink. In the late wake, the significant terms are the pressure diffusion, pressure–strain correlation and buoyancy. The production and turbulent diffusion terms are small as buoyancy suppresses the vertical turbulent motions. Among the more significant terms, the combination of pressure–strain correlation and buoyancy approximately balances with the pressure diffusion term. In addition, the terms seem to reflect the presence of gravitational waves in the late wake. Finally, we note that the colour bar range at
$Nt = 110$
in figure 17 is an order of magnitude larger than that in figure 16.
Figure 16. The budget of
$\langle u'v' \rangle$
in R50F50.
Figure 17. The budget of
$\langle u'w' \rangle$
in R50F50.
We summarise the budgets of the Reynolds shear stress in the early and late wakes in figures 18 and 19. In the early wake, both
$\langle u'v' \rangle$
and
$\langle u'w' \rangle$
extract energy from the mean flow through
$P_{12}$
and
$P_{13}$
, respectively, and transfer energy to the pressure fluctuation through
$\phi _{12}$
and
$\phi _{13}$
. Buoyancy enables the transfer between
$\langle u'w' \rangle$
and potential energy through
$B_{13}$
. The turbulent diffusion and pressure diffusion terms move energy away from and towards the regions where the two shear stresses are large. In the late wake,
$\langle u'v' \rangle$
still extracts energy from the mean flow through
$P_{12}$
and transfers energy to the pressure fluctuation through
$\phi _{12}$
. Both the turbulent diffusion and pressure diffusion terms move energy away from the regions where
$\langle u'v'\rangle$
is large.
Figure 18. A summary of the Reynolds shear stress budget in R50F50 at
$Nt=1$
(early wake).
Figure 19. A summary of the Reynolds shear stress budget in R50F50 at
$Nt=110$
(late wake).
5.4. The modelling of the pressure–strain correlation term
The ensemble-averaged data allow us to assess the existing closure models. In § 2, we demonstrate that eddy-viscosity type closures are unsuitable for stratified wake flows, necessitating the use of Reynolds stress models. In the context of Reynolds stress models, the dissipation term, the diffusion terms and the pressure–strain correlation term must be closed. Jain et al. (Reference Jain, Huang, Yang and Kunz2022a , Reference Jain, Huang, Li, Yang and Kunz2023) examined the dissipation closures and indicated that isotropic dissipation is not valid. Our data are consistent with these findings. Furthermore, the data suggest that it is valid to model the dissipation tensor as a diagonal one. Daly & Harlow (Reference Daly and Harlow1970) model the sum of the pressure, viscous and turbulent diffusion terms as diffusion terms. This is mostly valid as the sum of the three terms aligns with the mean gradient direction, at least according to the data presented here. In this subsection, we closely examine the pressure–strain correlation term and the closure in Speziale et al. (Reference Speziale, Sarkar and Gatski1991), which we will refer to as the SSG model.
Figures 20 and 21 compare the SSG model and DNS results at
$Nt = 1$
(early wake) and
$Nt = 110$
(late wake) in R50F50 for the deviatoric stress components
$\phi _{12}$
,
$\phi _{13}$
and
$\phi _{23}$
. The SSG model accurately predicts
$\phi _{12}$
and
$\phi _{13}$
at
$Nt = 1$
when given the correct inputs. However, for
$\phi _{23}$
, it is approximately zero in the DNS data, but the SSG model predicts structures that are not observed in the DNS results. The performance of the model deteriorates at
$Nt = 110$
. Although the model still predicts the magnitude of the terms, it fails to capture the spatial coherence. Specifically, the SSG-predicted
$\phi _{23}$
and the DNS data have opposite signs. The poor performance of the model in the late wake is not unexpected, as the model assumes near-equilibrium states, whereas the late wake is far from equilibrium.
Figure 20. The pressure–strain correlation of the deviatoric stress at
$Nt = 1$
. The top row shows the prediction of the SSG model. The bottom row shows the DNS. We use the same colour bar for the same term.
Figure 21. The pressure–strain correlation of the deviatoric stress at
$Nt = 110$
.
Figures 22 and 23 compare the SSG model and DNS results at
$Nt = 1$
(early wake) and
$Nt = 110$
(late wake) in R50F50 for the normal components
$\phi _{11}$
,
$\phi _{22}$
and
$\phi _{33}$
. While the SSG model correctly captures
$\phi _{11}$
as a major sink and
$\phi _{33}$
as a major source in the early wake, the spatial distributions of these terms are not accurately predicted. Furthermore, the SSG model predicts some coherent structure for
$\phi _{22}$
, whereas it is approximately zero in DNS. In the late wake, the model still correctly predicts
$\phi _{11}$
as a sink. However, instead of redistributing from
$\langle u'u'\rangle$
to
$\langle v'v'\rangle$
as observed in the DNS results, the SSG model redistributes from
$\langle u'u'\rangle$
and
$\langle v'v'\rangle$
to
$\langle w'w'\rangle$
. Finally, SSG’s prediction of
$\phi _{33}$
is almost an order-of-magnitude higher value than that in DNS.
Figure 22. The pressure–strain correlation of the normal stress at
$Nt = 1$
.
Figure 23. The pressure–strain correlation of the normal stress at
$Nt = 110$
.
6. Concluding remarks
This paper presents a detailed analysis of the Reynolds stresses and their budgets in temporally evolving stratified wake flows, derived from ensembles of DNS. By leveraging ensemble averaging, we mitigate statistical errors, enabling robust identification of the dominant terms in the Reynolds stress transport equations as well as their spatial structures. Our findings elucidate the processes governing the generation, dissipation and redistribution of Reynolds stresses, while shedding light on the influence of the Reynolds and Froude numbers. Additionally, we assess the validity of existing Reynolds stress models, revealing critical insights into their performance across different wake regimes.
Key findings from the budget analysis include the following. In the late wake, the streamwise Reynolds stress component
$ \langle u'u' \rangle$
exhibits two distinct peaks, driven by corresponding peaks in its production term. Viscous diffusion remains negligible until the late wake, where it becomes significant in the normal stress transport equations but continues to play a minor role in the shear stress budgets. The pressure–strain correlation redistributes energy from
$ \langle u'u' \rangle$
to
$ \langle w'w' \rangle$
in the early wake, and to
$ \langle v'v' \rangle$
in the late wake. A strong correlation between the pressure and turbulent diffusion terms is observed, despite no apparent direct connection between them. Regarding the modelling of the Reynolds stresses, while the validity of the Boussinesq hypothesis is generally questioned, the assumption of constant eddy viscosity
$ \nu _{t,12}$
appears reasonable. Modelling of the unclosed terms in the Reynolds stress transport equations is also examined. The sum of pressure, viscous and turbulent diffusion terms can be treated collectively as a diffusion term; however, neglecting the pressure diffusion component introduces significant errors. Additionally, while isotropy is not a valid assumption for the dissipation tensor, the data support modelling it as diagonal. Finally, the SSG model performs well in the early wake for
$ \langle \phi _{11} \rangle$
and
$ \langle \phi _{13} \rangle$
, but deteriorates for
$ \langle \phi _{23} \rangle$
and in the late wake.
The findings in this study underscore the need for continued investigation of the Reynolds stress budgets. In particular, the Reynolds number range investigated here is limited. Higher Reynolds numbers will lead to more prominent layering of shear and stronger emission of internal waves. These effects are known to influence Reynolds stress budgets, therefore further analyses will be necessary when higher Reynolds number data become available. Meanwhile, this study focuses on the budget analysis itself. Linking these budget terms to flow phenomena/structures, such as internal gravity waves, is also necessary.
Acknowledgements.
This work was sponsored by the US Office of Naval Research under contracts N000142012315, N0001419012057 and N000142312695, with Dr P. Chang as Technical Monitor.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Further results
In this appendix, we present the Reynolds stress budgets for R20F50 and R20F02. The results at these additional conditions are largely consistent with those presented in § 5, allowing us to extrapolate the observations – at least within the Reynolds number range investigated here.
Figures 24, 25 and 26 present the budgets of
$\langle u'u'\rangle$
,
$\langle v'v'\rangle$
and
$\langle w'w'\rangle$
for R20F50. The results are largely similar to those presented in figures 6, 8 and 9 for R50F50. Here, we highlight the differences. Unlike R50F50, the production term in the transport equation of
$\langle v'v'\rangle$
is positive at the wake centre at
$Nt=110$
.
Figure 24. The budget of
$\langle u'u' \rangle$
of R20F50.
Figure 25. The budget of
$\langle v'v' \rangle$
of R20F50.
Figure 26. The budget of
$\langle w'w' \rangle$
of R20F50.
Figure 27 and figure 28 present the budget terms in the transport equation of
$\langle u'u'\rangle$
for R20F50 and R20F02, respectively. Similar to the observations for case R50F50, the production term acts exclusively as a source term and never as a sink. For all cases and all time steps presented, the budget of
$\langle u'u'\rangle$
is primarily governed by the production, turbulent diffusion, pressure diffusion, and pressure-strain correlation, where the combined contributions of the turbulent diffusion, pressure diffusion, and pressure-strain correlation nearly balance the production term.
Figure 27. The budget of
$\langle u'v'\rangle$
of R20F50.
Figure 28. The budget of
$\langle u'v'\rangle$
of R20F02.
The significant terms in the transport equation of
$\langle u'w' \rangle$
are plotted in figures 29 and 30 for cases R20F50 and R20F02, respectively. While most budget terms are similar to the results of R50F50, a clear Froude number dependence of the buoyancy term emerges when comparing figure 30 with figures 17 and 29. At
$Nt=4$
, instead of being comparable to the production term, the buoyancy term of R20F02 is approximately an order of magnitude higher than the other terms, acting as a sink without stripe structures. In the late wake, for all cases presented here, the significant terms are the pressure diffusion, pressure–strain correlation and buoyancy. Again, the combination of the pressure–strain correlation and buoyancy approximately balances with the pressure diffusion term.
Figure 29. The budget of
$\langle u'w' \rangle$
of R20F50.
Figure 30. The budget of
$\langle u'w' \rangle$
of R20F02.
Open access
