3.1. Instantaneous wake behaviour

Figure 2(a–c) illustrates the instantaneous defect velocity ( $u_d = U_\infty -u_x$ ) at $x/D=3$ for R90F1.9 at three different time instances. The near wake exhibits two vertically elongated, ear-like lobes positioned atop a central structure. Specifically, at $tU_\infty /D=109.34$ (figure 2 a), both lobes are situated close to the centre-plane ( $y=0$ ). Subsequently, at $tU_\infty /D=109.76$ (figure 2 b), the lobes diverge in the spanwise direction, moving away from the centre-plane. Finally, at $tU_\infty /D=110.152$ (figure 2 c), the lobes converge, returning to a position closest to the centre-plane. This sequence demonstrates a spanwise flapping motion of the wake, characterised by a distinct temporal frequency, as will be shown later. The selected timestamps are chosen to accurately capture a complete cycle of this flapping motion. Notably, the wake structure maintains a configuration with predominantly lateral symmetry throughout the oscillation, resembling varicose modes with symmetric, in-phase deformations that have been found previously in shear flows, e.g. shear layers at the edge of a plane jet (Stanley & Sarkar Reference Stanley and Sarkar1997) and the wake of a cylindrical roughness element on a flat plate (Bucci et al. Reference Bucci, Puckert, Andriano, Loiseau, Cherubini, Robinet and Rist2018). For a more dynamic visualisation of this flapping motion, please refer to the supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10923.

At ${\textit{Fr}}=6$ , the wake exhibits a structural morphology qualitatively similar to that observed at ${\textit{Fr}}=1.9$ , characterised by two lobes connected to a central structure (figure 2 d–f and supplementary movie 2). A notable dissimilarity, however, lies in the reduced elongation of the lobes compared with R90F1.9. Nevertheless, an analogous flapping motion is evident, wherein the lobes undergo continuous inward and outward oscillations, maintaining a symmetric configuration.

The temporal evolution of the unstratified case (R90F $\infty$ ) deviates somewhat from the previously described stratified cases. Figure 2(g–i) depicts $u_d$ contours at $x/D=3$ for three different time instances. While a double-lobed structure remains discernible, the lobes are considerably less differentiated from the central structure. Moreover, the instantaneous wake structure exhibits asymmetry about the centre-plane ( $y=0$ ). At $tU_\infty /D=72.24$ (figure 2 g), the left half of the wake displays a larger extent compared with the right half. Conversely, this polarity of asymmetry reverses as the simulation progresses to $tU_\infty /D=84.76$ (figure 2 h), where the right half manifests a larger extent. Subsequently, the asymmetry oscillates back to the left half at $tU_\infty /D=98.24$ (figure 2 i). These timestamps were carefully selected to fully capture a complete cycle of these oscillations. Based on the timestamps, it can be inferred that this oscillating asymmetry is characterised by a low frequency. The unstratified wake also exhibits a high-frequency spanwise flapping of the double lobes, analogous to the stratified cases, as demonstrated in supplementary movie 3. Consequently, both the oscillating asymmetry (low frequency) and flapping (high frequency) modes are dominant in the R90F $\infty$ flow, and these modes will be discussed in subsequent sections.

To gain further insight into the high frequency flapping mode, contours of instantaneous vertical vorticity ( $\omega _z$ ) for the R90F1.9 wake, chosen as an example, are shown on the horizontal centre-plane ( $z=0$ ) in figure 3. Shear layer separation is evident on each side of the body at $x/D\approx 1.5{-}2.5$ . These shear layers later become unstable and rollup into Kelvin–Helmholtz (K–H) billows after $x/D\approx 2.5$ before breaking up into small-scale structures further downstream. Looking closely at $x/D=3$ at $tU_\infty /D=109.34$ , the K–H rollups have not yet reached that streamwise location and, thus, the $\omega _z$ imprint is the narrowest in the spanwise direction. This is the same time at which the $u_d$ double-lobed structure observed in figure 2(a) is the closest to the centre-plane ( $y=0$ ). Furthermore, at $tU_\infty /D=109.76$ , the K–H billows are somewhat centred around $x/D=3$ giving the widest spanwise imprint of $\omega _z$ . This timestamp also corresponds to the diverging double-lobed structure in figure 2(b). Therefore, it can be surmised that the aforementioned flapping mode is related to the shear layers becoming unstable after leaving the inclined spheroid. For a more dynamic visualisation of the shear layer instability, please refer to supplementary movie 4.

Figure 3. Instantaneous vertical vorticity ( $\omega _z$ ) contours plotted on the centre horizontal plane ( $z=0$ ) for R90F1.9 at three distinct timestamps similar to figure 2(a–c). The separated shear layers exhibit Kelvin–Helmholtz billows.

3.2. Point spectra

The spanwise fluctuation velocity ( $u^{\prime}_y$ ) is obtained from the computed $u^{\prime}_r, u^{\prime}_\theta$ components. Analysis of $u^{\prime}_y (t)$ with pointwise spectra in this subsection and $u^{\prime}_y (r, \theta , t)$ with SPOD spectra in the following subsection reveal dominant modes as elaborated here.

Figure 4. Point spectra of spanwise velocity fluctuation ( $u^{\prime}_y$ ) for (a) R90F1.9, (b) R90F6 and (c,d) R90F $\infty$ at different locations in the wake, where panel (d) is computed for a longer time series. Location of points A, B and C for each $\textit{Fr}$ case is marked atop the $u_d$ contours shown for $x/D=3$ . At $x/D=10$ , spectra shown in dashed red represent a point at the edge of the wake and dashed green represent a point on the vertical centre-plane.

Figure 4(a) shows temporal point spectra of the spanwise fluctuation velocity ( $u^{\prime}_y$ ) for R90F1.9 at different points on the cross-stream planes at $x/D=3$ and $10$ . At $x/D=3$ , point spectra show a peak at $St \approx 1.15$ for points at $\theta =0^\circ \text{(point A) and } 180^\circ$ (point B) and radial distance of $r/D=0.2$ from centre. Points A and B are located in the right and left shear layers. The spectrum at point C ( $\theta =270^\circ , r/D=0.5$ ) in the central structure at $x/D=3$ does not show any prominent peak. Hence, not much fluctuation activity is present in the central wake structure and the $St \approx 1.15$ peak is the dominant wake mode, as will be more rigorously established with SPOD analysis. At a downstream location of $x/D=10$ , the spectra peak at $St\approx 1$ for point $\theta =0^\circ \text{ and } r/D=0.2$ at the wake edge. However, this peak is weaker compared with that at $x/D=3$ , suggesting streamwise evolution of the coherent shear-layer mode. The spectra at further downstream locations show progressive diminishing of this peak until it disappears at $x/D \approx 20$ . The spectrum at a point on the centre-plane ( $\theta =270^\circ , r/D=0.2$ ) does not show any dominant frequency at $x/D=10$ .

Point spectra of $u^{\prime}_y$ for R90F6 (figure 4 b) show characteristics similar to R90F1.9. A peak at $St\approx 1.25$ is observed at $x/D=3$ for locations A ( $r/D=0.4,\theta =315^\circ$ ) and B ( $r/D=0.4,\theta =225^\circ$ ). Again, no peak is observed in the spectrum at point C ( $r/D=0.6,\theta =270^\circ$ ) in the centre-plane. At $x/D=10$ , a local peak is observed at $St\approx 0.75$ for a point at the wake edge ( $r/D=0.75,\theta =300^\circ$ ), whereas no dominant peak is observed in the vertical centre-plane ( $r/D=1,\theta =270^\circ$ ).

