2.1. Governing equations and numerical scheme

The wake of a disk is simulated by solving the 3-D 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 appropriately used 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 ${\it Re} = U_{\infty } D/\nu$ (where $\nu$ is the kinematic viscosity) is $5000$. Since the stratification is non-uniform, the Froude number defined by $Fr (z) = U_{\infty }/N(z) D$ is variable and its case-dependent behaviour is described in § 2.2. The minimum value of the Froude number ($Fr_{min}$) is set to $1$ for all the cases.

Figure 1. Schematic showing the simulation set-up and domain.

The filtered non-dimensional equations are as follows:

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

(2.2)\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}{{\it Re} } \frac{\partial}{\partial x_{j}}\left[\left(1 + \frac{\nu_{sgs}}{\nu}\right) \frac{\partial u_{i}}{\partial x_{j}}\right] - \frac{\rho_{d}}{Fr_{min}^{2}} \delta_{i3}, \end{gather}

(2.3)\begin{gather} \frac{\partial \rho}{\partial t} + u_{j}\frac{\partial \rho}{\partial x_{j}} = \frac{1}{{\it Re} Pr} \frac{\partial}{\partial x_{j}}\left[\left(1 + \frac{\kappa_{sgs}}{\kappa}\right)\frac{\partial \rho}{\partial x_{j}}\right], \end{gather}

where $u_{i}$ refers to the filtered velocities in the $x$, $y$ and $z$ directions for $i = 1, 2$ and $3$, respectively, $\nu _{sgs}$ and $\nu$ in (2.2) are the subgrid-scale kinematic viscosity and the kinematic viscosity, respectively, while $\kappa _{sgs}$ and $\kappa$ in (2.3) are the subgrid-scale diffusion coefficient and the diffusion coefficient, respectively. The Prandtl number defined by $Pr=\nu / \kappa$ is set to one. No major qualitative influence of $Pr$ was found in the wake study by de Stadler, Sarkar & Brucker (Reference de Stadler, Sarkar and Brucker2010) who varied $Pr$ between 0.2 and 7. The dynamic eddy viscosity model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) is employed to obtain $\nu _{sgs}$, following the implementation of Chongsiripinyo & Sarkar (Reference Chongsiripinyo and Sarkar2020), and the subgrid Prandtl number is also set to unity to obtain $\kappa _{sgs}$. The eddy viscosity model calculates the values of eddy viscosity in the vicinity of the solid boundary by following a similar procedure as used for the velocity field, i.e. by enforcing Dirichlet boundary conditions on virtual points on the solid boundary and then using a linear reconstruction scheme to modify the value at the forcing points near the solid boundary (Balaras Reference Balaras2004). The parameters used for non-dimensionalization are as follows: free-stream velocity ($U_{\infty }$) for velocity, disk diameter ($D$) for length, advection time ($D/U_{\infty }$) for time, dynamic pressure ($\rho U_{\infty }^{2}$) for pressure and characteristic change in background density ($-D (\partial \Delta \rho _{b} / \partial z) |_{max}$) for the flow induced density deviation.

The periodicity in the azimuthal direction is leveraged in solving the discretized pressure Poisson equation in the predictor step, reducing it to a pentadiagonal system of linear equations, which is then solved using a direct solver (Rossi & Toivanen Reference Rossi and Toivanen1999). The disk is represented by the immersed body method of Balaras (Reference Balaras2004) and Yang & Balaras (Reference Yang and Balaras2006). At the inlet boundary, a uniform stream of velocity ($U_{\infty }$) is imposed while an Orlanski-type convective boundary condition is used for the outflow (Orlanski Reference Orlanski1976). A Neumann boundary condition is imposed at the radial boundary of the domain for density as well as the three velocity components. To prevent spurious reflection of waves back into the domain, sponge layers are used at the inlet (streamwise length $10D$), outlet (streamwise length $10D$) and the cylindrical walls (radial length $15D$) of the domain to gradually relax the velocities and the density to their respective background values at the boundaries.

2.2. Density profiles

The profiles chosen for the variation of background density in the five cases are shown in figure 2. Except for the benchmark 111 case, the profiles have non-uniform $N(z)$ and $Fr(z) = U/N(z)D$. Profile 111 with linear stratification has a constant $Fr (z)$, which is the conventional body-based Froude number of $Fr = U_{\infty }/ND = 1$. The maximum value of the density gradient is the same for all four nonlinear profiles and corresponds to $Fr_{min} =1$. In profile 614 (figure 2b), a hyperbolic tangent function is used to bridge two linearly stratified regions (6 and 4 represent the value, rounded to an integer, of the local $Fr$ in the linearly stratified regions). Profile I1I (I standing for ‘Infinity’) is a complete hyperbolic tangent profile between two constant-density regions, therefore having infinite Froude number far above and below the body.

