Disk wakes in nonlinear stratification

The wake of a disk is simulated by solving the three dimensional, incompressible, unsteady form of the conservation equations for mass, momentum and density. A high-resolution large eddy simulation (LES) with the Boussinesq approximation for density effects is used. The disk 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}$ (Fig. 1). The density field $(\rho(\mathbf{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}(\mathbf{x},t))$ so that $\rho(\mathbf{x},t)=\rho_{0}+\Delta\rho_{b}(z)+\rho_{d}(\mathbf{x},t)$. The Reynolds number of the flow, defined as $\Rey=U_{\infty}D/\nu$ with $D$ being the diameter of the disk, 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 behavior is described in section 2.2. The minimum value of the Froude number ($Fr_{min}$) is set to $1$ for all the cases.

The filtered nondimensional equations are as follows

where $u_{i}$ refers to the filtered velocities in $x$, $y$, and $z$ directions for $i=1,2$, and $3$ respectively. $\nu_{sgs}$ and $\nu$ in Eq. (2) are the subgrid-scale kinematic viscosity and the kinematic viscosity, respectively, while $\kappa_{sgs}$ and $\kappa$ in Eq. (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 et al. (2010) who varied $Pr$ between 0.2 and 7. The dynamic eddy viscosity model (Germano et al., 1991) is employed to obtain $\nu_{sgs}$, following the implementation of Chongsiripinyo & Sarkar (2020), and the subgrid Prandtl number is also set to unity to obtain $\kappa_{sgs}$. The parameters used for nondimensionalization 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 discretised 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, 1999). The disk is represented by the immersed body method of Balaras (2004) and Yang & Balaras (2006). 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, 1976). 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.