Spectra of $u^{\prime}_y$ for the R90F $\infty$ wake at $x/D=3$ show a peak at $St\approx 1.28$ (figure 4 c) for points A ( $r/D=0.5,\theta =315^\circ$ ) and B ( $r/D=0.5,\theta =225^\circ$ ) on the double-lobed structure, see figure 2(g–i). This peak corresponds to the flapping mode. A significant power spectral density (PSD) level is also observed for low frequencies $St\sim O(0.01)$ at point C ( $r/D=0.5,\theta =270^\circ$ ) on the centre-plane. Taking a longer time series for the point spectra helped resolve those low frequencies, see figure 4(d). A peak is observed at $St\approx 0.04$ , which corresponds to the slow oscillating asymmetry exhibited by the unstratified wake, discussed in § 3.1. At $x/D = 10$ and beyond, no prominent peak is seen in the power spectra. Thus, shear-layer flapping and oscillating asymmetric modes are confined to the near wake.

3.3. SPOD eigenspectra and leading eigenmenodes

SPOD analysis was conducted for the flow on a 2-D plane at $x/D=3$ for the R90 series. For the unstratified case R90F $\infty$ , the analysed flow field spans $0\leqslant r/D \leqslant 10$ , comprising a total of $N_r^{\textit{SPOD}}=678$ radial grid points and $N_\theta =256$ points in the azimuthal direction. Stratified cases R90F6 and R90F1.9 contain $N_r^{\textit{SPOD}}=940$ grid points over $0\leqslant r/D \leqslant 11.5$ , and $N_\theta =128$ points in the azimuthal direction. Table 2 lists the parameters for the SPOD analysis across various $\textit{Fr}$ cases at ${\textit{Re}}=9\times 10^4$ . In this context, $N$ denotes the total number of snapshots, spaced by $\Delta t D/U_\infty$ between consecutive snapshots; $N_{\textit{freq}}$ represents the number of frequencies that are being resolved. The $N$ snapshots are divided into discrete sets of overlapping blocks with an overlap of $N_{\textit{ovlp}}$ snapshots. A total of $N_{\textit{blk}} = (N - N_{\textit{ovlp}})/(N_{\textit{freq}}-N_{\textit{ovlp}})$ SPOD modes are extracted at each frequency.

Table 2. SPOD parameters for ${\textit{Re}}=9\times 10^4$ wakes

Figure 5. (a) SPOD eigenspectra, $\lambda ^{(1)}, \lambda ^{(2)},\ldots , \lambda ^{(n)}$ , for all eigenmodes of R90F1.9 at a streamwise location of $x/D=3$ . (b) Real part of the leading SPOD mode for the spanwise velocity, ${\varPhi}_{v}^{(1)}(y,z,St;x/D)$ , corresponding to the peak in the eigenspectrum of $\lambda ^{(1)}$ .

Figure 5(a) presents the SPOD eigenspectra of the R90F1.9 wake at $x/D=3$ . The leading eigenmode ( $\lambda ^{(1)}$ ) is highlighted in blue. A prominent high-frequency peak is evident in the eigenspectrum of $\lambda ^{(1)}$ at $St\approx 1.15$ . Notably, the frequency of $St\approx 1.15$ is rarely associated with a vortex shedding (VS) mode in the literature, which typically has a frequency $St\sim \textit{O}(0.1)$ (Berger, Scholz & Schumm Reference Berger, Scholz and Schumm1990; El Khoury, Andersson & Pettersen Reference El Khoury, Andersson and Pettersen2012; Ohh & Spedding Reference Ohh and Spedding2024). To delve deeper into the nature of this high-frequency mode, the real part of the normalised (by $L_\infty$ norm) leading SPOD mode, ${\varPhi}_v^{(1)}(y,z,St;x/D)/{\Vert}{\varPhi}_v^{(1)}(y,z,St;x/D){\Vert}_\infty$ , where ${\varPhi}_v$ refers to the spanwise velocity $u_y$ , is presented in figure 5(b). From the plot, two distinct coherent regions/patches of opposite sign are discernible, situated on either side of the centre-plane ( $y=0$ ). The sign of these patches ensures that they are in opposition, meaning that if the right patch exhibits a positive spanwise velocity, the corresponding left patch will have a negative spanwise velocity of comparable magnitude. Consequently, both coherent patches simultaneously move outward or inward relative to the centre-plane. Furthermore, a striking similarity emerges between figures 2(a–c) and 5(b), as the oscillating portion of the wake aligns with the coherent patches in the eigenmode. Moreover, the size of the eigenmode coherent patches is comparable to the swept area of the oscillating double-lobes in the instantaneous wake. Therefore, this eigenmode corresponds to the flapping motion of the shear layers observed in the instantaneous defect velocity field ( $u_d$ ). It is worth noting that high-frequency modes with $St \sim O(1)$ have been identified in the near wake of a disk and attributed to shear layer instability (Berger et al. Reference Berger, Scholz and Schumm1990; Nidhan et al. Reference Nidhan, Chongsiripinyo, Schmidt and Sarkar2020).

Figure 6. Same as figure 5 for R90F6.

Similarly, the SPOD eigenspectra for all eigenmodes of the R90F6 wake at $x/D=3$ are depicted in figure 6(a). A prominent high-frequency peak at $St\approx 1.21$ is evident in the leading mode ( $\lambda ^{(1)}$ ). Furthermore, examining the shape of the leading SPOD mode for spanwise velocity ( ${\varPhi}_v^{(1)}(y,z,St;x/D)$ ) at the peak frequency (figure 6 b), it becomes evident that two patches of opposite signs are present in a spanwise symmetric manner, reminiscent of ${\textit{Fr}}=1.9$ . By juxtaposing $u_d$ contours in figure 2(d–f) with the leading SPOD mode, it can be deduced that the leading mode corresponds to the high-frequency flapping of the two lobes/shear layers that were observed in § 3.1.

Figure 7. Same as figure 5 for R90F $\infty$ . Here, panel (b) corresponds to the leading eigenmode at $St\approx 1.28$ peak and panel (c) corresponds to eigenmode of low-frequency mode at $St\approx 0.04$ .

Eigenspectra for R90F $\infty$ differ somewhat from the stratified cases, as illustrated in figure 7(a). While a high-frequency local maximum is observed for $\lambda ^{(1)}$ at $St\approx 1.28$ , similar to the aforementioned cases, a substantial amount of energy persists in the low-frequency modes for $St\lt 0.1$ . Similar to point spectra (figure 4 d), eigenspectra also show a peak at $St\approx 0.04$ , therefore, the eigenmode for this frequency was also diagnosed along with that for the high-frequency peak. Figures 7(b) and 7(c) depict the shape of the leading eigenmodes for spanwise velocity ( ${\varPhi}_v^{(1)}(y,z,St;x/D)$ ) at $St\approx 1.28$ (high-frequency peak) and $St\approx 0.04$ (low-frequency mode), respectively. At $St\approx 1.28$ , the leading eigenmode again exhibits a two-patch structure with a sign change between $y \gt 0$ and $y\lt 0$ . Comparing the shape and location of the patches with $u_d$ contours in figure 2(g–i), it can be inferred that this mode represents the flapping mode of the double-lobed structure. Notably, the low-frequency mode of $St\approx 0.04$ exhibits a prominent coherent structure on the centre-plane ( $y=0$ ). The sign of the $St \approx 0.04$ mode remains consistent across the structure, implying that the spanwise velocity in this region is uniformly in the same direction. This mode corresponds to the oscillating asymmetric mode, which modifies the size of both lobes simultaneously, increasing one and decreasing the other.