Figure 2. Variation of background density for (a) 111, (b) 614, (c) I1I, (d) $614-$, (e) $614+$, along with ( f) respective buoyancy frequency profiles. Solid black line at the centre of (a–e) represents the disk (not to scale). Profiles are summarized in table 1 and also discussed in text. For $614-$, the upper half of the disk is in a constant-$N$ region and the lower is in the pycnocline. For $614+$, the disk is in the bottom constant-$N$ region with its upper edge $0.5D$ below the pycnocline.

For the I1I and 614 profiles, the disk centre coincides with the centre of the density profile and, furthermore, the disk lies entirely within the nonlinearly stratified region. Profiles $614-$ and $614+$ are obtained after shifting profile 614 (and not the disk) vertically by $-2.5D$ and $3D$, respectively. The disk centre is located at $z = 0$ for all five profiles and the negative and positive signs in the names $614-$ and $614+$ indicate whether the profile is shifted downwards or upwards, respectively. The vertical shift in $614-$ (figure 2d) is chosen so that the upper half of the disk lies in the linearly stratified region with $Fr = 6$ and the lower half in the pycnocline. The vertical shift in $614+$ (figure 2e) is chosen so that the disk lies entirely in the linearly stratified region with $Fr = 4$ while being very close to the pycnocline, specifically, its upper edge is $0.5D$ below the pycnocline.

The profile I1I is given by the following function, which was also used by Ermanyuk & Gavrilov (Reference Ermanyuk and Gavrilov2002, Reference Ermanyuk and Gavrilov2003) and Nicolaou et al. (Reference Nicolaou, Garman and Stevenson1995) in their studies:

(2.4)\begin{equation} \Delta \rho_{b}(z)={-}\frac{\rho_{0} \varepsilon}{2} \tanh{\frac{2z}{\delta}} . \end{equation}

The central nonlinear region of the 614 profile is obtained from the I1I profile. The linearly stratified regions of profile 614 are added by matching the slopes at $z=1.25\delta$ and $z=-\delta$ so that the density profile is continuously differentiable in $z$. This results in Froude numbers of $6.13$ and $3.74$ above and below the pycnocline layer respectively, for profile 614.

The profiles $614-$ and $614+$ are shifted with respect to 614 and their central nonlinear region is as follows:

(2.5)\begin{equation} \Delta \rho_{b}(z)={-}\frac{\rho_{0} \varepsilon}{2} \tanh{\frac{2(z - s)}{\delta}} , \end{equation}

where $s$ denotes the shift and takes values $-2.5D$ and $3D$ for $614-$ and $614+$, respectively.

The local Froude number corresponding to this pycnocline region can be calculated using the local background buoyancy frequency

(2.6)\begin{equation} Fr(z) = \frac{U_{\infty}}{\sqrt{\varepsilon gD}} \sqrt{\frac{\delta}{D}} \left|\cosh{\frac{2(z-s)}{\delta}}\right| . \end{equation}

The minimum value of $Fr(z)$ occurs at $z=s$. Thus, the minimum Froude number for profiles 614 and I1I is

(2.7)\begin{equation} Fr_{min} = Fr(z=s=0) = \frac{U_{\infty}}{\sqrt{\varepsilon gD}} \sqrt{\frac{\delta}{D}} . \end{equation}

Equation (2.7) highlights two important non-dimensional parameters for a body moving through a pycnocline: $U_{\infty }/\sqrt {\varepsilon gD}$, which is the conventional Froude number defined using reduced gravity $\varepsilon g$, and $\delta / D$, which is the non-dimensional thickness of the pycnocline. In the simulations, $U_{\infty }/\sqrt {\varepsilon gD} = 1/\sqrt {2}$ and $\delta / D = 2$ for profiles 614 and I1I to obtain $Fr_{min} = 1$. Note that for 111, which has constant linear stratification throughout, $Fr_{min}=1$ still holds. Parameters for each case are given in table 1.

Table 1. Physical parameters of the simulated cases. For each case, ${\it Re} = 5000$ and $Pr = 1$. The domain with $L_{r} = 60$, $L_{\theta } = 2{\rm \pi}$, $L_{x}^{-} = 30$ and $L_{x+} = 102$ is discretized using $N_{r} = 479$, $N_{\theta } = 128$ and $N_{x} = 2176$ points.

