1. Introduction
1.1. Internal waves and their impact
Relative motion between a submerged body and its stratified environment gives rise to different types of internal gravity waves. Observations of these waves in the ocean and the atmosphere prompted the first studies of these waves in the geophysical context. Flow over mountains gives rise to lee waves that can lead to changes in properties of air, sometimes visualised as wave clouds. In the ocean, substantial energy is transferred from impinging mean currents and tides at bottom topography to internal waves, e.g. Garrett & Kunze (Reference Garrett and Kunze2007), which, through nonlinear processes break down to turbulence, e.g. Sarkar & Scotti (Reference Sarkar and Scotti2017), and the associated turbulent mixing provides an important control on ocean stratification and currents, e.g. Wunsch & Ferrari (Reference Wunsch and Ferrari2004). The motivating application – bodies moving in a stratified environment – of the present work is different. Here, the body generates lee waves, which are steady in the reference frame of the body, and the unsteady motions in the wake also generate waves (unsteady in the frame of the body). Characterisation of the space–time structure and energetics of these waves is of interest as is quantification of their relationship with the wake and the body generator.
The internal or body Froude number, 
${{Fr}} = U / ND$, is the important overall parameter that governs internal gravity waves and buoyancy effects on wake dynamics. Here, the density-stratified background is characterised by the buoyancy frequency 
$N$ where 
$N^2 = - (g/\rho _0) \partial \rho _b /\partial z$, the relative velocity between the body and the ambient by 
$U$ and the size of the body by 
$D$. In geophysical flows, 
${{Fr}}$ is generally less or much less than 
$O(1)$, while in engineering applications (underwater submersibles, aerial vehicles, wind farms, marine turbines) 
${{Fr}} \geqslant O(1)$.
1.2. Lee wave studies in the laboratory and with simulations
Laboratory experiments of stratified flow past a model hill (Hunt & Snyder Reference Hunt and Snyder1980) and model ridges (Castro, Snyder & Marsh Reference Castro, Snyder and Marsh1983) have characterised lee waves through extensive flow visualisation. Chomaz, Bonneton & Hopfinger (Reference Chomaz, Bonneton and Hopfinger1993) experimentally studied and classified the change of near-wake structure of a sphere wake in uniform stratification as 
${{Fr}}$ is decreased from greater-than to less-than 
$O(1)$ values. In a companion study, Bonneton, Chomaz & Hopfinger (Reference Bonneton, Chomaz and Hopfinger1993) described wave characteristics on horizontal planes above the wake using isopyncal displacement inferred from fluorescent dye images. Non-uniform stratification was treated by Robey (Reference Robey1997) who, using experiments and numerics, found the wave field of a sphere moving at the base of a thermocline to be composed of discrete modes that give rise to horizontal-plane patterns reminiscent of Kelvin wakes on the air–sea boundary. Recently, Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) investigated the wave field of bodies of different shapes and determined how the wavelength and velocity amplitude of lee waves varied with 
${{Fr}}$. They also proposed a model for lee wave amplitude that involves modelling the drag as a point force in the Navier–Stokes equations.
Various numerical simulations have also described these lee waves and their effects on the flow at the body and in the wake. For a sphere at 
${{Re}} =200$, Hanazaki (Reference Hanazaki1988) found that lee waves suppress separation on decreasing 
${{Fr}}$ when 
${{Fr}} > 1$ but induce separation on decreasing 
${{Fr}}$ when 
${{Fr}} < 1$. Lee waves lead to oscillatory modulation of the wake width and velocity (Pal et al. Reference Pal, Sarkar, Posa and Balaras2017) and also, near the body, lead to a local minimum of the wake defect velocity at 
$Nt = {\rm \pi}$ corresponding to the half-wavelength of the lee wave (Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020). Here, 
$Nt = Nx/U = x/({{Fr}} D)$ converts distance behind the body to time in buoyancy units. Lee waves have a stronger effect on the wake of a slender body, i.e. large value of aspect ratio (
$AR$) given by 
$L/D$ where 
$L$ is the length of the body. When the half-wavelength of the lee wave matches 
$L$, a condition that is realised at a critical value of 
${{Fr}}_{c} = AR/{\rm \pi}$, the wave leads to maximal contraction of the separated flow at the trailing edge. Thus, in a large-eddy simulations (LES) study (Ortiz-Tarin, Chongsiripinyo & Sarkar Reference Ortiz-Tarin, Chongsiripinyo and Sarkar2019) of flow past a 
$4:1$ spheroid at 
${{Re}} = 10^4$, the wake which was turbulent at high 
${{Fr}}$ exhibited strong suppression of wake turbulence at 
${{Fr}} = 1$. The problem of a disk at 
${{Re}} = 5000$ in nonlinear stratification with minimum 
${{Fr}} =1$ was studied by Gola et al. (Reference Gola, Nidhan, Ortiz-Tarin, Pham and Sarkar2023). Linear theory (Voisin Reference Voisin1991, Reference Voisin2003, Reference Voisin2007) was used with an equivalent body comprising of the disk and its separation bubble and found to predict the wave amplitude of a single test case (
${{Fr}} =1$) with linear stratification but not the nonlinearly stratified cases. In the present work, lee waves in linear stratification will be characterised over a wide range of 
${{Fr}}$ and the applicability of linear theory will be evaluated.
1.3. Evolution of stratified turbulent wakes
Turbulent stratified wakes exhibit a stark contrast in the evolution of their deficit velocity and length scales relative to their unstratified counterpart. Although the present focus is on internal waves, previous work in this area is briefly summarised for completeness.
Several experimental studies (Lin & Pao Reference Lin and Pao1979; Lin, Boyer & Fernando Reference Lin, Boyer and Fernando1992; Spedding, Browand & Fincham Reference Spedding, Browand and Fincham1996; Spedding Reference Spedding1997; Bonnier & Eiff Reference Bonnier and Eiff2002) have reported significant buoyancy effects on the evolution of wake thickness and/or wake deficit in the flow behind a horizontally moving sphere. Spedding (Reference Spedding1997) found that the decay of the wake deficit velocity is characterised by three stages: the near wake regime; the non-equilibrium (NEQ) regime, where the wake decay rate is reduced as it adjusts to the surrounding stratification; and the quasi-two-dimensional regime, where pancake vortices form and the wake decay speeds up. The transition between stages occurred at nominally fixed values of 
$Nt$. In another laboratory study, Spedding (Reference Spedding2001) showed that stratification biases turbulence towards a faster decay of the vertical fluctuations as compared with the horizontal fluctuations.
Numerical investigations of the problem started with a temporal model (Gourlay et al. Reference Gourlay, Arendt, Fritts and Werne2001) where the wake is simulated in a frame moving with the body, a streamwise periodic domain is used and flow statistics evolve in time. The model does not include the wake generator and is instead initialised with a flow field that approximates experimentally observed mean velocity and turbulence intensities at some distance from the body. Since the boundary layer on the body and the small scales in the near wake do not need to be resolved, the temporal model has the advantage that the simulated flow can progress into the late wake with reasonable computational cost. The limitations are that the absence of the wake generator excludes lee waves and their influence on the wake as well as vortex shedding from the body. Also, temporal models are sensitive to initial flow conditions (Redford, Castro & Coleman Reference Redford, Castro and Coleman2012). The alternative approach of body-inclusive simulations will be adopted here.
The multistage decay of the wake has been studied numerically using temporal models (Gourlay et al. Reference Gourlay, Arendt, Fritts and Werne2001; Dommermuth et al. Reference Dommermuth, Rottman, Innis and Novikov2002; Brucker & Sarkar Reference Brucker and Sarkar2010; Diamessis, Spedding & Domaradzki Reference Diamessis, Spedding and Domaradzki2011; de Stadler & Sarkar Reference de Stadler and Sarkar2012; Zhou & Diamessis Reference Zhou and Diamessis2019). The anisotropy created by stratification suppressing vertical fluctuations leads to decreased turbulent production in the wake (Brucker & Sarkar Reference Brucker and Sarkar2010), leading to the longer lifespan of the stratified wake, a result that is supported by Redford, Lund & Coleman (Reference Redford, Lund and Coleman2015) in their direct numerical simulations (DNS) of a weakly stratified turbulent wake. The buoyancy-induced decrease in the correlation coefficient between vertical and horizontal fluctuations is a generic feature of stratified shear flows and was identified by Jacobitz, Sarkar & VanAtta (Reference Jacobitz, Sarkar and VanAtta1997) in homogenous shear flow and also by others in following studies of a variety of shear flows.
In recent times, body-inclusive simulations, e.g. Orr et al. (Reference Orr, Domaradzki, Spedding and Constantinescu2015), Pal et al. (Reference Pal, Sarkar, Posa and Balaras2017), Ortiz-Tarin et al. (Reference Ortiz-Tarin, Chongsiripinyo and Sarkar2019) and Chongsiripinyo & Sarkar (Reference Chongsiripinyo and Sarkar2020), have been conducted to increase the realism of simulations. Such an approach captures flow separation, vortex shedding from the body and also the lee waves leading to more accurate representation of the near and intermediate wake. Body inclusive simulations are akin to laboratory experiments by virtue of including the body, while also alleviating the aforementioned limitations of temporal-model simulations. Besides, body-inclusive simulations lead to statistically steady data so that numerical techniques like spectral proper orthogonal decomposition can be fruitfully used to extract wake structures that are statistically space–time coherent, e.g. Nidhan et al. (Reference Nidhan, Chongsiripinyo, Schmidt and Sarkar2020) and Nidhan, Schmidt & Sarkar (Reference Nidhan, Schmidt and Sarkar2022). Reaching beyond 
$x/D = O(100)$ in body-inclusive simulations is computationally too expensive but this issue can be overcome by a hybrid simulation technique that continues the flow to larger 
$x/D$ (equivalently 
$Nt$) in a follow-up simulation that uses a temporal model (Pasquetti Reference Pasquetti2011) or a spatial model (VanDine, Chongsiripinyo & Sarkar Reference VanDine, Chongsiripinyo and Sarkar2018) with a coarser grid that needs to resolve only the large length scales of the far wake.
1.4. Wake waves and their link to wake structures
Wake waves are internal waves forced by the unsteady fluctuations in the turbulent wake. The first report of stratified wake waves was by Gilreath & Brandt (Reference Gilreath and Brandt1985) who observed ‘short, random’ unsteady waves in their experiments with a streamlined body in towed and self-propelled modes. Besides uniform stratification, they also considered a pycnocline to find solitary internal wave propagation at low 
${{Fr}}$. Bonneton et al. (Reference Bonneton, Chomaz and Hopfinger1993) distinguished wake waves from lee waves and, by spectral analysis of density-probe measurements, linked these waves to coherent wake structures. Robey (Reference Robey1997) linked the wake waves to the size of large-scale structures in the wake. Towed sphere experiments by Brandt & Rottier (Reference Brandt and Rottier2015) showed that lee waves dominate in the 
${{Fr}} \lesssim 1$ regime while wake waves dominate in the 
${{Fr}} \gtrsim 1$ regime. They also provide a qualitative estimate of the total potential energy that was found to scale as 
${{Fr}}^{2}$. The recent work of Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) considered wake waves in addition to lee waves and quantified the dependence of their wavelength and velocity amplitude on 
${{Fr}}$.
Abdilghanie & Diamessis (Reference Abdilghanie and Diamessis2013) conducted a comprehensive study of the internal wave field of a stratified turbulent wake using temporal simulations for a wide range of 
${{Re}}$ and found that, at higher 
${{Re}}$, there is a broader range of wave propagation angles (
$40^{\circ }\unicode{x2013}55^{\circ }$), at least at early 
$Nt$, relative to lower values of 
${{Re}}$. Brucker & Sarkar (Reference Brucker and Sarkar2010) studied towed and self-propelled wakes using temporal simulations at 
${{Re}} = 5 \times 10^{4}$ and various 
${{Fr}}$. They concluded that internal wave flux dominated the turbulent dissipation in the wake kinetic energy budget for 
$20 < Nt < 75$. This was also supported by Rowe, Diamessis & Zhou (Reference Rowe, Diamessis and Zhou2020) who, in their temporal-model simulations over a broad range of 
${{Re}}$ and 
${{Fr}}$, found that internal wave radiation is an important sink for wake kinetic energy after 
$Nt = 10$. Indeed, the energetic importance of internal waves for turbulence in stratified shear flows is a more general result in view of such demonstrations for a shear layer (Pham, Sarkar & Brucker Reference Pham, Sarkar and Brucker2009; Watanabe et al. Reference Watanabe, Riley, Nagata, Onishi and Matsuda2018) and a boundary layer (Taylor & Sarkar Reference Taylor and Sarkar2007).
1.5. Objectives of present study
Body-inclusive simulations capture lee waves as well as wake waves and present an opportunity to characterise and contrast their properties. This motivates the present waves-focussed study of a disk wake in linear stratification at 
${{Re}} = 5000$.
Our results on scaling laws for 
${{Fr}}$ dependency of internal waves in turbulent flow past a disk will be compared with previous laboratory experiments and, for wake waves, also with body exclusive temporal model simulations. The skill of linear theory in predicting lee waves will be assessed. The internal wave field in cases from a recent study of disk wakes at an order of magnitude larger 
${{Re}} = 50\,000$ that focussed on wake decay laws will also be analysed.
The specific questions that we intend to answer as stratification is varied parametrically for a turbulent wake are as follows. (i) How does linear theory applied to an equivalent body, an approach which was shown to predict body-generated lee waves in a turbulent disk wake at a single 
${{Fr}} =1$, work over the range of 
${{Fr}}$ considered here? For example, how do linear theory results, 
$w \propto {{Fr}}^{-1}$ scaling of the amplitude and 
$\varLambda = 2{\rm \pi} {{Fr}} D$ scaling of wavelength, compare with the simulation results? (ii) How do internal wave structural properties (phase angle, wavelength, spatial distribution) change as a function of 
${{Fr}}$? Are there scaling laws for the 
${{Fr}}$ dependence? (iii) There is fluctuation energy in the turbulent wake and also outside it in the form of internal waves. How are kinetic and potential energy partitioned between the wake and wave field? How does the energy in wake waves compare with that in lee waves? (iv) Are results for the two cases at higher 
${{Re}} = 50\,000$ consistent with trends exhibited in the parametric study at the lower value of 
${{Re}} = 5000$?
A disk is used as the wake generator in body-inclusive LES, the numerics of which are laid out in § 2. Lee waves and their dependence on 
${{Fr}}$ is described in § 3. Characteristics of wake-generated internal gravity waves are described in § 4. Energy partition of the entire flow field is done to calculate the energies associated with lee waves, wake waves, mean wake and turbulent wake, and their variation with 
${{Fr}}$ is analysed in § 5. The influence of Reynolds number is briefly assessed in § 6. To conclude, a summary and a discussion are presented in § 7.
2. Methodology
2.1. Governing equations and numerical scheme
The wake of a disk is simulated for five different cases of stratification by solving the three-dimensional, incompressible, unsteady form of the conservation equations for mass, momentum and density. A high-resolution LES with the Boussinesq approximation for density effects is used. The disk, with diameter 
$D$ and thickness 
$0.01D$, is immersed perpendicular to a flow with velocity 
$U_{\infty }$. The equations are numerically solved in cylindrical coordinates but both Cartesian 
$(x,y,z)$ and cylindrical 
$(r,\theta,x)$ coordinates are used as appropriate in the discussion. Here, 
$x$ is streamwise, 
$y$ is spanwise and positive along 
$\theta =0^{\circ }$ and 
$z$ is vertical and positive along 
$\theta =90^{\circ }$ (figure 1). The density field 
$(\rho (\boldsymbol {x},t))$ is split into a constant reference density 
$(\rho _{0})$, the variation of the background 
$(\Delta \rho _{b}(z))$ and the flow induced deviation, 
$(\rho _{d}(\boldsymbol {x},t))$ so that 
$\rho (\boldsymbol {x},t)=\rho _{0}+\Delta \rho _{b}(z)+\rho _{d}(\boldsymbol {x},t)$. The Reynolds number of the flow, defined as 
${{Re}} = U_{\infty } D/\nu$ is 
$5000$. The different levels of stratification are quantified by the Froude number, 
$Fr = U_{\infty }/ ND$, which takes values of 
$0.5, 1, 1.5, 2$ and 
$5$ for the five cases. Here, 
$N$ is the buoyancy frequency given by 
$N = \sqrt {(-g/\rho _{0}) (\partial \Delta \rho _{b} (z) / \partial z) }$.
(a) Schematic showing the simulation set-up and domain. (b) Radial and streamwise grids visualised on the 
$y = 0$ plane. (c) Radial and azimuthal grids visualised on the 
$x = 0$ plane. (Every fifth gridline is shown, with blue lines representing the sponge region. Disk (not to scale) is shown in yellow.)
The filtered non-dimensional equations, using 
$D$, 
$U_{\infty }$ and 
$\rho _{0}$ as the characteristic length, velocity and density scales, respectively, are as follows
$$\begin{gather} \frac{\partial u_{i}}{\partial x_{i}} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial u_{i}}{\partial t} + u_{j}\frac{\partial u_{i}}{\partial x_{j}} ={-}\frac{\partial p}{\partial x_{i}} + \frac{1}{{\textit{Re}}} \frac{\partial}{\partial x_{j}}\Big[\left(1 + \frac{\nu_{sgs}}{\nu}\right) \frac{\partial u_{i}}{\partial x_{j}}\Big] - \frac{\rho_{d}}{Fr^{2}} \delta_{i3}, \end{gather}$$
$$\begin{gather}\frac{\partial \rho}{\partial t} + u_{j}\frac{\partial \rho}{\partial x_{j}} = \frac{1}{{\textit{Re}} Pr} \frac{\partial}{\partial x_{j}}\Big[\left(1 + \frac{\kappa_{sgs}}{\kappa}\right)\frac{\partial \rho}{\partial x_{j}}\Big], \end{gather}$$
where 
$u_{i}$ refers to the filtered non-dimensional velocities in the cartesian coordinate system for 
$i = 1, 2$ and 
$3$, respectively. The Smagorinsky eddy viscosity model with dynamic procedure (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) is employed to obtain 
$\nu _{sgs}$ wherein, the coefficient 
$C$ in the expression 
$\nu _{sgs} = C \tilde {\varDelta }^{2}|\tilde {S}|$ (
$\tilde {\varDelta }^{3}$ being the cell volume and 
$\widetilde {|S|}$ being the strain rate magnitude) is obtained dynamically using the least squares approach of Lilly (Reference Lilly1992) and the Lagrangian averaging method of Meneveau, Lund & Cabot (Reference Meneveau, Lund and Cabot1996). The Lagrangian averaging method takes the cumulative average of 
$C$ as a function of time, with more weight assigned to recent times in flow evolution. The Prandtl number (
$Pr=\nu / \kappa$) as well as the subgrid Prandtl number (
$Pr=\nu _{sgs} / \kappa _{sgs}$) is set to unity. Second-order central finite difference is used in space. A low storage third-order Runge–Kutta scheme for the advective terms and a Crank–Nicolson term for the viscous terms are used to advance the solution in time. The disk is represented using the immersed boundary method of Balaras (Reference Balaras2004) and Yang & Balaras (Reference Yang and Balaras2006).
2.2. Domain and statistics
The domain extends from 
$x/D = -L_{x}^{-} = -30$ to 
$x/D = L_{x}^+ =102$ in the streamwise direction and from 
$r/D = 0$ to 
$r/D = L_{r} =60$ in the radial direction (figure 1b,c). The number of grid points in the streamwise, azimuthal and radial directions are 
$N_{x} = 2176$, 
$N_{\theta } = 256$ and 
$N_{r} = 479$. The disk surface is resolved into 47 312 triangles to identify the fluid–solid interface. The streamwise and radial grids in the flow solver are structured and employ uniform stretching as distance from the disk increases. There are 73 points along the disk in the radial direction. Maximum values of 
$\Delta x/\eta$ and 
$\Delta r / \eta$, where 
$\eta = (\nu ^{3} / \varepsilon )^{1/4}$ is the Kolmogorov length scale, are used to ascertain adequate grid resolution. The chosen grid has 
$(\Delta x / \eta )_{max} = 4.84$ and 
$(\Delta r / \eta )_{max} = 5.43$ across all cases. The grid in the azimuthal direction is homogenous. Owing to the high grid resolution, the subgrid viscosity (
$\nu _{sgs}$) is not large. The instantaneous maximum value of subgrid viscosity fluctuates between 
$\nu _{sgs} / \nu = 2.5$ and 
$\nu _{sgs} / \nu = 5$ after statistical stationarity is achieved. The location of this maximum value is 
$0.5 < r/D < 1.5$ and 
$0< x/D < 2$. Beyond 
$x/D = 10$, 
$\nu _{sgs}/\nu < 1$. Statistics are collected by temporal averaging after the initial transient has subsided and statistical steady state has been reached. The time interval for averaging is 
$140 D / U_{\infty }$. For any dependent variable in the simulation, 
$\langle. \rangle$ represents the time average and superscript 
$'$ represents the fluctuation about that average. From this section onwards, all the variables will appear in their dimensional form and any non-dimensionalisation will be explicitly shown, e.g. 
$x/D$ or 
$w/U_{\infty }$.
2.3. Boundary conditions
For boundary conditions, the inlet (
$x/D = -30$) has a uniform free stream of magnitude 
$U_{\infty }$ in the positive 
$x$ direction, the outlet (
$x/D = 102$) has an Orlanski-type convective boundary condition (Orlanski Reference Orlanski1976), and the radial boundary (
$r/D = 60$) has a Neumann boundary condition for density as well as velocity. The centreline boundary (
$r = 0$) is handled by taking the average of the two symmetrically located ghost cells over the centreline to define the radial and azimuthal velocity components at the centreline. To prevent spurious reflection of waves back into the domain, sponge layers are used at the three boundaries of the computational domain. The radial sponge starts at 
$r/D = 45$ and goes up to the radial boundary 
$r/D = 60$, the inlet sponge is from 
$x/D = -30$ to 
$x/D = -20$ and the outlet sponge is from 
$x/D = 90$ to 
$x/D = 102$. Each of the sponge layers use a quadratic damping function of the form 
$f_{damp} = C ( (x - x_{spng})/(L_{x} - x_{spng}) )^{2} (\phi _{0} - \phi (x))$, where the variable 
$\phi (x)$ is relaxed to the target value 
$\phi _{0}$ in the sponge layer given by 
$x_{spng} < x < L_{b}$. For example, in the radial sponge, where 
$x$ in 
$f_{damp} (x)$ is equivalent to 
$r$, the sponge begins at 
$r_{spng}/D = 45$ and extends up to 
$L_{r}/D = 60$, and the strength of the damping function is 
$C = 5$.
3. Body-generated lee waves
In stratified environments, the wave component of the flow field consists of body-generated lee waves and wake-generated internal waves. The former is steady in time with respect to the body while the latter is unsteady. Figure 2 illustrates these two types of waves for the five values of Froude number by plotting the vertical velocity at the top half of the vertical centreplane. The steady lee waves are seen as the large band of waves spanning most of the domain while the unsteady waves are seen as a radially thinner band of irregular waves superposed on the lee waves. Gola et al. (Reference Gola, Nidhan, Ortiz-Tarin, Pham and Sarkar2023) studied the steady lee waves generated by a disk at 
${{Re}} = 5000$ focussing on how nonlinear stratification (pycnocline) as well as the relative position of the pycnocline with respect to the disk modifies the flow. The minimum background Froude number for the pycnocline cases was 
$Fr =1$ as was the constant value of 
$Fr$ in the linearly stratified baseline case.
Instantaneous contours of vertical velocity showing body-generated lee waves (steady in the simulation frame) on the 
$\theta = 90^{\circ }$, vertical centreplane for the five cases.
For the steady lee waves, Gola et al. (Reference Gola, Nidhan, Ortiz-Tarin, Pham and Sarkar2023) found that the wavelength in all cases and the amplitude in the linearly stratified 
${{Fr}} =1$ case could be deduced from the linear theory (with the value of 
$N$ in the far field where waves propagate) given by Voisin (Reference Voisin1991, Reference Voisin2003, Reference Voisin2007) and by modelling the disk plus its separation bubble as one body, specifically a Rankine ovoid. The potential flow past an ovoid is given by the superposition of a doublet (source/sink of strength 
$\pm m$ separated by distance of 
$2a$) and the uniform incident flow. The expression for the vertical velocity deduced by Ortiz-Tarin et al. (Reference Ortiz-Tarin, Chongsiripinyo and Sarkar2019) after adapting the linear theory to an ovoid-shaped body is
\begin{align} w(x,y,z) &\sim{-}\frac{mN}{{\rm \pi} U_{\infty} r_{xyz}} \sin{\theta} \cos{\psi} (1+\cot^{2}{\psi}\cos^{2}{\theta})^{{1}/{2}} \sin{\left(\frac{Na}{U_{\infty}} \cos{\psi} |\sin{\theta}| \right)} \nonumber\\ &\quad \times\sin{\left(\frac{N}{U_{\infty}} r_{xyz} |\sin{\theta}| \right)} , \end{align}
where 
$r_{xyz}=\sqrt {x^{2} + y^{2} + z^{2}}$, 
$\psi = \arctan (\sqrt {y^{2}+z^{2}}/x)$ and 
$\theta = \arctan {(z/y)}$ is the azimuthal angle in figure 1. Specialisation to the 
$y = 0$ plane leads to
\begin{equation} w(x,y=0,z) \sim{-}\frac{mN}{{\rm \pi} U_{\infty} r_{xz}} \cos{\psi} \sin{\Big(\frac{Na}{U_{\infty}} \cos{\psi}\Big)} \sin{\Big(\frac{N}{U_{\infty}} r_{xz}\Big)} , \end{equation}
where 
$r_{xz}=\sqrt {x^{2}+z^{2}}$ and 
$\psi = \arctan (z/x)$. The theoretical result, (3.2), will be subsequently compared with the simulation results.
As in the 
${{Fr}} =1$ case tested by Gola et al. (Reference Gola, Nidhan, Ortiz-Tarin, Pham and Sarkar2023), 
$m$ and 
$a$ are calculated here from the potential-flow solution for a Rankine ovoid of length 
$L_{ro}$ and cross-sectional diameter 
$D_{ro}$ fitted using the streamline 
$\langle u \rangle = 0.95 U_{\infty }$:
\begin{equation}
(L_{ro}^2 - 4a^2 )^2 =
\left( \frac{8ma}{{\rm \pi} U_{\infty}} \right) L_{ro} , \quad
D_{ro}^2 = \left( \frac{8ma}{{\rm \pi} U_{\infty}} \right)
\frac{1}{\sqrt{4a^2 + D_{ro}^2}}.
\end{equation} It can be deduced from the argument of the last sine term in (3.2) that the wavelength of the steady lee waves is 
$\varLambda / D = 2{\rm \pi} {{Fr}}$. The wavelength of the lee waves obtained from the simulations along the line 
$x=z, y=0$ is shown in figure 3, wherein it follows the trend 
$\varLambda = 6.7 {{Fr}}$, which is close to the asymptotic result of 
$\varLambda / D = 2 {\rm \pi}{{Fr}}$.
Wavelength of body-generated lee waves.
Note that (3.1) does not apply to 
${{Fr}} < O(0.1)$ when the flow is mostly around the disk and the region involved in wave generation shrinks from the crest to the dividing streamline at a depth of approximately 
$U_\infty /N$ from the crest. Thus, the singular limit 
${{Fr}} \rightarrow 0$ of (3.1) is inadmissible.
Figure 4 shows the vertical velocity 
$w$ for the five cases on the line 
$z = x, y =0$ as a function of distance from the disk centre, 
$s = \sqrt {x^{2} + z^{2}}$. Excellent agreement of the simulation results with linear theory, (3.2), is seen. Also plotted alongside in figure 4(b) are the defect velocity (
$U_{d} = U_{\infty } - \langle u \rangle$) contours. There is a case-dependent difference in the separation bubble. The length of the separation bubble increases with increasing 
${{Fr}}$ leading to a different set of values for 
$m$ and 
$a$ in (3.1) for each case. Figure 5 shows the comparison of the vertical velocity contours on the 
$r\unicode{x2013}\theta$ plane obtained from the simulation and linear theory (3.1) for 
$Fr = 1$ at 
$x/D = 10$ and 
$30$. The contours are seen to span a larger area in the 
$r\unicode{x2013}\theta$ plane as distance from the body is increased. Again, good agreement between the simulation and theory is found. The amplitude of the lee wave (maximum or minimum at 
$x/D = 20, r/D > 10$), which is plotted in figure 6(a), reveals excellent agreement between theory and simulation. The 
$w$-amplitude of the lee wave in (3.1) is proportional to the product of 
$m$ and 
$N$ so that the change in the size of the separation region behind the disk also plays a part in the resulting wave amplitude as is reflected by the dependence of 
$m$ in (3.1) on 
${{Fr}}$ (figure 6b). This results in a scaling that is different from 
$w/U_{\infty } \sim {{{Fr}}}^{-1}$, which would be the case if 
$m$ were assumed constant across 
${{Fr}}$.
Body-generated lee waves. (a,c,e,g,i) Vertical mean velocity plotted along the line 
$z=x, y=0$ on the top half of the vertical plane (
$s = \sqrt {x^{2} + z^{2}}$). (b,d,f,h,j) Defect velocity (
$U_{d}$) contours showing the increasing size of the fitted Rankine ovoid (yellow dashed line) with increasing 
${{Fr}}$.
Body-generated lee waves shown at two downstream locations for 
${{Fr}} =1$: vertical velocity contours compared with linear theory at 
$x/D = 10$ (a–c) and 
$30$ (d–f). Line plots compare 
$w$ at 
$\theta = 90^{\circ }$.
(a) Vertical mean velocity magnitude of the lee waves at local maxima/minima near 
$x/D = 20$ in figure 4(a). (b) Variation of 
$m$ and 
$a$ (obtained by fitting the resulting ovoid to the mean separation bubble) with 
${{Fr}}$.
Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) obtained the mean horizontal divergence magnitude (
$|\partial \langle w \rangle / \partial z|$) of the internal wave field in an experimental study of waves emitted by towing a sphere, 
$6:1$ cylinder and 
$6:1$ spheroid in a stratified fluid. Measurements of 
$u(x,y)$ and 
$v(x,y)$ on a horizontal plane were used to infer 
$|\partial \langle w \rangle / \partial z|$. For comparison, the mean horizontal divergence magnitude at the local maxima/minima for 
$w$ at 
$x/D = 10, r/D = 3$ for the disk simulations is plotted along with the data from Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) in figure 7. Note that the chosen locations of the maxima/minima in the simulations are somewhat different from those in Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018), reflecting the difference between simulation and experimental methodologies. The simulation results for the disk are comparable to the experimental results of Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) for other axisymmetric bodies.
Comparison of mean divergence magnitude of the lee-wave field between the present disk simulation and the laboratory results reported by Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018).
4. Wake-generated internal gravity waves
The wake-generated internal waves are manifested in the simulation reference frame as the unsteady component of the wave field. To visualise these waves without their steady counterpart, instantaneous contours of the fluctuation (
$\partial w' / \partial z$) in vertical divergence on the top half of the vertical plane is plotted in figure 8 for all the five cases. In a time series of these contours, these waves are seen to advect downstream (see Supplementary material ‘Movie1.mp4’ available at https://doi.org/10.1017/jfm.2024.717), unlike the steady lee waves. When viewed in a reference frame moving with the free stream (see Supplementary material ‘Movie2.mp4’), these waves move slowly towards the left while growing in the vertical. In other words, the wave phase speed is less than the free stream velocity magnitude.
Instantaneous contours of vertical derivative of vertical velocity fluctuation showing wake-generated internal gravity waves on the 
$\theta = 90^{\circ }$ plane for the five cases.
Across the five cases, three features of these internal waves can be noted: (i) the inclination angle with the vertical (
$\varTheta$) is in a narrow range (
$35^{\circ }\unicode{x2013}40^{\circ }$), specifically 
$35^{\circ }, 39^{\circ }, 35^{\circ }, 40^{\circ }$ and 
$37^{\circ }$ in the order of increasing 
${{Fr}}$; (ii) the wavelength increases with increasing 
${{Fr}}$; (iii) the vertical extent of propagation with respect to the body decreases with increasing 
${{Fr}}$. The first feature, the clustering of the waves around the same vertical angle, was explained using the argument of maximising the vertical component of group velocity for a fixed horizontal wavenumber by Taylor & Sarkar (Reference Taylor and Sarkar2007). The vertical group velocity for internal gravity waves in stratified flow without rotation is given by
\begin{equation} c_{gz} ={-}\frac{N}{|\boldsymbol{k}|} \cos{\varTheta} \sin{\varTheta} ={-}\frac{N}{|k_{x}|} |\cos{\varTheta}| \cos{\varTheta} \sin{\varTheta}, \end{equation}
where 
$|k|$ is the wavenumber magnitude and 
$k_{x}$ is the horizontal wavenumber. Note that for our simulations, 
$\varTheta$ is an obtuse angle because the waves propagate in the negative 
$x$ direction as discussed earlier. It can be shown that the expression on the right-hand side of (4.1) is maximum when 
$\varTheta \approx 145^{\circ }$ or 
$35^{\circ }$ from the vertical, which is close to the angles observed in the simulations.
It is evident from figure 8 that the wavelength of the wake waves increases with 
${{Fr}}$. This dependence is quantified by computing the wavelength (
$\lambda$), normalised by 
$D$, from wave packets in the instantaneous 
$\partial w' / \partial z$ fields shown in figure 8. The wave packets are chosen to be 4–5 wavelengths long in the NEQ region away from the wake to avoid turbulence but also close enough to get a strong clean signal (this corresponds to the region 
$5 < z < 15, 5 < Nt < 30$ depending on 
${{Fr}}$). Figure 9 shows 
$\lambda$ plotted against 
${{Fr}}$. There is reasonable agreement with the scaling 
$\lambda \propto {{Fr}}^{1/3}$ (dashed line) found by Abdilghanie & Diamessis (Reference Abdilghanie and Diamessis2013) and Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018). The horizontal width (
$L_{Hk}$) of the kinetic energy profile of the disk wake varies as 
$L_{Hk} \propto x^{1/3}$ (Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020). Assume that the energetic eddies responsible for wave generation scale with 
$L_{Hk}$ and that the internal wave radiation becomes significant only when 
$Nt$ reaches a specific value. Then, the distance at which the internal wave radiation commences scales as 
${{Fr}}$, the value of 
$L_{Hk}$ at that location varies as 
${{Fr}}^{1/3}$ and, therefore, 
$\lambda \propto {{Fr}}^{1/3}$. Indeed, it will be shown in § 5.3 that the wake-wave potential energy (as a fraction of turbulent kinetic energy of the wake) becomes significant at 
$Nt \approx 5$ and then does not change appreciably until 
$Nt \approx 30$ when it decreases.
Wavelength of wake-generated internal gravity waves.
Lastly, to explain the trend in vertical extent of the waves, the contours of the r.m.s. buoyancy term 
$b_{rms} = g \rho _{d_{rms}} / \rho _{0}$ are plotted in figure 10(a–e), thereby providing a good visualisation of the envelope of the wake-generated internal waves. When the vertical extent of this envelope (quantified by a small value of 
$b_{rms} = 0.0015 U_{\infty } ^{2} / D$) is plotted against 
$Nt$ in figure 10(f), a good collapse is seen, showing that the vertical extent of propagation is the same for these waves in buoyancy time units. For example, 
$x = 40$ corresponds to 
$Nt = 80$ for 
$Fr = 0.5$ and 
$Nt = 8$ for 
$Fr = 5$, accounting for the difference in wave envelope between these cases at the same 
$x$.
(a–e) Wake wave envelope visualised by the root mean square (r.m.s.) of buoyancy term for the five cases. (f) Contour lines of 
$b_{rms} = 0.0015 U_{\infty }^{2}/D$ plotted with 
$Nt$.
5. Wake and wave energetics
A turbulent stratified flow has both kinetic and potential energy with each having a mean and a fluctuating component. Furthermore, the mean and the fluctuation can be divided into a turbulent wake component and a wave part surrounding the wake component. In this section, the variation of the energy partition among its components will be diagnosed. It will be shown that, across the five levels of stratification, there is a systematic variation of mean relative to turbulent components as well as potential versus kinetic energy. The definitions of various energy components are as follows:
mean kinetic energy,
\begin{equation} MKE = \tfrac{1}{2} ((\langle u \rangle-U_{\infty})^{2} + \langle v \rangle^{2} + \langle w \rangle^{2}); \end{equation}mean potential energy,
\begin{equation} MPE = \frac{g^{2} \langle \rho_d \rangle^{2}}{2 \rho_{0}^{2} N^{2}}; \end{equation}turbulent kinetic energy,
\begin{equation} TKE = \tfrac{1}{2} (\langle u'^{2} \rangle + \langle v'^{2} \rangle + \langle w'^{2} \rangle); \end{equation}turbulent potential energy,
\begin{equation} TPE = \frac{g^{2} \langle \rho_d'^{2} \rangle}{2 \rho_{0}^{2} N^{2}} . \end{equation}
The various plots in this section compare the evolution of these energy components, integrated over the 
$r\unicode{x2013}\theta$ plane, as a function of streamwise distance (
$x$). Let 
$A$ be the complete 
$r\unicode{x2013}\theta$ plane in the simulation domain at any 
$x$. If 
$C$ is a closed curve on 
$A$ enclosing the area 
$A_{C}$, the following three integrals of a quantity 
$E(r, \theta, x)$ can be defined as a function of 
$x$ (see figure 11 for two examples of the curve 
$C$ chosen to separate the wave from the wake field): Total integral over the complete 
$r\unicode{x2013}\theta$ plane,
\begin{equation} E_{TI}(x) = \int_{A} E(r, \theta, x) r \,{\rm d}r \,{\rm d}\theta; \end{equation}
inside integral over the area bounded by 
$C$ on 
$r\unicode{x2013}\theta$ plane (note that 
$C$ varies as a function of 
$x$),
\begin{equation} E_{II}(x) = \int_{A_{C}(x)} E(r, \theta, x) r \,{\rm d}r \,{\rm d}\theta; \end{equation}
outside integral, which is computed over the area external to that bounded by 
$C$,
\begin{equation} E_{OI}(x) = \int_{A - A_{C}(x)} E(r, \theta, x) r \,{\rm d}r \,{\rm d}\theta; \end{equation}
so that, 
$E_{TI}(x) = E_{II}(x) + E_{OI}(x)$.
Examples of curve 
$C$ (here an ellipse shown in dashed yellow line) used to separate the wave from the wake contribution: (a) mean kinetic energy contours for 
$Fr = 1.5$ at 
$x/D = 20$ as separated by 
$C$ based on half-width of wake 
$U_{d}$ profile; (b) turbulent kinetic energy contours for 
$Fr = 1.5$ at 
$x/D = 20$ as separated by 
$C$ based on half-width of wake turbulent kinetic energy profile.
5.1. Mean kinetic energy and mean potential energy
The curve 
$C$ in (5.6), and (5.7) will be chosen to be at the edge of the wake so that the outside integral serves as a surrogate for the wave energy and the inside integral for the wake energy. Specifically, an ellipse with its semimajor and semiminor lengths equal to the horizontal and vertical half-width of the mean defect velocity profile is chosen for each 
${{Fr}}$ case; these lengths vary as a function of 
$x$. Thus, using (5.5), (5.6) and (5.7), the wake and wave contributions based on the inside and outside regions of the ellipse can be calculated.
Figure 12(a–c) shows the various integrals of 
$MKE$ for all five cases. For cases with lower 
${{Fr}}$, 
$MKE_{TI}$ (figure 12a) is initially low because of shorter separation zones (see figure 4). However, 
$MKE_{TI}$ also exhibits slower decay with increasing 
$x$ for the higher-
${{Fr}}$ cases, reflecting the fact that stratification prolongs wake lifetime, so that the 
${{Fr}}$ variation in 
$MKE_{TI}$ eventually reverses at later time. There is significant mean kinetic energy in the outside wave region (figure 12c), which can be attributed to the lee waves, except at the highest simulated 
${{Fr}} = 5$ case with weak lee waves. Moving to mean potential energy shown in the bottom row, comparison of figure 12(e) and figure 12(f) shows that the outside contribution (
$MPE_{OI}$ ) owing to steady lee waves dominates over the wake contribution (
$MPE_{II}$) over the entire streamwise extent of the flow. Evidently, the turbulent wake does not lead to significant mean distortion of the isopycnals although there is fluctuating distortion leading to turbulent kinetic energy as will be shown shortly. Beyond 
$x = 10$, both 
$MKE_{TI}$ and 
$MPE_{TI}$ decrease with increasing 
${{Fr}}$. The decrease is quite substantial from 
${{Fr}} = 2$ to 
${{Fr}} = 5$. Note that the energy values are normalised using 
$U_\infty$ and 
$D$. Thus if the increase of 
${{Fr}}$ is due to a higher flow speed under the same background 
$N$, the 
$Fr$-related decrease in the dimensional value of energy would be less.
Area-integrated mean energies for all five cases. Total integral (
$TI$) of the mean energy is shown as well as its partition into wake or inside integral (
$II$) and wave or outside integral (
$OI$).
To summarise, the distribution of mean energy associated with stratified flow past a disk is such that almost all of the mean potential energy component resides outside the turbulent wake in the lee wave field. At 
${{Fr}} = 2$ or less, most of mean kinetic energy resides in the lee wave field in contrast to the higher 
${{Fr}} = 5$ case where it is the wake that carries most of the mean kinetic energy.
5.2. Turbulent kinetic energy and turbulent kinetic energy
Figure 13 shows the various fluctuation energies for all five cases. The fluctuation energy can be divided into the turbulent energy of the wake and the outside fluctuation energy carried by wake-generated internal waves. This is achieved by defining an ellipse with its semimajor and semiminor axes equal to the horizontal and vertical half-width of the turbulent kinetic energy profile, respectively. The various contributions to the turbulent component are plotted similar to what was shown for the mean component in figure 12. Looking at the 
$TKE_{TI}$ profiles in figure 13(a), it is evident that the general 
${{Fr}}$-dependence seen in figure 12(a) has been flipped, i.e. the total 
$TKE$ is larger for the cases with larger 
${{Fr}}$. Furthermore, the cases with higher 
${{Fr}}$ also have stronger turbulence as shown by the interior integral 
$TKE_{II}$ in figure 12(b). The exception is the case of 
$Fr = 0.5$. Here, both 
$TKE_{II}$ (interior wake turbulence) and 
$TKE_{OI}$ (external to the wake) flatten at 
$x \approx 20, Nt \approx 40$ and then increase slightly before decaying again. The increase of wake fluctuation energy in the late NEQ/early quasi-two-dimensional regime of the 
${{Fr}} = 0.5$ case is a result of the wake structure becoming increasingly two-dimensional and its flapping in the horizontal, similar to what was found for 
${{Fr}} < 1$ sphere wakes at 
${{Re}} = 3700$ by Pal et al. (Reference Pal, Sarkar, Posa and Balaras2016).
Area integrated fluctuation energies for all five cases. The total integral of the fluctuating energy is shown as well as its partition into wake or inside integral and wave or outside integral.
The 
$TPE_{TI}$ for 
$Fr = 0.5$ also shows slight stalling at 
$x/D \approx 10$ before decaying again at 
$x/D \approx 30$, although this is likely the result of internal wave activity in the NEQ regime as can be seen from increase in external fluctuation energy, 
$TPE_{OI}$ in figure 13(f). Comparing figure 13 with figure 12 shows that 
$Fr = 5$ is qualitatively different than the other four cases, e.g. it has a larger potential energy in the fluctuations than in the mean.
5.3. Ratio of kinetic and potential energies
Figure 14 shows the mean and fluctuation energy ratios (
$MPE_{TI}/MKE_{TI}$ and 
$TPE_{TI}/TKE_{TI}$), respectively, plotted as a function of 
$x/D$ as well as 
$Nt$. For the mean energies, the ratio 
$MPE_{TI}/MKE_{TI}$ approaches a constant value after 
$Nt \approx 5$. This constant 
$O(1)$ value has a small variation for the lower 
${{Fr}}$ cases, e.g between 
$0.9$ at 
${{Fr}} = 0.5$ to 0.6 at 
${{Fr}} = 2$. However, there is a large decrease of the ratio to 0.2 at 
${{Fr}} = 5$, pointing to a decrease in the importance of the lee wave field at high 
${{Fr}}$. Insofar as the fluctuation energy ratio, 
$TPE_{TI}/TKE_{TI}$, it increases as the wake evolves to approach a plateau at 0.6–0.7 in the region 
$5 < Nt < 30$ (figure 14d). This plateau lies within the NEQ regime where the wake adjusts to the background stratification and the activity of the wake-generated internal waves is the highest.
Ratio of potential and kinetic energy for all the five cases. This ratio is plotted for the mean and turbulent components as a function of 
$x/D$ (a,c) and 
$Nt$ (b,d).
5.4. Comparison of body lee wave energy with wake wave energy
Previous laboratory measurements (Brandt & Rottier Reference Brandt and Rottier2015; Meunier et al. Reference Meunier, Le Dizès, Redekopp and Spedding2018) at selected spatial locations have shown that, while the body wave field is more energetic than the wake wave field at lower 
${{Fr}}$, there is a crossover to dominance of wake wave energy at sufficiently large 
${{Fr}}$. The wave energy is computed from the outside integral terms, i.e. 
$E_{body \ wave} = MPE_{OI} + MKE_{OI}$ and 
$E_{wake \ wave} = TPE_{OI} + TKE_{OI}$. Figure 15 compares the downstream evolution of 
$E_{body \ wave}$ and 
$E_{wake \ wave}$. Unlike the behaviour in the lower-
${{Fr}}$ cases, the energy in the wake-generated internal waves at at 
${{Fr}} = 5$ surpasses the energy in steady lee waves. Thus, in the disk wake, a crossover to the dominance of wake waves occurs at a value of 
${{Fr}}$ between 2 and 5.
(a) Total energy in body lee waves (b). Total energy in wake-generated internal waves.
6. Comparison with disk wakes at 
$Re = 50\,000$
In this section, internal waves from the data (Chongsiripinyo & Sarkar Reference Chongsiripinyo and Sarkar2020) for a disk wake at a higher 
${{Re}} = 50\,000$ and with 
${{Fr}} = 2, 10$ are quantified and the results compared with the 
${{Fr}}$ dependencies discussed in the previous sections.
6.1. Structure of body-generated lee waves and wake-generated internal gravity waves
Lee waves for 
${{Re}} = 5000$ and 
${{Re}} = 50\,000$ at 
${{Fr}} = 2$ are compared in figure 16 by plotting vertical velocity (
$w$) contours on the 
$\theta = 90^{\circ }$ plane and also the value of 
$w$ along the line 
$z=x,y=0$. The influence of 
${{Re}}$ on the wavelength and spatial distribution of the wave field is negligible. The lee waves for the 
${{Re}} = 50\,000$ case at 
${{Fr}} = 10$ (not shown here) are weaker than at 
$Fr =2$ and have longer wavelength, in agreement with the linear theory result of § 3.
Lee wave comparison for 
$Re = 5000$ and 
$Re = 50\,000$ at 
$Fr = 2$: (a,b) instantaneous vertical velocity contours; (c) vertical velocity on the line 
$z=x, y=0$ (
$s = \sqrt {x^{2} + z^{2}}$).
The influence of 
${{Re}}$ on wake waves is illustrated by figure 17(a,b) that shows 
$\partial w' / \partial z$ on the vertical centre half-plane (
$\theta = 90^{\circ }$). Similarity is observed in the inclination of the waves (
$\varTheta \approx 40^{\circ }$) as well as their wavelength (
$\lambda /D \approx 4.5$). The 
${{Re}} = 50\,000$ wake is more turbulent than at 
${{Re}} = 5000$ as suggested visually by figure 17(b) with enhanced small-scale features in the wake core relative to 17(a).
Wake-generated internal wave comparison for (a) 
$Re = 5000$ and (b) 
$Re = 50\,000$ at 
$Fr = 2$.
The wake wave amplitude, measured by 
$(\partial w' / \partial z)_{rms}$, is compared between 
${{Re}} = 5000$ and 
${{Re}} = 50\,000$ in figure 17(c,d). Proximity to the wake leads to a higher difference in wave amplitude, as is observed at 
$z/D = 2$ relative to 
$z/D = 4$. The scales of wake structures that give rise to wake waves are more broadband for 
${{Re}} = 50\,000$ as compared with 
${{Re}} = 5000$ and at higher 
${{Re}}$, waves are less affected initially by viscosity.
Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) report an empirical scaling of 
${{Re}}^{0.4}$ for the wake wave amplitude and Abdilghanie & Diamessis (Reference Abdilghanie and Diamessis2013) report that wave momentum flux varies as 
$\langle uw \rangle \sim {{Re}}^{0.25}$. Since the increase of wave amplitude with 
${{Re}}$ depends on streamwise location in the present simulations, we refrain from proposing a power law dependence on 
${{Re}}$.
6.2. Plane-integrated wave energy
Analysis of the 
${{Re}} = 5000$ results showed that when 
$Fr$ is increased, the total energy in wake waves progressively increases while the total energy in body waves decreases. Those results also showed that at 
$Fr = 5$, the wake wave energy overtook the body wave energy. Figure 18(a,b) shows the downstream variation of body wave energy (
$E_{body \ wave} = MPE_{OI} + MKE_{OI}$) and wake wave energy (
$E_{wake \ wave} = TPE_{OI} + TKE_{OI}$), after including the additional results for 
${{Re}} = 50\,000$ at 
$Fr = 2, 10$ to the 
${{Re}} = 5000$ cases. The streamwise evolution of 
$E_{body \ wave}$ at 
${{Fr}} = 2$ is qualitatively similar between the two values of 
${{Re}}$ although the values are slightly larger at the higher 
${{Re}}$. The 
${{Fr}} = 10$ case in figure 18(a) has the weakest stratification and also the smallest body wave energy among all the simulations, which is consistent with the trend at 
${{Re}} = 5000$. Turning to the wake wave energy as a function of 
$x/D$, the 
${{Fr}} = 2$ case shows consistently higher wave energy at 
${{Re}} = 50\,000$ relative to the lower-
${{Re}}$ case.
(a) Total energy in body lee waves. (b) Total energy in wake waves. (c) Ratio of total turbulent potential energy and turbulent kinetic energy.
When the ratio 
$TPE_{TI}/TKE_{TI}$ is plotted against 
$Nt$ in figure 18(c), now also including the data for 
$Re = 50\,000$, the two new cases also show that the ratio increases as a function of 
$Nt$ until reaching a maximum value of 
$0.6\unicode{x2013}0.7$ in the NEQ region, 
$5 < Nt< 30$. This implies that the balance between the turbulent energies for a stratified flow evolves in a similar fashion not only across different 
$Fr$ for a fixed 
$Re$, but also for higher 
$Re$, pointing towards a potential universality of this statistic for stratified flows.
7. Summary and discussion
Body-inclusive LES are used to study the effect of changing stratification levels in flow past a disk. The focus is on internal gravity waves. Systematic investigation of internal waves over a range of 
${{Fr}}$ is new insofar as body-inclusive simulations and adds to our knowledge base from previous experiments, temporal-model simulations and theory. Linear stratification at five different levels is considered: 
${{Fr}} = 0.5, 1, 1.5, 2$ and 
$5$ and the flow is quantified up to a relatively long streamwise distance of 
$x/D = 80$. The Reynolds number based on the free stream velocity (
$U_\infty$) and disk diameter (
$D$) is fixed at 
$5000$ and, thus, variability across these cases can be attributed to 
${{Fr}}$ alone. The dataset of Chongsiripinyo & Sarkar (Reference Chongsiripinyo and Sarkar2020) for a disk at 
${{Re}} = 50\,000$ at 
$Fr =2,10$ is also analysed for studying Reynolds number dependence.
The questions posed in § 1.5 are answered and discussed below.
7.1. Prediction of lee waves
With regards to the steady lee waves, the non-dimensional amplitude (
$w$) of the vertical velocity decreases with increasing 
${{Fr}}$. All five cases showed excellent agreement with the linear theory for steady lee waves by Voisin (Reference Voisin2007) when an equivalent body including the separation bubble was used to calculate the potential flow solution. As to the separation bubble, its length differed considerably among the 
${{Fr}} = 0.5, 1$ and 
$1.5$ cases but not so much between 
${{Fr}} = 2$ and 
${{Fr}} = 5$. According to linear theory (3.3a,b), 
$w/ U_\infty \propto Fr ^ {-1}$ if 
$m$ is constant and the change of 
$m$ adds some variability. The wavelength (
$\varLambda$) was found to increase linearly with 
${{Fr}}$ according to 
$\varLambda / D \sim 6.7 {{Fr}}$, close to the asymptotic relation of 
$\varLambda / D = 2{\rm \pi} {{Fr}}$. The amplitude of mean 
$\partial w /\partial z$, which was obtained by Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) in laboratory experiments, was compared with the present simulation results and good mutual agreement was found. The theory proposed by Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) involves a forcing of the momentum conservation equation (to model the drag) as the leading-order contribution to wave generation while the theory of Voisin (Reference Voisin2007) that was adapted for the present application is different and involves a volumetric forcing of the continuity equation. Our adaptation of the theory involves the inclusion of the separation bubble, whose size is related to and generally increases with the drag force, as part of the wave forcing. Thus, although not equivalent, both the theory used here and that used by Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018) lead to higher wave amplitude for blunter body shapes.
7.2. Properties of wake waves
Unlike the body-attached lee waves, the wake waves advect downstream with respect to the disk. In a reference frame moving with the free stream, the wake waves move upstream consistent with the notion that the wake has a velocity deficit with respect to the free stream. A narrow band (
$\varTheta = 35^{\circ }\unicode{x2013}40^{\circ }$ with respect to the vertical) is found for the wave propagation angles for all the five 
${{Fr}}$ cases. Narrow-band emission from turbulent shear flows has been found previously, e.g. from a shear layer (laboratory experiments by Sutherland & Linden (Reference Sutherland and Linden1998) and DNS by Pham et al. (Reference Pham, Sarkar and Brucker2009)), a turbulent Ekman layer (Taylor & Sarkar Reference Taylor and Sarkar2007) and temporal-model simulations of a wake (Abdilghanie & Diamessis Reference Abdilghanie and Diamessis2013). The wave propagation angles observed here are consistent with the analytical result of 
$35^\circ$ that follows from the vertical group velocity maximisation (equivalently, viscous decay minimisation) criterion given by Taylor & Sarkar (Reference Taylor and Sarkar2007). It is worth noting that Abdilghanie & Diamessis (Reference Abdilghanie and Diamessis2013) found that the narrow band in their lower 
${{Re}}$ cases expanded to a wider band of inclination angles (
$26^{\circ }\unicode{x2013}50^{\circ }$) at higher 
${{Re}}$; a possible reason is the reduced effect of viscosity on the waves.
Other characteristics of the wake waves were studied as well. The wavelength of these waves increases with an increase in 
${{Fr}}$, agreeing well with the 
$\lambda \propto {{Fr}}^{1/3}$ scaling found by Abdilghanie & Diamessis (Reference Abdilghanie and Diamessis2013) and Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018). Also, the wavelength remains approximately constant until the end of the test domain at 
$x = 80$. The horizontal width of the turbulent kinetic energy profile increases as 
$x^{1/3}$ and the wave potential energy as a fraction of turbulent kinetic energy becomes significant at 
$Nt = Nx/U_\infty \approx 5$ for all five cases. Thus, the streamwise distance at which the internal wave burst becomes significant scales as 
${{Fr}}$, the length scale of the energetic flow structures (same order as the turbulent kinetic energy profile thickness) at that point scales as 
${{Fr}} ^{1/3}$, and so does the corresponding wavelength of the radiated internal wave.
The envelope of the wake waves is quantified by choosing the r.m.s. value of the buoyancy at the envelope boundary to be a small fixed value. The width of the thus defined envelope at a given 
$x$ expands with decreasing 
${{Fr}}$. Nevertheless, self-similar evolution is found when the wave envelope is plotted as a function of the buoyancy time scale (
$Nt$). The use of 
$Nt$ in place of 
$x/D$ is necessary because the evolution time scale of wake waves is a buoyancy effect that depends on stratification (
$1/N$) and not on advection (
$D/U_{\infty }$).
7.3. Wave and wake energetics
Because of the presence of both a turbulent wake and internal gravity waves external to the wake, it makes sense to partition the energies in the flow into their respective interior (wake) and exterior (wave) contributions. This is done by calculating the area integrals of the energy which are then segregated into their wake and wave parts. The interior wake part of the energy is defined to be inside an ellipse whose dimensions are determined by the 
$x$-dependent values of the vertical and horizontal half-widths of the wake and the exterior wave part is the remaining portion. The wake–wave decomposition is applied to the kinetic energy (both the mean component, 
$MKE$, and the turbulent component, 
$TKE$) and the potential energy (mean component, 
$MPE$, and turbulent component, 
$TPE$).
Comparison of the body lee wave energy (outside integrals of mean energies) and wake wave energy (outside integrals of fluctuation energies) shows that lee waves are a dominant feature for lower 
${{Fr}} \sim O(1)$ while unsteady wake waves become more significant as 
${{Fr}}$ is increased. For the disk wake, 
${{Fr}} = 2$ is below the crossover and 
${{Fr}} =5$ is above the crossover point at which both types of waves have similar energetic importance. The crossover from lee-wave dominance to wake-wave dominance at sufficiently large 
${{Fr}}$ has been seen at selected measurement stations in previous experimental studies, e.g. Brandt & Rottier (Reference Brandt and Rottier2015) and Meunier et al. (Reference Meunier, Le Dizès, Redekopp and Spedding2018).
For the mean flow, it is found that although mean kinetic energy gets its contribution from both the wake defect and the lee waves, mean potential energy is primarily in the lee waves with 
${{Fr}} = 0.5$ having the highest value. Beyond the near wake, the wake contribution to mean kinetic energy increases with a decrease in 
${{Fr}}$, i.e. strengthening stratification, consistent with stratification prolonging the lifetime of the wake. The ratio of total mean potential energy to total mean kinetic energy reaches a case-dependent constant value beyond 
$Nt \approx 5$ that decreases with increasing 
${{Fr}}$. The ratio takes higher values of 0.6 –1 at 
${{Fr}}$ between 0.5 and 2 – the regime where lee waves dominate over wake waves.
For the fluctuation energies, 
${{Fr}} = 5$ has the highest wake turbulent kinetic energy. Notably, the ratio of turbulent potential energy and turbulent kinetic energy exhibits a universal trend based on all the cases considered in this study, where it is seen to have a constant value of 0.6 to 0.7 for 
$5 < Nt < 40$, which coincides with the NEQ regime, wherein the wake-generated internal wave activity is the highest as part of the adjustment of the wake turbulence to the background stratification.
7.4. Effect of increasing Reynolds number
Using the dataset of Chongsiripinyo & Sarkar (Reference Chongsiripinyo and Sarkar2020), it was established that the increase in Reynolds number has little effect on lee wave properties such as wavelength and amplitude. However, the amplitude of wake waves, measured by pointwise values of 
$(\partial w' / \partial z)_{rms}$, increases with 
${{Re}}$. The relative enhancement of wave energy is higher at locations closer to the wake. The increase of wake wave amplitude with increasing 
${{Re}}$ is consistent with previous results (Abdilghanie & Diamessis Reference Abdilghanie and Diamessis2013; Meunier et al. Reference Meunier, Le Dizès, Redekopp and Spedding2018). Importantly, the ratio of cross-sectional areas of turbulent potential energy and turbulent kinetic energy exhibits some universality. Specifically, the peak value of this ratio takes the value 0.6–0.7 in the NEQ regime for all 
${{Re}}$ and 
$Fr$ cases in this study.
7.5. Limitations and outlook
Internal waves are examined in this study with body-inclusive simulations for 
${{Fr}}$ between 0.5 and 10. It is of interest to expand the parameter space to 
${{Fr}} > 10$, a regime with initially weak buoyancy effects, and test the present scaling laws on wake wave structure and energetics. At such high values of 
${{Fr}}$, a large streamwise domain with 
$x/D = {{Fr}} Nt$ is required to capture the regime (
$5 < Nt < 40$) that has energetic internal waves. The robustness of the present results at higher 
${{Re}}$ also deserves further study. A pure body-inclusive approach with the long streamwise domains required at high 
${{Fr}}$ or the fine grid spacing required at high 
${{Re}}$ would incur excessively high computational cost. The hybrid approach, which combines a body-inclusive simulation with a body-exclusive simulation (either spatial or temporal), is a good option to extend the 
${{Fr}}$ and 
${{Re}}$ range of simulations.
Supplementary material
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.717.
Funding
The authors gratefully acknowledge the support of Office of Naval Research, grant no. N00014-20-1-2253.
Declaration of interests
The authors report no conflict of interest.
Open access