It is noted that there is significant energy at the low-frequency cutoff of the spectrum in figure 7(a). The capture of such a very low frequency (VLF) mode, which is often seen in laboratory wakes, presents a well-known challenge to numerical simulations because of the necessity of a very long simulation time. For example, Zhu & Morrison (Reference Zhu and Morrison2021) found increased energy at low frequencies ( $St\lt 0.01$ ) in the SPOD eigenspectra for the $m=\pm 1$ mode, but were unable to find a specific VLF peak because of the insufficient length of the time series. They inferred VLF dynamics – a multistable wake that exhibits a slow rotation of the vortex shedding plane (Rigas et al. Reference Rigas, Oxlade, Morgans and Morrison2014). Investigation of the VLF mode is beyond the scope of this paper.

3.4. Oscillating force on the body

Dominant modes of the flow, presented in previous subsections, impart significant forces on the wake generator. Force coefficients in this section are defined as $C_i=F_i/(1/2\rho U^2_\infty A)$ , where $A=\pi D^2/4$ and $F_i$ is the net force on the body in the $i$ -direction. Figure 8(a) shows the time evolution of the net spanwise force coefficient ( $C_y$ ) on the body for R90F $\infty$ . A cyclic forcing can be observed with an amplitude of $|C_y|\approx 0.025$ with a mean value $\langle C_y\rangle \approx 0$ . The amplitude of the sideways force is significant, ${\sim} 10\,\%$ of the mean drag ( $\langle C_x\rangle =0.238$ ) and lift ( $\langle C_z\rangle =0.264$ ). The sideways force (harmonic oscillation superposed on a mean) is significantly different from the ${\textit{Re}}=3\times 10^4$ wake, where the sideways force has a non-zero mean and no temporal variability.

Figure 8. (a) Time evolution of total spanwise force on the body, (b) Spectra of total force in the spanwise, vertical and streamwise directions for R90F $\infty$ .

The spectrum deduced from the $C_y$ time series exhibits a dominant peak at $St\approx 0.04$ (figure 8 b). This is the same frequency as the leading eigenmode observed for the asymmetric oscillatory mode. This clearly shows that the flow field imparts a cyclic spanwise force on the body through its leading eigenmode. Furthermore, no high-frequency peak is observed for any of the force coefficients. Since the high-frequency SPOD mode at $St \approx 1.28$ had two patches of opposite phase on either side of the body, their force contributions cancel out and this oscillatory mode imparts no net force to the body.

Since the R90F1.9 and R90F6 wakes are found to be spanwise symmetric at every time instant, no net oscillating spanwise force is observed. Also, no dominant peaks are exhibited in the spectra (not shown for brevity) of the time-varying force.

4.1. Mean defect velocity ( $U_d$ )

Figure 9. Contours of mean defect velocity ( $U_d$ ) for two streamwise locations (a,b,e, f,i,j) $x/D=3$ and (c,d,g,h,k,l) $x/D=10$ are shown for (a–d) ${\textit{Fr}}=\infty$ , (e–h) ${\textit{Fr}}=6$ and (i–l) ${\textit{Fr}}=1.9$ at (a,c,e,g,i,k) ${\textit{Re}}=9\times 10^4$ and (b,d, f,h,j,l) ${\textit{Re}}=3\times 10^4$ (adapted from NJOS25). Isopycnals are overlaid on the plots (red) for stratified cases. Radial extent ( $r/D$ ) of domain is $r/D=1$ unless explicitly mentioned.

The wake is characterised by a momentum deficit relative to an inbound flow and quantification of the mean defect velocity, defined as $U_d = U_\infty - U_x$ , is of importance. Figure 9 displays $U_d$ contours at two streamwise locations $x/D=3\text{ and }10$ for ${\textit{Fr}}=\infty$ , $6$ and $1.9$ wakes at ${\textit{Re}}=9\times 10^4$ and ${\textit{Re}}=3\times 10^4$ (adapted from NJOS25). For the stratified cases, isopycnals (lines of constant density) are superimposed on the $U_d$ contours.

In the post-separation region at $x/D=3$ , the time-averaged wakes for all $\textit{Fr}$ cases at ${\textit{Re}}=9\times 10^4$ exhibit symmetry about the centre-plane ( $y=0$ ), as depicted in figures 9(a), 9(e) and 9(i). It is worth emphasising that, although the instantaneous wake of R90F $\infty$ was found to be laterally asymmetric (figure 2 g–i), the mean stays symmetric at all streamwise locations. This behaviour contrasts sharply with the asymmetric wakes observed at R30F $\infty$ and R30F6 (figure 9 b, f). Furthermore, all wakes at ${\textit{Re}}=9\times 10^4$ display a double-lobed structure with a lobe emanating from each lateral side of the body. For R90F $\infty$ , the lobes are well rounded (figure 9 a). In contrast, the stratified cases R90F6 and R90F1.9 demonstrate elongation of the lobes in the vertical direction, with this stretching increasing as stratification intensifies (figure 9 e,i). Thus, similar to the lower- ${\textit{Re}}$ case of NJOS25, buoyancy inhibits rollup and streamwise vortex formation. At the same $x/D$ location, significant distortion of the isopycnals is also observed within the wake for both stratified cases. The wake edges are bounded by isopycnals.

At downstream locations of $x/D=10$ , the double-lobed structure appears to be preserved for the R90F $\infty$ case, as shown in figure 9(c). Two distinct local maxima for $U_d$ are also observed at these locations, similar to R30F $\infty$ (figure 9 d). These lobes constitute a streamwise vortex pair, discussed extensively by NJOS25. Furthermore, higher vertically downward drift of the wake structure is observed for R90F $\infty$ as compared with R30F $\infty$ . The stratified cases R90F6 and R90F1.9 exhibit a loss of the double-lobed structures by $x/D=10$ , as shown in figures 9(g) and 9(k). Both stratified wakes show a slight increase of the width relative to the height, demonstrating the anisotropic effects of stratification. The isopycnals tend to flatten out after $x/D=3$ for both cases. For ${\textit{Fr}}=6$ , the flattening process is slower compared with ${\textit{Fr}}=1.9$ due to weaker buoyancy, as R90F6 exhibits slightly higher distortions of isopycnals at $x/D=10$ . A key ${\textit{Re}}$ effect is that the wakes at ${\textit{Re}}=9\times 10^4$ are more mixed at $x/D =10$ relative to ${\textit{Re}}=3\times 10^4$ as is reflected by the spatial contours. Due to this increased mixing, the $U_d$ magnitude also decays faster at higher ${\textit{Re}}$ , as shown by the difference in contour levels at $x/D=10$ .

Figure 10. Contours of mean defect velocity ( $U_d$ ) for R90F1 at four streamwise locations: (a) $x/D=3$ ; (b) $x/D=5$ ; (c) $x/D=10$ and (d) $x/D=20$ . Isopycnals are overlaid on the plots (red). Radial extent of domain is $r/D=1$ for each panel.