2.3. Domain and grid

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. The number of grid points employed for discretizing the domain is $N_{x}=2176$, $N_{\theta }=128$ and $N_{r}=479$ in the streamwise, azimuthal and radial directions, respectively, resulting in approximately $130$ million grid points. The disk is resolved into 47 312 triangles. The LES grid is non-uniform in the streamwise and radial directions and is designed to have high resolution. In terms of the Kolmogorov length ($\eta =(\nu ^{3} / \epsilon _{k}^{T})^{1/4}$), the maximum value of $\Delta x / \eta$ is $4.39$ and that of $\Delta r / \eta$ is $5.03$. Also, the turbulent dissipation rate used for the calculation of $\eta$ includes the resolved-scale dissipation ($\epsilon = 2\nu \langle s_{ij}' s_{ij}' \rangle$) as well as the subgrid-scale dissipation ($\epsilon _{sgs} = - \langle \tau _{sgs_{ij}}' s_{ij}' \rangle$ with $\tau _{sgs_{ij}}= -2 \nu _{sgs} S_{ij}$), where $s_{ij}'$ and $S_{ij}$ are the fluctuating and mean strain rates, respectively. See § 2.4 for averaging details. Most of the turbulent dissipation resides in the resolved scales.

2.4. Statistics

Each of the simulations is approximately $260$ non-dimensional time units ($D / U_{\infty }$) long, which is approximately $2.5$ flow throughs. Statistics are collected after the initial transience has subsided and statistical steady state is reached (which takes approximately 120 non-dimensional time units). Averaging is done over a total of $140$ non-dimensional time units after statistical steady state is reached. For any dependent variable in the simulation, $\langle {\cdot } \rangle$ is used to represent the time average and $'$ is used to represent the fluctuation about that average, over the multiple ensembles collected during the statistical steady state.

3.1. Steady lee waves

Disturbances in stratified environments generate internal waves. These can be steady lee waves generated by the body or unsteady internal waves generated by turbulence in the wake of body. In this section, the characteristics of the body generated lee waves on the vertical centre plane for 111 and 614 will be analysed using linear asymptotic theory. Note that since I1I is essentially unstratified away from the pycnocline layer, it does not show any lee waves.

Instantaneous contours of vertical velocity ($w$) for 111 (top half) and 614 (top and bottom halves) on a vertical plane passing through the centreline are plotted in figure 3(a–c). The amplitude of steady lee waves decays moving away from the body. The lee waves observed in 111 are symmetric with respect to $z =0$ while 614 has two different sets of lee waves (top and bottom) owing to the two different local values of $Fr(z)$ ($6.13$ above and $3.74$ below) in the linearly stratified regions above and below the pycnocline layer. The waves in the 614 case have a much smaller amplitude than in the 111 case and a larger wavelength. (Note that we use $w$ instead of $\langle w \rangle$ to visualize the steady lee waves to give a complete picture of the flow field. Unsteady internal waves, that are discussed in § 4.3 can be seen in figure 3(a–c) as distortions in the wave field near the wake that grow as $x/D$ increases.)

Figure 3. Instantaneous vertical velocity at $T = 200 D/U_{\infty }$ showing internal waves (a–c) and their amplitudes (d–f) along the dashed arrow labelled by $d/D$. (a,d) – 111 ($90^{\circ }$ plane); (b,e) – 614 ($90^{\circ }$ plane); (c, f) – 614 ($270^{\circ }$ plane). The instantaneous velocity contains both the steady lee wave, which is dominant in the far field, and the unsteady wake generated waves.

Linear theory, given by Voisin (Reference Voisin1991, Reference Voisin2003, Reference Voisin2007) is used to analytically predict the wavelength and amplitude of the lee waves. The theory computes a wave function, $\chi$ from a linearized set of inviscid equations involving a source term $q$ on the right-hand side of the continuity equation to model the moving body using a Green's function approach for large times ($Nt\gg 1$), which in our case translates to $x/D \gg Fr$. Once $\chi$ is known, $w$ is calculated. The expression for $w$ for the case of a Rankine ovoid as derived by Ortiz-Tarin et al. (Reference Ortiz-Tarin, Chongsiripinyo and Sarkar2019) and employed for a $4:1$ spheroid wake can be reduced to give the wave pattern on the vertical centreline plane as

(3.1)\begin{equation} w(x,y=0,z) \sim{-}\frac{mN}{{\rm \pi} U_{\infty} r_{xz}} \cos{\psi} \sin{\left(\frac{Na}{U_{\infty}} \cos{\psi}\right)} \sin{\left(\frac{N}{U_{\infty}} r_{xz}\right)} , \end{equation}