The profiles chosen for the variation of background density in the five cases are shown in Fig. 2. Except for the benchmark 111 case, the profiles have nonuniform $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 $\mbox{{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 (Fig. 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.

For the I1I and 614 profiles, the disk center coincides with the center 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 center 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-$ (Fig. 2d) is chosen so that the upper half of the disk lies in the linearly stratified region with $\mbox{{Fr}}=6$ and the lower half in the pycnocline. The vertical shift in $614+$ (Fig. 2e) is chosen so that the disk lies entirely in the linearly stratified region with $\mbox{{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 (2002, 2003) and Nicolaou et al. (1995) in their studies:

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:

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:

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

Eq. 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 nondimensional 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.

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 LES grid is non uniform in the streamwise and radial direction 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}^{\prime}s_{ij}^{\prime}\rangle$) as well as the subgrid scale dissipation ($\epsilon_{sgs}=-\langle\tau_{sgs_{ij}}^{\prime}s_{ij}^{\prime}\rangle$ with $\tau_{sgs_{ij}}=-2\nu_{sgs}S_{ij}$), where $s_{ij}^{\prime}$ and $S_{ij}$ are the fluctuating and mean strain rates respectively. Most of the turbulent dissipation resides in the resolved scales.

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 center plane for 111 and 614 will be analyzed 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 centerline are plotted in Fig. 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.

Linear theory, given by Voisin (1991, 2003, 2007) 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>>1$), which in our case translates to $x/D>>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. (2019) and employed for a 4:1 spheroid wake can be reduced to give the wave pattern on the vertical centerline plane as:

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}$:

Note that Ortiz-Tarin et al. (2019) calculated $m$ and $a$ using the potential flow solution of a Rankine oval, which is a 2-dimensional solution. Eq. (9) is based on the solution for a Rankine ovoid (which is 3-dimensional), and serves as a correction to Ortiz-Tarin et al. (2019). Note that, as per Eq. 8, the wave amplitude decreases with decreasing $N$ or increasing Fr.

The justification for using the expression for a Rankine ovoid instead of a disk is as follows. 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.

Fig. 3(d) contrasts the wave amplitude between simulation and the theoretical estimate of Eq. 8. The comparison is along the dashed line with an arrow in Fig. 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 Fig. 3(e) and $Fr=3.74$ for Fig. 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 $\mbox{{Fr}}=1$ but the lee waves form and propagate in the linearly stratified regions where the Froude number is larger ($\mbox{{Fr}}=6.13$ and $\mbox{{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)$, WKB theory cannot be used and we do not proceed further with linear analysis.

All cases show a strong effect of buoyancy on the wake since $\mbox{{Fr}}=1$ at the disk is common to them but there are differences among the cases as elaborated below. Fig. 4(a) shows the mean centerline 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., 2017). In contrast to 111, the modulation in 614 and I1I, although present, is small. The wavelength of modulation for 111 is $2\pi DFr$ as in the linear-theory result (Eq. 3.1), and as noted in previous work.

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 (2020) for a disk at a higher Reynolds number of $5\times 10^{4}$. The decay law transitions to $x^{-0.45}$ around $x/D=50(Nt=50)$ signaling the beginning of the Q2D 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 to 614 and I1I. The mean velocity contours on a vertical streamwise cut (Fig. 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 (Fig. 4a).

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 Fig. 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 Fig. 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 (Section 3.2), again pointing towards the transition to Q2D regime.

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 centerline, 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.

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 large central value of buoyancy frequency ($\mbox{{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 centered and shifted stratification profiles show Kelvin wake wave patterns. However, as demonstrated here, the wake wave structure differs qualitatively from the centered 614 wake when the stratification profile is shifted with respect to the disk center 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:

Here, $\phi(z)$ is the vertical displacement eigenfunction, and $k$ and $\Omega(k)$ are the wavenumber and angular frequency corresponding to the dispersion relation. The TGE is numerically solved as an eigenvalue problem for $\phi(z)$ by treating $\Omega$ as an eigenvalue and $k$ as a parameter, eg. by Robey (1997). The boundary conditions used are

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

Fig. 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 Fig. 9(c) to the approximation given by Barber (1993):

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 Fig. 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 Eq. 16 provides a very good approximation to the dispersion relation, it will be used for the remaining analysis of this section.

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

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$. Fig. 9(d) and 9(e) respectively show mode 1 and mode 2 wave patterns, respectively, plotted using Eq. 17.

Fig. 10 shows the instantaneous radial velocity contours for 111, 614, $614-$ and $614+$ plotted on a half-horizontal plane passing through the respective pycnocline center. (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 Fig. 9e) but not mode 1. In contrast, the phase lines for $614-$ and $614+$ are more complex. In the region $0<z<10$ they resemble the mode 1 pattern of Fig. 9(d) but, for $z>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 center travels along the center of the pycnocline layer to displace fluid symmetrically in the upward and downward direction, corresponding to the antisymmetric mode 2 eigenfunction for displacement and therefore generates mode 2 waves. Any vertical offset from the pycnocline center 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 (1997). 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 center of the pycnocline. (Note that the wavepattern for I1I is the same as that of 614 because they have the same hyperbolic tangent function modeling the nonlinearity.)

Kelvin wake waves are steady in the frame of reference of the moving body like the lee waves but constitute a distinct family of waves and appear only when the stratification is nonuniform. Fig. 11 shows instantaneous contours of radial velocity on horizontal planes for the centered profiles 111, 614, and I1I at different vertical locations. At $z/D=-2.5$ (Fig. 11 (a-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 centerline at $z/D=0$, (e.g. the previously shown Fig. 10a), and so Fig. 11(a) has the horizontal imprint of only lee waves. However, for 614 and I1I, the contours at $z/D=-2.5$ in Fig. 11(b,c) are mostly the imprints of the Kelvin wake waves in Fig. 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 $\mbox{{Fr}}=3.74$ (Fig. 11(e,h,k)). For I1I, once the Kelvin wake waves disappear, lee waves are not observed since the profile is unstratified outside the pycnocline region (Fig. 11(f,i,l)). Turning back to case 111, we see the lee waves weakening as $z/D$ is increased (Fig. 11(d,g,j)) but their pattern is still clear at $z/D=-30$.

The internal gravity waves generated by the turbulent wake are unsteady. These waves are visualized in Fig. 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 a characteristic inclination angle which suggests a narrow frequency band in the far field according to the dispersion relationship, Eq. 10, for the intrinsic (in a frame where the background has zero velocity) wave frequency. Fig. 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 $\Omega/N=0.78\pm 0.03$.

For profile $614+$, the $\partial w^{\prime}/\partial z$ contours on the entire vertical plane ($\theta=90^{\circ}$ and $270^{\circ}$) are plotted in Fig. 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 ($\mbox{{Fr}}=3.74$). The inclination angles in this region of $z/D<0$ cluster around $\Theta=44^{\circ}$ which gives $\Omega/N\approx 0.72$ (Note that $N$ here corresponds to the local $\mbox{{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. $\Theta\approx 0^{\circ}$, which implies that only the high-frequency near-$N$ ($N$ here again corresponds to the local $\mbox{{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.

Wake characteristics, namely, the mean defect velocity and turbulent kinetic energy, show qualitatively different behavior for the shifted pycnocline cases $614-$ and $614+$ relative to the centered 614 profile. Fig. 14 compares the streamwise evolution of centerline mean defect velocity ($U_{0}/U_{\infty}$) among all five simulated cases. Recall that for $614+$, the disk is entirely in the $\mbox{{Fr}}=3.74$ region with its upper edge in the pycnocline and, for $614-$, half of the disk is in the $\mbox{{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 (Fig. 14). $614-$ has a weaker effective stratification relative to $614+$ and, therefore, has a somewhat lower $U_{0}$.

The TKE contours plotted in Fig. 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 $\mbox{{Fr}}=1$ at the disk center. 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+$.