At the strong stratification of ${\textit{Fr}}=1$ , buoyancy significantly alters the wake structure and dynamics. NJOS25 observed two distinct patches in the mean defect velocity ( $U_d$ ) contours for the R30F1 wake (see their figure 11) near separation ( $x/D=3$ ), with a smaller secondary defect region present above the larger primary region. This structure, characterised by two local maxima, evolved downstream into distinct layers. Figure 10 displays $U_d$ contours for the R90F1 wake at four streamwise locations: $x/D=3, 5, 10$ and $20$ . At $x/D=3$ , two distinct $U_d$ patches are observed, similar to R30F1. At both values of ${\textit{Re}}$ , the edges of the $U_d$ contours appear to align consistently with the isopycnals throughout the domain, resembling a jigsaw fit. It can be inferred that buoyancy strongly influences the flow topology in the vertical plane leading to a secondary wake above a primary wake. Progressing further downstream, this double-layered structure of the wake is maintained throughout the domain. Overall, the $\textit{Fr} = 1$ wake remains double-layered at ${\textit{Re}}=9\times 10^4$ although the boundaries of the two layers are not as sharp as at ${\textit{Re}}=3\times 10^4$ . Both R90F1 and R30F1 wakes are found to be transitional throughout the domain, i.e. they do not become fully turbulent, and therefore, are not discussed further in this study.

4.2. Vertical velocity for spheroid at pitch

When a slender object, a prolate spheroid in this study, is placed at a moderate pitch angle with respect to the oncoming flow, a substantial vertical velocity is imparted to the fluid by the body. Furthermore, a counter-rotating vortex pair is observed in the wake and it induces a significant vertical velocity between the vortices. Figure 11 compares the mean vertical velocity ( $U_z$ ) between a moderate value of pitch angle ( $\alpha = 10^\circ$ , R90F $\infty$ ) and zero pitch ( $\alpha = 0^\circ$ , R90F $\infty$ A0). In the pitched case R90F $\infty$ (figure 11 a–c), a region of negative $U_z$ is observed in the centre plane ( $y=0$ ) at all streamwise locations. These regions of negative $U_z$ are flanked by two weaker positive $U_z$ regions on either side at $x/D=10$ and $20$ . This is consistent with a $U_z$ field generated by a pair of counter-rotating vortices, with a positive vortex on the right and a negative vortex on the left, see figure 12. Furthermore, when compared with the mean defect velocity $U_d$ (figure 9 a,c), $U_z$ exhibits a similar magnitude to $U_d$ at all streamwise locations and, therefore, becomes significant to the evolution of wake turbulence statistics.

Figure 11. Contours of mean vertical velocity ( $U_z$ ) at three streamwise locations $x/D=3$ , $x/D=10$ and $x/D=20$ . Panels (a)–(c) corresponds to the $\alpha =10^\circ$ case R90F $\infty$ and panels (d)–( f) to the $\alpha =0^\circ$ case R90F $\infty$ A0. Radial extent ( $r/D$ ) of the domain is shown on each panel.

The unpitched R90F $\infty$ A0 case in figure 11(d–f) exhibits a significantly different $U_z$ field compared to the moderate pitched case. When the flow leaves the body at $x/D=3$ , the magnitude of $U_z$ is approximately 10 times smaller than that in the pitched case. Due to the absence of a counter-rotating vortex pair, no region of negative $U_z$ is observed in the centre-plane. Furthermore, as the flow progresses downstream, $U_z$ remains 10 times smaller compared with the pitched case.

4.3. Buoyancy effects on streamwise vorticity

The shedding of a streamwise vortex pair from a body and its vertical drift has been well studied for an homogeneous fluid in the literature. Recently, Ohh & Spedding (Reference Ohh and Spedding2024) and NJOS25 studied the vertical drift of the vortex wake for an inclined 6 : 1 spheroid at $\alpha =10^\circ$ . Unlike the unequal-strength vortices formed by the asymmetric separation reported by NJOS25, the vortex pair at the present higher ${\textit{Re}} = 9\times 10^4$ (figure 12) is symmetric. Since the downward drift of the vortex pair and its modelling was discussed by NJOS25, it is not further considered here.

Figure 12. Contours of mean streamwise vorticity ( $\langle \omega _x\rangle$ ) for R90F $\infty$ at three streamwise locations: (a) $x/D=5$ ; (b) $x/D=10$ and (c) $x/D=20$ . Radial extent ( $r/D$ ) of domain is shown on each panel.

Figure 13. Contours of (a–d) mean streamwise vorticity ( $\langle \omega _x\rangle$ ) and (e–h) baroclinic torque ( $\omega _{BT}$ ) for R90F6 at four streamwise locations $x/D=5$ , $x/D=10$ , $x/D=20$ and $x/D=30$ . Isopycnals are overlaid on the $\omega _{BT}$ contours (black). Radial extent ( $r/D$ ) of the domain is shown on each panel.

Buoyancy is found to substantially affect the evolution of streamwise vorticity even for the relatively high $\textit{Fr} = 6$ . Figure 13(a–d) shows the mean streamwise vorticity ( $\langle \omega _x\rangle$ ) for R90F6 at four locations. Additionally, figure 13(e–h) shows the corresponding streamwise component of the mean baroclinic torque ( $-(1/{Fr})^2\partial \langle \rho \rangle /\partial y$ ) denoted succinctly by $\omega _{BT}$ . At $x/D=5$ , the primary vortex pair looks very similar to R90F $\infty$ (figure 12 a), with positive on the right and negative on the left. The induced velocity displaces the vortex pair vertically downwards. However, $\omega _{BT}$ has patches opposing the existing vorticity field with negative towards the right and positive towards the left. The opposing torque has two notable effects. First, the strength of the primary vortex pair is reduced so that it is substantially weaker at $x/D =10$ and is decimated by $x/D = 20$ . Second, its cumulative action generates two new secondary vortex pairs, one above and another below the centreline. For instance, compare $\omega _{BT}$ ( $x/D = 5)$ in the top half of figure 13(e) with $\langle \omega _x\rangle$ ( $x/D = 10$ ) in the top half of figure 13(b) and, similarly, the bottom halves of $\omega _{BT}$ ( $x/D = 10)$ and $\langle \omega _x\rangle$ ( $x/D = 20$ ).

Figure 14. Schematic for estimating the streamwise component of baroclinic torque ( $\omega _{BT}=-(1/{Fr})^2\partial \langle \rho \rangle /\partial y$ ).

In the near wake, the initial value of $\omega _{BT}$ not only depends on the stratification, but also on the geometry of the wake generator. We estimate the initial $\omega _{BT}$ as follows. Superscript $^*$ will distinguish dimensional quantities (when used) from their non-dimensional counterparts. Consider two fluid particles upstream of the spheroid (figure 14), where green has a vertical offset of $\delta z^*$ and a spanwise offset of $\delta y^*$ , and assume that the particles take the shown simplified paths. Initially, blue is vertically offset from green by $-\delta z^* = -6D \sin (\alpha )$ and is heavier than green by

(4.1) \begin{equation} \delta \rho ^* = -6D \sin (\alpha )C, \end{equation}

where $C$ is the dimensional background density gradient in the vertical. Reverting back to non-dimensional variables with $D$ for length and $CD$ for the density variation, it follows that

(4.2) \begin{equation} \delta \rho = -6 \sin (\alpha ){\textrm {sgn}}(C), \end{equation}

where $\textrm {sgn}(C) = -1$ . After travelling their simplified pathlines, both particles eventually reach the same vertical level, while maintaining the initial spanwise offset ( $\delta y^* = D/2$ ) of green relative to blue or $\delta y = 1/2$ , leading to a spanwise non-dimensional density gradient,