where $r_{xz}=\sqrt {x^{2}+z^{2}}$ and $\psi = \arctan (z/x)$. Also, $m$ and $a$ are calculated from the potential flow stagnation point solution for a Rankine ovoid of length $L_{ro}$ and cross-sectional diameter $D_{ro}$

(3.2a,b)\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}

Note that Ortiz-Tarin et al. (Reference Ortiz-Tarin, Chongsiripinyo and Sarkar2019) calculated $m$ and $a$ using the potential flow solution of a Rankine oval, which is a 2-D solution. Equation (3.2a,b) is based on the solution for a Rankine ovoid (which is three-dimensional), and serves as a correction to Ortiz-Tarin et al. (Reference Ortiz-Tarin, Chongsiripinyo and Sarkar2019). Note that, as per (3.1), the wave amplitude decreases with decreasing $N$ or increasing $Fr$.

The Rankine ovoid was chosen since it is close to the observed shape of the disk plus separation bubble and it has a simple potential flow solution – a 3-D source and a 3-D sink – needed for the linear theory of wave generation. The ovoid dimensions are based on the observed dimensions of the separation bubble. To apply linear theory to a bluff body such as a disk, the separation bubble created by the disk should also be included as a part of the extended body involved in wave generation. The complete extended body for 111, I1I and 614 cases resembles an ovoid with $L_{ro} \approx 2.5D$ and $D_{ro} \approx 1.5D$. It should be noted that the exact body shape is insignificant as far as the wavelength of the waves is considered, but the wave amplitude depends on the body shape.

Figure 3(d) compares the wave amplitude between simulation and the theoretical estimate of (3.1). The comparison is along the dashed line with an arrow in figure 3(a), which is at an angle of $45^{\circ }$ to the $x$-axis. The theory, when used with the Rankine ovoid approximation for the disk and the separated flow behind it, is able to accurately capture the amplitude variation specially moving away from the body for the 111 case which has linear stratification with constant $N$.

A similar analysis is performed for the lee waves of 614 after taking the Froude number in the analysis to be that of the linearly stratified region, i.e. $Fr = 6.13$ for figure 3(e) and $Fr = 3.74$ for figure 3( f). This procedure results in the correct prediction of the wavelength of the propagating lee waves. Theory is able to capture the order of magnitude of the significantly reduced (relative to 111) wave amplitude but the theoretical estimate of the wave amplitude is an under-prediction. The disk moves through a local stratification corresponding to $Fr = 1$ but the lee waves form and propagate in the linearly stratified regions where the Froude number is larger ($Fr = 6.13$ and $Fr = 3.74$). Therefore, the true wave amplitude corresponds to $Fr$ somewhat lower than that used in the theory and, thus, larger than the theoretical estimate. Since the wavelength is comparable to the variability scale of $N(z)$, Wentzel–Kramers–Brillouin (WKB) theory cannot be used and we do not proceed further with linear analysis.

3.2. Mean defect velocity

All cases show a strong effect of buoyancy on the wake since $Fr = 1$ at the disk is common to them but there are differences among the cases as elaborated below. Figure 4(a) shows the mean centreline defect velocity for the three cases plotted against the streamwise distance. The 111 case shows strong oscillatory modulation of the wake and its defect velocity, similar to the sphere wake (Pal et al. Reference Pal, Sarkar, Posa and Balaras2017). In contrast to 111, the modulation in 614 and I1I, although present, is small. The wavelength of modulation for 111 is $2 {\rm \pi}D Fr$ as in the linear theory result (3.1), and as noted in previous work.

Figure 4. (a) Mean defect velocity at centreline. (b) Mean streamwise velocity contours on vertical plane passing through centreline. Dashed yellow line represents the separation bubble ($U=0$).

All cases show a similar decay law of $U_0 \propto x^{-0.18}$ in the NEQ regime, which was also observed by Chongsiripinyo & Sarkar (Reference Chongsiripinyo and Sarkar2020) for a disk at a higher Reynolds number of $5 \times 10^{4}$. Around $x/D = 45 (Nt = 45)$, the decay law transitions to a faster decay rate, signalling the beginning of the Q2-D regime. A notable difference in the defect velocity profiles occurs at $x/D \sim 4$, where the initial dip in the profile for 111 is larger in magnitude and occurs earlier than that of 614 and I1I, owing to the stronger lee wave in 111 as compared with 614 and I1I. The mean velocity contours on a vertical streamwise cut (figure 4b) show that the lee wave modulation in the wake for 111 is stronger and its separation bubble is smaller, which is consistent with the earlier dip in the profile of defect velocity (figure 4a).