(4.3) \begin{equation} \frac {\partial \langle \rho \rangle }{\partial y} \sim \frac {-\delta \rho }{\delta y} = 12 \sin (\alpha ) \textrm {sgn}(C)\, . \end{equation}

The baroclinic torque $\omega _{BT}$ originating from the motion of the blue/green particle pair is

(4.4) \begin{equation} \omega _{BT} = \frac {-1}{{Fr}^2}\frac {\partial \langle \rho \rangle }{\partial y} = -\frac {12 \sin (\alpha )}{{Fr}^2} \textrm {sgn}(C) . \end{equation}

For the case R90F6, where $\alpha =10^\circ$ and ${\textit{Fr}}=6$ , the above-mentioned estimate leads to $\omega _{BT} = 0.06$ . Notwithstanding the assumption of simplified pathlines, we arrive at a reasonable estimate since the simulation result (figure 13 e) has $\omega _{BT} = O(0.1)$ . The sign is also correctly estimated for the $y\lt 0$ half-plane considered here. Since, $\omega _{BT}$ flips sign if a similar argument is followed for a blue particle that originates at $y^* = D/2$ instead of $y^* = - D/2$ , the final result for the initial value of non-dimensional baroclinic torque is written as

(4.5) \begin{equation} \omega _{BT} = \pm \frac {12 \sin (\alpha )}{{Fr}^2} . \end{equation}

By $x/D=10$ , a secondary vortex pair is formed directly above the primary vortex pair, but has the opposite sense of rotation (figure 13 b). Correspondingly, its motion reverses to vertically upward. By $x/D =20$ , the primary vortex has been replaced by two secondary vortex pairs, one above and the other below the centreline. Each secondary vortex pair drifts upward, $\omega _{BT}$ reverses sign too and acts to oppose the secondary vorticity. Thus, the vortex pair experiences a vertical oscillation.

It is worth noting that, at the lower ${\textit{Re}}$ of NJOS25 where the $\textit{Fr} = 6$ case exhibited asymmetry in the vortex pair formed after flow separation, buoyancy effects on the unequal-strength vortices resulted in intertwined patches of opposite-signed baroclinic torque that weakened the amplitude and coherence of the generated secondary vorticity. Consequently, the vertical oscillation of the wake found by Ohh & Spedding (Reference Ohh and Spedding2024) and NJOS25 is stronger in the present R90F6 case where the secondary vortex pair is more coherent.

5.1. Streamwise evolution of TKE

Figure 15(a) depicts the streamwise evolution of area-integrated TKE ( $\{E^T_K\}=\int _A E^T_K \,{\rm d}A$ ) for the R90 series. The area integral is computed over a circular region of radius $r=4D$ . All $\textit{Fr}$ cases exhibit a peak in TKE immediately after leaving the body at $x/D\approx 3$ with comparable values. Beyond this location, a monotonic decay is observed. This location also corresponds to the breakup of wake into small-scale structures. The unstratified R90F $\infty$ wake displays a power-law decay of $\{E^T_K\}\sim x^{-2/3}$ , whereas the weakly stratified case R90F6 decays faster, following $\{E^T_K\}\sim x^{-1.1}$ . It is important to note that these fits are empirical, do not represent self-similar power laws and are employed here for differentiating among the various cases. The strongly stratified case R90F1.9 exhibits the lowest TKE of all $\textit{Fr}$ cases at ${\textit{Re}}=9\times 10^4$ throughout the evolution, with the decay slowing down towards the end of the domain. For this case, buoyancy-induced oscillations can also be observed in $\{E^T_K\}$ after $x/D\approx 10$ .

Figure 15. Area-integrated TKE evolution in the streamwise direction at (a) ${\textit{Re}}=9\times 10^4$ and (b) ${\textit{Re}}=3\times 10^4$ for $\textit{Fr} = \infty , 6$ and $1.9$ . All cases are for $\alpha =10^\circ$ . Dotted lines denote empirical curve fit.

Area-integrated TKE evolution for the R30 wakes is depicted in figure 15(b) for all $\textit{Fr}$ cases. R30F $\infty$ and R30F6 wakes are indistinguishable until approximately $x/D\approx 7$ or $Nt\approx 1.1$ , consistent with the similarities observed in mean defect velocity (figure 9 b, f). After $Nt \approx 1.1$ , buoyancy effects in the ${\textit{Fr}}=6$ wake become progressively important as will be elaborated within the framework of the TKE transport equation in the following subsections. The maxima of TKE for both cases occur at $x/D\approx 5$ and is slightly delayed when compared with ${\textit{Re}}=9\times 10^4$ . These peak values correspond to locations of transition to turbulence when the flow, after laminar boundary-layer separation, breaks up into small-scale structures (see figure 4 of NJOS25). At low ${\textit{Re}}$ , the relatively stronger viscous forces delay transition and, correspondingly, the peak in TKE. Notably, both $\textit{Fr}$ cases exhibit faster power-law decay of $\{E^T_K\}$ compared with their higher- ${\textit{Re}}$ counterparts. The evolution of the unstratified R30F $\infty$ wake is $\{E^T_K\}\sim x^{-1}$ instead of $x^{-2/3}$ at ${\textit{Re}} = 9 \times 10^4$ and, for the ${\textit{Fr}} = 6$ cases, $\{E^T_K\}\sim x^{-1.6}$ for R30F6 instead of $x^{-1.1}$ for R90F6.

Unlike all other cases in figure 15, TKE for R30F1.9 exhibits a slow and gradual increase after flow separation. For this case, TKE peaks quite far downstream from the body at $x/D\approx 12$ . Subsequently, TKE levels remain relatively low as the flow progresses downstream.

5.2. TKE budget terms

In the preceding subsection, it was noted that wakes at the two ${\textit{Re}}$ levels exhibit distinct TKE decay rates. Moreover, these decay rates are significantly influenced by the background stratification. To further understand these ${\textit{Re}}$ and $\textit{Fr}$ trends, we quantify the following TKE transport equation:

(5.1) \begin{equation} U_i \frac {\partial E^T_K}{\partial x_i} + \frac {\partial T_i}{\partial x_i} = P - \varepsilon + B , \end{equation}

where $T_i$ is the turbulent transport, $P$ is turbulent production, $\varepsilon$ is turbulent dissipation rate and $B$ is turbulent buoyancy flux defined as

(5.2) \begin{equation} P = -\big\langle u^{\prime}_{i}u^{\prime}_j\big\rangle \frac {\partial U_i}{\partial x_j}, \qquad \varepsilon = 2\nu \big\langle s^{\prime}_{\textit{ij}}s^{\prime}_{\textit{ij}}\big\rangle - \big\langle \tau ^{{\prime}s}_{\textit{ij}}s^{\prime}_{\textit{ij}} \big\rangle , \qquad B = -\frac {g}{\rho _0}\big\langle \rho ^{\prime}u^{\prime}_z \big\rangle . \end{equation}

Here, $s_{\textit{ij}}=(\partial _ju_i+\partial _iu_j)/2$ is the strain rate tensor and $\tau ^s_{\textit{ij}}=-2\nu _ss_{\textit{ij}}$ is the subgrid stress tensor. The subgrid term in the TKE transport equation is found to be relatively small in the wake at both ${\textit{Re}}$ – a consequence of the good grid resolution employed here.

Production $P$ acts as a source in the transport equation, extracting energy from the MKE reservoir and generating an equivalent amount of TKE. Dissipation $\varepsilon$ serves as a sink of TKE. Buoyancy $B$ facilitates energy exchange between TKE and TPE (turbulent potential energy), with a positive buoyancy flux converting TPE into TKE and vice versa. Unlike the aforementioned source and sink terms, turbulent transport ( $T_i$ ) is responsible for spatial redistribution of TKE and is given by

(5.3) \begin{equation} T_i = \frac {1}{2}\big\langle u^{\prime}_{i}u^{\prime}_{j}u^{\prime}_{j}\big\rangle + \big\langle u^{\prime}_{i}p^{\prime}\big\rangle - 2\nu \big\langle u^{\prime}_{j}s^{\prime}_{\textit{ij}}\big\rangle - \big\langle u^{\prime}_{j}\tau ^{{\prime}s}_{\textit{ij}}\big\rangle . \end{equation}

In the wake, turbulent transport redistributes energy, primarily in the $y{-}z$ plane. Consequently, the contribution of the area integral of $\partial T_i/ \partial x_i$ to the area-integrated TKE balance is negligible.

Figure 16. Area-integrated turbulent (a,c) production and (b,d) dissipation evolution in streamwise direction at (a,c) ${\textit{Re}}=9\times 10^4$ and (b,d) ${\textit{Re}}=3\times 10^4$ for $\textit{Fr} = \infty , 6$ and $1.9$ . All cases are for $\alpha =10^\circ$ . Dotted lines denote empirical curve fit.

The streamwise evolution of area-integrated production, $\{P\}$ , is depicted in figure 16(a) for the R90 series. The  $\{P\}$ has already reached its peak for all $\textit{Fr}$ cases by $x/D=3$ , followed by a consistent decline. Stratification generally tends to decrease $\{P\}$ . For example, $\{P\}\sim x^{-1.96}$ in the R90F $\infty$ wake, while the R90F6 wake exhibits $\{P\}\sim x^{-2.65}$ , which is a faster decay. For R90F1.9, production initially exhibits a rapid decay compared with other cases until $x/D\approx 7$ , beyond which, the decay rate slows down and $\{P\}$ becomes comparable to the other cases.

Area-integrated dissipation, $\{\varepsilon \}$ , in the R90 series also shows a decay after reaching its peak at $x/D\approx 3$ , as illustrated in figure 16(c). For the R90F $\infty$ wake, $\{\varepsilon \}\sim x^{-1.85}$ , with the decay rates increasing as the stratification level increases, indicating a damping effect due to buoyancy.

The R30 cases are shown in figures 16(b) and 16(d). Buoyancy affects the evolution of $\{P\}$ and $\{\varepsilon \}$ in R30F6 relative to R30F $\infty$ similar to its influence for their R90 equivalents. For example, $\{P\}$ for both values of $\textit{Fr}$ experiences a peak at $x/D\approx 4$ and subsequently decays with rates similar to the corresponding R90 cases. Moving to ${\textit{Re}}$ , its effect on the wake TKE is weak if the $\textit{Fr}$ value is fixed at $\infty$ or $6$ . Although the initial values of $\{P\}$ and $\{\varepsilon \}$ are smaller and their initial evolution is different for the ${\textit{Re}}=3\times 10^4$ wakes, their eventual decay rates are similar to their higher- ${\textit{Re}}$ counterparts.

The R30F1.9 wake stands out from all the other cases discussed in this section due to the absence of small-scale structures, resulting in a relatively small level of velocity fluctuations in the near wake. Both $\{P\}$ and $\{\varepsilon \}$ exhibit the lowest initial values among all cases and, after an initial increase, remain nearly constant until $x/D\approx 12$ or $Nt\approx 6$ . This region also corresponds to the gradual increase in TKE, as depicted in figure 15(b). Beyond this point, both quantities undergo a decay, with production values being comparable and dissipation slightly higher than that of the R90F1.9 wake.

Figure 17. Ratio of area-integrated production and dissipation with respect to streamwise coordinate at (a) ${\textit{Re}}=9\times 10^4$ and (b) ${\textit{Re}}=3\times 10^4$ . All cases are for $\alpha =10^\circ$ .

The ratio $\{P\}/\{\varepsilon \}$ (figure 17) is a direct measure of the relative significance of production and dissipation in the TKE balance. When $\{P\}/\{\varepsilon \}\lt 1$ , dissipation consumes more energy than production contributes, leading to a decrease in TKE. For all cases at both Reynolds numbers, $\{P\}/\{\varepsilon \}$ exceeds 1, while the flow leaves the body at $x/D=3$ , resulting in a streamwise increase of TKE. However, this trend is short-lived, as $\{P\}/\{\varepsilon \}$ soon dips below 1. This occurrence initiates the decay of TKE, which was previously seen in figure 15. For instance, all cases at ${\textit{Re}}=9\times 10^4$ exhibit $\{P\}/\{\varepsilon \}\lt 1$ from $x/D\approx 4.5$ onwards, accompanied by the onset of TKE decay at the same location. Another example is the R30F1.9 wake, which exhibits an extended region of $\{P\}/\{\varepsilon \}\gt 1$ until $x/D\approx 12$ , corresponding to the gradual TKE increase before the decay commences.

The relationship between the power laws of TKE and dissipation in some of the simulated wakes can be understood as follows. When the TKE budget is predominantly dominated by dissipation, the area-integrated TKE transport equation simplifies to

(5.4) \begin{equation} U_x \frac {\partial \{E^T_K\}}{\partial x} = - \{\varepsilon \} . \end{equation}

Assuming a power law for TKE as $\{E^T_K\}\sim x^{-p}$ and assuming $U_x\sim U_\infty$ remains relatively unchanged, (5.4) leads to $\{\varepsilon \} \sim x^{-p-1}$ . The above-mentioned relationship between TKE and dissipation power-law exponents is typically observed in this study when $\{P\}/\{\varepsilon \}$ falls below approximately 0.5. This correlation can be seen for the wake R90F6, where $\{E^T_K\}\sim x^{-1.1}, \{\varepsilon \}\sim x^{-2.1}$ and $p\approx -1.1$ . Similarly, it can be observed for R30F $\infty$ , where $\{E^T_K\}\sim x^{-1}, \{\varepsilon \}\sim x^{-1.85}$ and $p\approx -1$ , and also R30F6, where $\{E^T_K\}\sim x^{-1.6}$ , $\{\varepsilon \}\sim x^{-2.47}$ and $p\approx -1.6$ . However, $\{P\}/ \{\varepsilon \}$ does not become as small as 0.5 in the R90F $\infty$ wake and the dissipation does not follow the $x^{-p-1}$ scaling.

Figure 18. (a,b) Area-integrated turbulent buoyancy flux $\{B\}$ and (c,d) ratio of area-integrated buoyancy and dissipation for ${\textit{Re}}=9\times 10^4$ and ${\textit{Re}}=3\times 10^4$ , respectively. All cases are for $\alpha =10^\circ$ . Note that the $\{ B \}$ -scale in panel (b) is an order of magnitude smaller than that in panel (a).

We now move to a discussion of the direct contribution of the buoyancy flux to the TKE balance in the stratified wakes. The streamwise evolution of the area-integrated buoyancy flux, $\{B\}$ , is depicted in figure 18(a) for ${\textit{Re}}=9\times 10^4$ . It is evident that the R90F1.9 wake exhibits a stronger (larger amplitude) buoyancy flux compared with the R90F6 wake. For a given vertical displacement of the turbulent motions ( $\Delta z_t$ ) in a background with buoyancy frequency $N^2$ , the resulting buoyancy is $ B \sim g/\rho _0 \, ( \partial \rho _b/\partial z) \Delta z_t = -N^2 \Delta z_t$ . Since the initial downward $\Delta z_t$ at $x/D =3$ , just after flow separation, is set by the inclination of the spheroid and the Reynolds number, the difference among cases is small as is suggested by comparing among the rows of figure 2. Therefore, the higher $N^2$ case R90F1.9 has a higher buoyancy flux.