3.3. Turbulent kinetic energy

The nonlinearly stratified cases exhibit stronger levels of turbulent kinetic energy (TKE). This difference is illustrated by the contours of TKE at $x/D = 20$ in figure 5(a), plotted for the three stratification profiles. The pycnocline cases have a higher TKE (by 50 %–100 %) relative to the 111 case over the entire wake length as shown by the area-averaged values of TKE plotted in figure 5(b). The area averaging at any streamwise location is done in the region covered by a circle of radius $3D$, therefore containing most of the turbulent zone. The TKE is substantially larger in I1I and 614 relative to 111, especially in the NEQ regime. For example, I1I has almost twice the TKE of 111 at $x/D = 45$. It is worth noting the slight increase in TKE for the case 111 at $x/D \approx 50$ which is where the decay law for the mean defect velocity changes (§ 3.2), again pointing towards the transition to Q2-D regime.

Figure 5. (a) The TKE contours at $x/D = 20$. (b) Area-averaged values of TKE as a function of streamwise distance.

To explain the relatively high turbulence, we inspect the terms in the TKE budget, especially the lateral production (the production by vertical Reynolds shear stress is small) and the wave flux. Although the hyperbolic tangent profile is close to a linear profile near the centreline, the weaker stratification regions above and below are not as effective in suppressing near-wake turbulence stresses, which results in higher TKE production further downstream. Also, the weaker far-field stratification in 614 and I1I relative to 111 does not allow the full frequency range of waves generated in the $Fr = 1$ central region to escape into the far field. The reduction of wave energy flux traps fluctuations in the wake and results in a TKE increase. The two quantities, lateral TKE production and wave flux, are diagnosed in the following subsections. Although both contribute, it is found that the increase in TKE production is somewhat larger than the decrease in wave flux.

3.3.1. The TKE production

Contours of lateral production $P_{xy}$, which is the dominant component of TKE production for all three cases, are plotted in figure 6(a) at $x/D = 20$. The increased $P_{xy}$ in 614 and I1I is partly a result of the higher lateral Reynolds shear stress $\langle -u'v'\rangle$.

Figure 6. (a) Lateral production contours at ${x}/{D} = 20$. (b) Area-averaged values of lateral production as a function of streamwise distance.

Area-averaged values of the lateral production over a circle of radius $3D$ are shown in figure 6(b). The lateral production is consistently larger in the 614 and I1I cases over $10 < x/D < 50$ and is as much as twice at some locations.

3.3.2. Energy flux of wake generated unsteady waves

The dispersion relation for internal gravity waves in a stratified non-rotating medium

(3.3)\begin{equation} \varOmega = N \cos {\varTheta} , \end{equation}

where $\varTheta$ is the angle of the phase line with the vertical, limits the maximum allowable frequency to $N$. Since weaker stratification supports a narrower frequency band of wave propagation, 614 and I1I can be expected to allow a smaller wave energy flux to radiate through the pycnocline layer as compared with 111. This is indeed the case as can be seen by the contours of the radial wave flux $\langle p' u_{r} ' \rangle$ in figure 7(a), where it is evident that the wave flux is restricted within the pycnocline layer ($-2 < z/D < 2$) instead of being radiated away as in 111. Note that the means of the radial velocity and pressure are subtracted when computing the flux so as to discard the component from the steady lee waves.

Figure 7. (a) Radial waveflux contours at ${x}/{D} = 20$. (b) Total waveflux integrated over rectangular perimeter outside the turbulent wake as a function of streamwise distance.

A quantitative comparison among the cases is performed by computing the line integral of the internal wave flux over a rectangular perimeter. Specifically, the integral

(3.4)\begin{equation} \oint_{C} \frac{1}{\rho_{0}} \langle {\boldsymbol u}' p' \rangle \boldsymbol{\cdot} {\boldsymbol n} \, {\rm d}l , \end{equation}

which is a function of $x/D$, is calculated using a rectangle $C$ with vertical sides at $y/D = \pm 4$ and top/bottom sides at $z/D = \pm 2.5$ (thereby containing the pycnocline layer). Figure 7(b) shows that the flux of wave energy radiated away from the wake of 111 is much larger than 614 and I1I in the early part of the NEQ regime. However, at $x/D = 40$ in the late NEQ stage, the fluxes for 614 and I1I overtake the flux of 111. Both of these trends can be explained by the schematic shown in figure 8(a). The wake generated internal waves that are trapped inside undergo complex wave–wave interaction after getting vertically restricted, resulting in sideward escape of the wave flux. The sideward escape can be verified by using an appropriate scaling argument for the vertical as well as lateral wave flux at the top and side boundaries of the rectangle, respectively,

(3.5a,b)\begin{equation} \langle p'w' \rangle \sim p_{rms} w_{rms},\quad \langle p'v' \rangle \sim p_{rms} v_{rms}. \end{equation}

Here, the subscript rms stands for the root mean square value. The values of $w_{rms}$ at the point $P_{top}$, $v_{rms}$ at the point $P_{side}$ and $p_{rms}$ at $P_{top}$ and $P_{side}$ are shown in figure 8(b–e), respectively. Comparison of figures 8(d) and 8(e) shows that, although pressure fluctuations at the top boundary are highest initially for 111 among all cases and boundaries, the pressure fluctuations for 614 and I1I increase with time at the side boundaries because of the wave trapping and the resulting sideward escape. The velocity fluctuations at the boundaries also show similar trends. The sideward escape of the wave energy takes some time to manifest which, in the simulation frame, translates to streamwise distance, hence causing the slight delay in the increase of integrated wave flux for 614 and I1I. Figure 8(b–e) is consistent with the following result: most of the wave flux for 614 and I1I that escapes the wake core does so in the form of lateral wave flux $\langle p'v' \rangle$ unlike 111 where the vertical wave flux $\langle p'w' \rangle$ is the major contributor.

Figure 8. (a) Schematic showing trapping of internal waves for 614 and I1I compared with 111, (b) $w_{rms}$ at $P_{top}$, (c) $v_{rms}$ at $P_{side}$, (d) $p_{rms}$ at $P_{top}$, (e) $p_{rms}$ at $P_{side}$.

4.1. Taylor–Goldstein equation and Kelvin wake waves

In cases involving pycnoclines, there is a family of steady waves which resemble a Kelvin ship wave pattern, on horizontal planes inside the pycnocline layer. Intuitively, a stratification profile with a large central value of buoyancy frequency ($Fr \leq O(1)$) that changes rapidly in the vertical (over $\delta /D \leq O(1)$) to a small value at the edges can be thought of as a smoothed density jump and thus one can expect a pattern resembling that generated on the air–water interface by a moving ship. Both centred and shifted stratification profiles show Kelvin wake wave patterns. However, as demonstrated here, the wake wave structure differs qualitatively from the centred 614 wake when the stratification profile is shifted with respect to the disk centre, as in $614-$ and $614+$ .

The dispersion relation of air–water surface gravity waves leads to the waveforms of air–water Kelvin ship waves. To obtain an appropriate dispersion relation for the nonlinear continuously stratified profiles, the Taylor–Goldstein equation is numerically solved

(4.1)\begin{equation} \frac{{\rm d}^{2} \phi (z)}{{\rm d} z^{2}} + k^{2} \left(\frac{N^{2}(z)}{\varOmega^{2}} - 1\right)\phi (z) = 0. \end{equation}

Here, $\phi (z)$ is the vertical displacement eigenfunction, and $k$ and $\varOmega (k)$ are the wavenumber and angular frequency corresponding to the dispersion relation. The Taylor–Goldstein equation is numerically solved as an eigenvalue problem for $\phi (z)$ by treating $\varOmega$ as an eigenvalue and $k$ as a parameter, e.g. Robey (Reference Robey1997). The boundary conditions used are

(4.2)\begin{equation} \phi(z={-}L_{r}) = \phi(z=L_{r}) = 0. \end{equation}

Also, based on Sturm–Liouville theory, the eigenfunctions are normalized using

(4.3)\begin{equation} \int_{z={-}L_{r}}^{z=L_{r}} \phi_{m} N^{2}(z) \phi_{n} \,{\rm d}z = \delta_{mn}. \end{equation}

Figures 9(a) and 9(b), respectively, show the mode 1 and mode 2 eigenfunctions for different values of $k$. The numerically calculated dispersion relation is also compared in figure 9(c) with the approximation given by Barber (Reference Barber1993)

(4.4)\begin{equation} \varOmega(k) = \frac{c_{p_{0}}k}{1 + \dfrac{c_{p_{0}}k}{N_{max}}}. \end{equation}

Here, $c_{p_{0}}$ is the limiting long-wave phase speed of a given mode at $k=0$, which is obtained by extrapolation, i.e. $c_{p_{0}}$ is the slope at the origin of the mode-specific curve in figure 9(c). For the chosen profiles, $c_{p_{0}}/U_{\infty }$ for mode 1 and mode 2 are $2.23$ and $0.68$, respectively. Since (4.4) provides a very good approximation to the dispersion relation, it will be used for the remaining analysis of this section.