Moreover, $\{B\}$ for both wakes exhibits significant modulation by lee waves and oscillates with a wavelength equal to the wavelength ( $\lambda \approx 2\pi {Fr}$ ) that has been previously reported for body generated lee waves. Figure 18(c) shows the evolution of $\{B\}/\{\varepsilon \}$ for ${\textit{Re}}=9\times 10^4$ . For R90F6, the contribution of buoyancy is relatively small to the TKE budget until $x/D\approx 20$ or $Nt\approx 3$ , beyond which it contributes up to half of the dissipation value. Whereas, for R90F1.9, buoyancy becomes important close to the body, starting at $x/D\approx 6$ or $Nt\approx 3$ ( $Nt$ same as for R90F6). Along with high buoyancy induced fluctuations, buoyancy occasionally becomes comparable to dissipation, converting a significant amount of TKE into TPE and vice versa. Therefore, $B$ directly affects the R90F1.9 wake for most of the domain.

Turning to the lower- ${\textit{Re}}$ series (figure 18 b), the R30F6 wake displays a comparable level of $\{B\}$ as its higher ${\textit{Re}}$ counterpart R90F6 (note the smaller $\{B\}$ scale in panel a relative to panel b). Interestingly, the R30F1.9 wake exhibits a smaller amplitude of $\{B\}$ than R30F6, despite having a higher $N^2$ . The R30F1.9 near wake is quasi-laminar and lacking small-scale structures (see figure 4c of NJOS25), it experiences smaller density and velocity fluctuations (analogous to its low level of TKE and $\{B\}$ ) until $x/D\approx 12$ . Consequently, the buoyancy flux is small in this case. Figure 18(d) shows the evolution of the ratio $\{B\}/\{\varepsilon \}$ in the R30 series. For both stratified cases, $\{B\}/\{\varepsilon \}$ is smaller than for their higher ${\textit{Re}}$ counterparts and shows similar buoyancy-induced oscillations. The R30F1.9 wake shows slightly higher $\{B\}/\{\varepsilon \}$ towards the end of domain, but is less than $0.5$ for the most part.

The contrast in the TKE balance between ${\textit{Fr}}=\infty$ and $6$ can be summed up at both ${\textit{Re}}$ as follows. Production for ${\textit{Fr}}=6$ is lower than ${\textit{Fr}}=\infty$ , bringing in less TKE into the stratified case budget. Dissipation for ${\textit{Fr}}=6$ is also lower than ${\textit{Fr}}=\infty$ , taking less TKE out of the stratified budget. Buoyancy flux for ${\textit{Fr}}=6$ is not sufficient to directly influence the stratified budget. Also, just after the flow separates from the spheroid, TKE is similar for both wakes. Taking the aforementioned features into account, it follows that the reduced production at ${\textit{Fr}}=6$ is the key reason for faster TKE decay compared with ${\textit{Fr}}=\infty$ . Although the buoyancy flux is small, it plays a significant role in indirectly affecting the turbulence statistics, as will be discussed in the next subsection.

5.3. Components of turbulent production

The production term in (5.2) can be further expanded into

(5.5) \begin{equation} P = -\big\langle u^{\prime}_{x}u^{\prime}_{j}\big\rangle \frac {\partial U_x}{\partial x_j}-\big\langle u^{\prime}_{y}u^{\prime}_{j}\big\rangle \frac {\partial U_y}{\partial x_j}-\big\langle u^{\prime}_{z}u^{\prime}_{j}\big\rangle \frac {\partial U_z}{\partial x_j}. \end{equation}

Figure 19. (a) Area-integrated values of components $\{P_{xy}\}$ , $\{P_{\textit{xz}}\}$ and total production $\{P\}$ , and (b) the fractional contributions of the dominant components, $\{P_{\textit{ij}}\}/\{P\}$ , for the R90F $\infty$ A0 wake.

The previous section established that the area-integrated production, $\{P\}$ , exhibits a faster decay rate at ${\textit{Fr}}=6$ compared with ${\textit{Fr}}=\infty$ at both Reynolds numbers and leads to lower TKE. To delve deeper into the role of buoyancy in this behaviour, the contribution of all dominant production component terms, denoted as $P_{\textit{ij}}=-\langle u^{\prime}_{i}u^{\prime}_{j}\rangle \partial _jU_i$ , are quantified. A spheroid wake that develops in the $x$ direction is a laterally thin shear flow implying that $U_y, U_z \ll U_x$ and $|\partial U_x/\partial x|\ll |\partial U_x/\partial y|,|\partial U_x/\partial z|$ . Consequently, $P$ is dominated by transverse terms, specifically $P_{xy}=-\langle u^{\prime}_{x}u^{\prime}_{y}\rangle \partial _yU_x$ and $P_{\textit{xz}}=-\langle u^{\prime}_{x}u^{\prime}_{z}\rangle \partial _zU_x$ . In the present simulations of an unstratified spheroid wake at zero angle of attack (R90F $\infty$ A0), it is confirmed that $P_{xy}$ and $P_{\textit{xz}}$ indeed constitute the dominant components of production, as illustrated in figure 19(a). Similar to $\{P\}$ , both $P_{xy}$ and $P_{\textit{xz}}$ exhibit a monotonic decline after reaching a peak at $x/D\approx 5$ . Each of $P_{xy}$ and $P_{\textit{xz}}$ contributes approximately $50\,\%$ of the total $P$ . Note that the total contribution of these two terms near body at $x/D=4$ is more than 100 % of $P$ . This is because of the significant negative contribution by $P_{xx}=-\langle u^{\prime}_{x}u^{\prime}_{x}\rangle \partial _xU_x$ in the near wake caused by rapid acceleration of the mean flow near the separation region.

The situation is different for a pitched $6:1$ spheroid, despite a moderate $\alpha =10^\circ$ . Owing to the pitch and flow separation into a counter-rotating vortex pair, significant vertically downward velocity is imparted to the fluid, as was illustrated in § 4 (figure 11). Therefore, in addition to $P_{xy}$ and $P_{\textit{xz}}$ , terms related to mean vertical velocity, such as $P_{zy}=-\langle u^{\prime}_{z}u^{\prime}_{y}\rangle \partial _yU_z$ and $P_{zz}=-\langle u^{\prime}_{z}u^{\prime}_{z}\rangle \partial _zU_z$ , also contribute crucially to the total production. Since the relative behaviour of production components for ${\textit{Fr}}=\infty$ and ${\textit{Fr}}=6$ is qualitatively similar for both ${\textit{Re}}$ cases, except during an initial transition period, the following discussion is limited to the higher ${\textit{Re}}=9\times 10^4$ cases for brevity.

Figure 20. Area-integrated production components $\{P_{xy}\}$ , $\{P_{\textit{xz}}\}$ , $\{P_{zy}\}$ , $\{P_{zz}\}$ and total production $\{P\}$ for (a) R90F $\infty$ and (b) R90F6 wakes. Ratio of production components and total production, $\{P_{\textit{ij}}\}/\{P\}$ , for (c) R90F $\infty$ and (d) R90F6.