Figure 9. (a) Mode 1 and (b) mode 2 eigenfunctions for different values of wavenumber in (4.1). (c) Numerically calculated dispersion relation compared with approximation given by Barber (Reference Barber1993). (d) Mode 1 and (e) mode 2 waveforms as calculated using (4.5a,b).

Using the dispersion relation, phase and group velocities associated with each mode can be calculated using $c_{p} = \varOmega / k$ and $c_{g} = {\rm d}\varOmega / {\rm d}k$, and then substituted in the expressions given by Keller & Munk (Reference Keller and Munk1970) to obtain the modal wave patterns corresponding to each eigenfunction

(4.5a,b)\begin{equation} x = \dfrac{\varphi U_{\infty} \left(1 - \dfrac{c_{p} c_{g}}{U_{\infty}^{2}}\right)}{k(c_{p}-c_{g})},\quad y = \dfrac{\varphi c_{g} \Bigg(1 - \dfrac{c_{p}^{2}}{U_{\infty}^{2}}\Bigg)^{{1}/{2}}}{k(c_{p}-c_{g})} , \end{equation}

where $(x,y)$ corresponds to a locus of points on the horizontal plane parametrically given in terms of $k$ for a constant value of phase, $\varphi$. Figures 9(d) and 9(e) respectively show mode 1 and mode 2 wave patterns, plotted using (4.5a,b).

Figure 10 shows the instantaneous radial velocity contours for 111, 614, $614-$ and $614+$ plotted on a half-horizontal plane passing through the respective pycnocline centre. (For 111, the plot is on $\theta = 0^{\circ }$ plane.) The waves in 614 have constant-phase lines which resemble mode 2 (shown earlier in figure 9e) but not mode 1. In contrast, the phase lines for $614-$ and $614+$ are more complex. In the region $0 < y/D < 10$ they resemble the mode 1 pattern of figure 9(d) but, for $y/D > 10$, they resemble the mode 2 pattern. The manifestation of any mode by a disturbance in the pycnocline layer depends on the location of the disturbance with respect to the pycnocline layer. In 614, the disk centre travels along the centre of the pycnocline layer to displace fluid symmetrically in the upward and downward directions, corresponding to the antisymmetric mode 2 eigenfunction for displacement and therefore generates mode 2 waves. Any vertical offset from the pycnocline centre modifies the antisymmetric pattern of the displacement so as to also involve some contribution from the mode 1 waveform. The presence of a mode 2 wave pattern in the central region of the lateral plane and a mode 1 pattern otherwise is in agreement with the visualizations in the experiment of Robey (Reference Robey1997). For profile $614+$, the contours look similar to that of $614-$ but with slightly lower intensity because of the disk being further away from the centre of the pycnocline. (Note that the wave pattern for I1I is the same as that of 614 because they have the same hyperbolic tangent function modelling the nonlinearity.)

Figure 10. Instantaneous contours of radial velocity on a half-pycnocline-centre plane at $T = 200 D/U_{\infty }$: (a) 111, (b) $614-$, (c) 614 and (d) $614+$. Dashed lines correspond to mode 2 waves which were shown in figure 9(e).

4.2. Distinction between lee waves and Kelvin wake waves

Kelvin wake waves are steady in the frame of reference of the moving body like the lee waves but constitute a new family of waves in terms of their distinct features and appear only when the stratification is non-uniform. Figure 11 shows instantaneous contours of radial velocity on horizontal planes for the centred profiles 111, 614 and I1I at different vertical locations. At $z/D=-2.5$ (figure 11a–c), the contours have superficial resemblance but they represent fundamentally different waves. For 111, since there is no nonlinear gradient, Kelvin wake waves are not formed on the horizontal plane passing through the domain centreline at $z/D=0$, (e.g. the previously shown figure 10a), and so figure 11(a) has the horizontal imprint of only lee waves. However, for 614 and I1I, the contours at $z/D=-2.5$ in figure 11(b,c) are mostly the imprints of the Kelvin wake waves in figure 10(c). For 614, Kelvin wake waves disappear farther away from the region of nonlinear stratification and the contours transition to those of lee waves corresponding to $Fr=3.74$ (figure 11e,h,k). For I1I, once the Kelvin wake waves disappear, lee waves are not observed since the profile is unstratified outside the pycnocline region (figure 11f,i,l). Turning back to case 111, we see the lee waves weakening as $z/D$ is increased (figure 11d,g,j) but their pattern is still clear at $z/D = -30$.

Figure 11. Instantaneous contours of radial velocity for the centred cases on vertically offset horizontal planes at $T = 200 D/U_{\infty }$ for 111 (a,d,g,j), 614 (b,e,h,k) and I1I (c, f,i,l). The vertical offset of each plane increases from (a–c) to (j–l).

4.3. Wake generated internal gravity waves

The internal gravity waves generated by the turbulent wake are unsteady. These waves are visualized in figure 12 for case 111 by contours of $\partial w' / \partial z$ on the $\theta = 90^{\circ }$ vertical plane. The phase lines of the radiated waves in the far field are seen to cluster around a characteristic inclination angle, which suggests a narrow frequency band in the far field according to the dispersion relationship, (3.3), for the intrinsic (in a frame where the background has zero velocity) wave frequency.

Figure 12. Instantaneous contours of vertical derivative of vertical velocity fluctuation at $T = 200 D/U_{\infty }$ showing wake generated internal gravity waves on the $\theta = 90^{\circ }$ plane for 111. Black lines are inclined at an angle of $39^{\circ }$ with the vertical.

Figure 12 also shows solid black lines plotted at an angle of $39^{\circ }$ from the vertical, which also seems to be the angle for the waves ($39^{\circ } \pm 2^{\circ }$). This gives $\varOmega / N = 0.78 \pm 0.03$.

For profile $614+$, the $\partial w' / \partial z$ contours on the entire vertical plane ($\theta = 90^{\circ }$ and $270^{\circ }$) are plotted in figure 13. Phase lines with a dominant inclination angle are seen in the negative $z/D$ region where the waves with downward group velocity travel in an entirely linear stratification ($Fr = 3.74$). The inclination angles in this region of $z/D < 0$ cluster around $\varTheta = 44^{\circ }$, which gives $\varOmega / N \approx 0.72$ (note that $N$ here corresponds to the local $Fr = 3.74$). The upward moving waves generated by the disk enter a region of nonlinear stratification where there is significant small-scale variability. As the waves exit the pycnocline layer from the top, they have near-zero inclination with respect to the vertical, i.e. $\varTheta \approx 0^{\circ }$, which implies that only the high-frequency near-$N$ ($N$ here again corresponds to the local $Fr = 6.13$) part of the wake generated waves escapes into the upper region of weak stratification. Thus, a significant portion of the wake generated waves is trapped within the pycnocline leading to the small-scale variability in that region.

Figure 13. Instantaneous contours of vertical derivative of vertical velocity fluctuation at $T = 200 D/U_{\infty }$ showing wake generated internal gravity waves on the $\theta = 90^{\circ }$ and $180^{\circ }$ planes for $614+$.

4.4. Wake characteristics

Wake characteristics, namely, the mean defect velocity and TKE, show qualitatively different behaviour for the shifted pycnocline cases $614-$ and $614+$ relative to the centred 614 profile. Figure 14 compares the streamwise evolution of centreline mean defect velocity ($U_0/U_\infty$) among all five simulated cases. Recall that, for $614+$, the disk is entirely in the $Fr = 3.74$ region with its upper edge in the pycnocline and, for $614-$, half of the disk is in the $Fr = 6.13$ region and half in the pycnocline. Since, for both $614-$ and $614+$, the disk wake evolves in relatively weaker stratification relative to the $Fr =1$ stratification seen by the disk in 614, there is a faster rate of decay of $U_0$ ($\propto x^{-0.35}$) in these cases (figure 14). $614-$ has a weaker effective stratification relative to $614+$ and, therefore, has a somewhat lower $U_0$. It is also worth mentioning that, since the simulation domain is not long enough in the streamwise direction, the limited data prevent discerning of the $x^{-0.76}$ decay law in the late Q2-D regime, which is otherwise observed in Spedding (Reference Spedding1997) and can also be deduced from the drag coefficient as shown by Meunier & Spedding (Reference Meunier and Spedding2004).

Figure 14. Evolution of centreline mean defect velocity for the five simulated cases.

The TKE contours plotted in figure 15 also show that the buoyancy effect is weaker for the shifted profiles; $614-$ and $614+$ show far higher TKE than the other three cases that have $Fr = 1$ at the disk centre. Since the bottom half of the disk of $614-$ is inside the pycnocline layer, the wake core is not symmetric about the horizontal axis and appears to be somewhat compressed from the bottom where the effective stratification is higher. The TKE profile at $x/D = 20$ is more symmetric for $614+$.

Figure 15. The TKE contours at $x/D = 20$ for all five cases; (a) 111, (b) 614, (c) I1I, (d) 614−, (e) 614+.