Figure 20(a) illustrates the evolution of the dominant production components alongside the total production for R90F $\infty$ . Similar to the total production, all the production components also reach a peak near $x/D\approx 3$ followed by a monotonic decline. Notably, $P_{zy}$ and $P_{zz}$ consistently demonstrate higher levels compared with $P_{xy}$ and $P_{\textit{xz}}$ for $x/D\gtrsim 5$ . As illustrated by figure 20(c), beyond $x/D\approx 5$ , $P_{zy}$ contributes approximately 35 % and $P_{zz}$ contributes approximately 30 % to the total $P$ . In contrast, $P_{xy}$ and $P_{\textit{xz}}$ both contribute only approximately 20 % each. These results highlight the importance of spanwise and vertical gradients of $U_z$ to TKE production. This is in stark contrast to the zero pitch case R90F $\infty$ A0 discussed previously (figure 19) where only the gradients of streamwise velocity ( $U_x$ ) are important.

The $P_{xy}$ for R90F6 exhibit nearly the same levels as those for R90F $\infty$ (figure 20 b), while $P_{\textit{xz}}$ shows a slight reduction. However, $P_{zy}$ and $P_{zz}$ show a much faster decay after $x/D\approx 4$ compared with their unstratified wake counterparts. Eventually, for the stratified wake, $P_{xy}$ and $P_{\textit{xz}}$ components become the dominant ones after $x/D\approx 7$ and remain so until the end of the domain. The ratio $\{P_{\textit{ij}}\}/\{P\}$ (figure 20 d) for both $P_{xy}$ and $P_{\textit{xz}}$ consistently increases for $x/D\gtrsim 5$ , with $P_{xy}$ contributing approximately 55 % and $P_{\textit{xz}}$ contributing approximately 40 % of the total $P$ towards the end of the domain. In contrast, $P_{zy}$ and $P_{zz}$ show reduced contributions, each contributing only approximately 15 % of the total $P$ towards the end of the domain. Since $P_{xy}$ is similar at both ${\textit{Fr}}=6$ and ${\textit{Fr}}=\infty$ , it follows that the reduction in $P_{\textit{xz}}$ , $P_{zy}$ and $P_{zz}$ contributions is responsible for the faster decay of the total production in the stratified wake.

Similar to R90F $\infty$ A0, the sum of all dominant production terms in both R90F $\infty$ and R90F6 wakes is slightly greater than 100 % of the total $P$ . This is due to the small negative contribution of certain production terms like $P_{xx}, P_{yy}$ and $P_{yz}$ .

5.4. Indirect effects of buoyancy flux on the turbulent statistics

Quantification of the Reynolds stresses ( $R_{\textit{ij}}=\langle u^{\prime}_{i}u^{\prime}_{j} \rangle$ ) involved in the dominant $P_{\textit{ij}}$ terms leads to further insights into the differential decay of production components under stratification. Figure 21 illustrates the streamwise evolution of the area-integrated absolute Reynolds stress components: $\{|R_{xy}|\}, \{|R_{yz}|\}, \{|R_{\textit{xz}}|\}, \{|R_{zz}|\}$ for R90F $\infty$ and R90F $6$ wakes. Since $R_{\textit{ij}}$ can be either positive or negative, we computed the area integral of its magnitude to capture the local impact of the Reynolds stress. It is evident that all $\{|R_{\textit{ij}}|\}$ components in both R90F $\infty$ and R90F6 cases commence with remarkably similar values after the flow separates from the body at approximately $x/D\approx 3$ . The $\{|R_{xy}|\}$ for both cases exhibit a close resemblance until $x/D\approx 20$ , beyond which a slightly faster decay is observed for the stratified case (figure 21 a). Since $\{|R_{xy}|\}$ undergoes minimal changes with stratification, it can be inferred that buoyancy has a negligible influence on altering $P_{xy}$ .

Figure 21. Area-integrated absolute Reynolds stress (a) $\{|R_{xy}|\}$ , (b) $\{|R_{yz}|\}$ , (c) $\{|R_{\textit{xz}}|\}$ and (d) $\{|R_{zz}|\}$ for R90F $\infty$ and R90F6.

The other Reynolds stresses, $\{|R_{yz}|\}, \{|R_{\textit{xz}}|\}$ and $\{|R_{zz}|\}$ , experience a faster decay when subjected to stratification of ${\textit{Fr}}=6$ compared with the unstratified case. For instance, $\{|R_{\textit{xz}}|\}$ for R90F6 begins to decay faster than R90F $\infty$ at $x/D\approx 4$ , with a significant drop observed after $x/D\approx 20$ (figure 21 c). This aligns with the observation of $P_{\textit{xz}} \lt P_{xy}$ in the ${\textit{Fr}}=6$ wake throughout the domain, unlike the ${\textit{Fr}}=\infty$ wake where $P_{\textit{xz}} \approx P_{xy}$ was observed (figure 20). Similarly, $\{|R_{yz}|\}$ and $\{|R_{zz}|\}$ also decay much faster for ${\textit{Fr}}=6$ than ${\textit{Fr}}=\infty$ for $x\gtrsim 4$ (figure 21 b, d). These results lead to a reduced contribution of $P_{yz}$ and $P_{zz}$ to the total production, as evident in figure 20.

One crucial point to note is that only the Reynolds stresses with $u^{\prime}_z$ , such as $R_{iz}=\langle u^{\prime}_i u^{\prime}_z \rangle$ , experience a reduction in the stratified case. In contrast, stratification appears to have a negligible impact on other Reynolds stresses. This phenomenon can be attributed to the role of buoyancy in the transport equations for $R_{\textit{ij}}$ . The $R_{\textit{ij}}$ transport equation for an incompressible flow without buoyancy, e.g. Pope (Reference Pope2000), needs to be supplemented by the following buoyancy term on its right-hand side,

(5.6) \begin{equation} B_{\textit{ij}} = -\frac {g}{\rho _0}\big(\big\langle \rho ^{\prime}u^{\prime}_{i}\big\rangle \delta _{jz} + \big\langle \rho ^{\prime}u^{\prime}_{j}\big\rangle \delta _{iz}\big). \end{equation}

The $B_{\textit{ij}}$ term specifically affects the Reynolds stress terms with $u^{\prime}_z$ , leading to a reduction in $\{|R_{\textit{xz}}|\}, \{|R_{yz}|\} \text{ and } \{|R_{zz}|\}$ . Consequently, this reduction also impacts the associated production components $P_{\textit{xz}}, P_{yz} \text{ and }P_{zz}$ . Note that $B_{zz} = 2B$ , where $B$ is the buoyancy term in the TKE equation.

It was shown in figure 18(a) that $\{B_{zz}\}$ ( $=2\{B\}$ ) remains negative throughout the evolution and, since it appears on the right-hand side of the $R_{zz}$ transport equation, would tend to decrease $R_{zz}$ . Moreover, $\{|R_{zz}|\}$ begins to decay rapidly at the location where $\{B_{zz}\}$ peaks, approximately $x/D\approx 4$ . This marks the threshold beyond which buoyancy starts influencing wake turbulence. The $B_{\textit{xz}}$ and $B_{yz}$ terms are not presented here for brevity.

To summarise, the results of this section show that, although buoyancy is not a significant sink of TKE, it plays an indirect role in significantly altering it. Buoyancy reduces turbulent production by suppressing its specific components that depend on the three Reynolds stress components that involve vertical velocity fluctuations ( $u^{\prime}_z$ ).