2 Problem formulation, numerical details and validation

Motivated by the unanswered question regarding near-wake turbulence when $Fr$ is low but $Re$ is not, we use DNS to investigate the flow past a sphere at $Re=3700$ and $Fr\in \{0.025,1\}$ . The three-dimensional Navier–Stokes equations are solved in a cylindrical coordinate system on a staggered grid using an immersed boundary method (IBM) (Balaras Reference Balaras2004; Yang & Balaras Reference Yang and Balaras2006) for representing the sphere.

The simulation parameters, domain size and grid distribution for the different cases are given in table 1. High resolution is used at the sphere surface (20 points across the boundary layer thickness at the point of maximum wall shear stress) and in the wake. The radial grid spacing is $\unicode[STIX]{x0394}r\simeq 0.0016$ in the cylindrical region ( $r<0.65$ ) that encloses the sphere, the azimuthal direction has 128 points, and $\unicode[STIX]{x0394}x\simeq 0.0016$ near the surface. The grid has mild stretching, radially and streamwise, away from the body. The IBM results and the grid resolution to resolve the flow have been successfully validated in the unstratified case against both previous simulations and laboratory experiments. Figure 1(a,b) shows that the variations of the surface pressure coefficient, $C_{p}$ , and the surface shear stress, $(\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D70C}U^{2})Re^{0.5}$ ( $\unicode[STIX]{x1D70F}$ is the shear stress and $\unicode[STIX]{x1D70C}$ is the reference density), as a function of the azimuthal angle, match well with results in the available literature. Table 2 shows that key characteristics of the near-body flow such as the Strouhal number ( $St$ ) of the dominant shedding frequency, the separation angle ( $\unicode[STIX]{x1D711}_{s}$ ), the coefficient of drag ( $\overline{C_{d}}$ ) and the pressure coefficient $(\overline{C_{pb}})$ at the rearward stagnation point also match with previously reported values.

Figure 1.

Validation of the unstratified wake: (a) pressure coefficient, $C_{p}$ , (b) normalized drag coefficient, $(\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D70C}U^{2})Re^{0.5}$ . Here, $\unicode[STIX]{x1D703}$ is the azimuthal angle, with $\unicode[STIX]{x1D703}=0$ corresponding to the forward stagnation point.

Table 1.

Simulation parameters. The sphere is located at $(0,0,0)$ . The substantial domain size in the radial and upstream directions, along with the sponge region, eliminates the spurious reflection of internal waves.

3 Results and discussion

Figure 2 shows the downstream evolution of the turbulent kinetic energy (TKE) integrated over cross-stream ( $x_{2}{-}x_{3}$ ) planes for cases with different $Fr$ . It should be noted that $x_{3}$ denotes the vertical coordinate, the horizontal directions are $x_{1}$ (streamwise) and $x_{2}$ (lateral), and the sphere centre is at the origin. All statistics are computed after the initial transient by time averaging over an interval of $1.5L_{x}/U$ which is sufficient to obtain converged statistics. Buoyancy in a stratified wake has been found to suppress turbulence in previous studies, and, accordingly, the TKE decreases when $Fr$ decreases from 1 to 0.8 to 0.5. However, the trend reverses when $Fr$ decreases to 0.25 and beyond: the TKE increases with decreasing $Fr$ . The value of the TKE in the $Fr=0.25$ case increases to a level comparable to the $Fr=0.8$ case, and a further decrease of $Fr$ to 0.21 leads to values of the TKE larger than in the unstratified case. Subsequent reduction in $Fr$ beyond 0.21 leads to progressive augmentation of the TKE.

Figure 2.

Evolution of integrated TKE in the streamwise direction. The plotted quantity is the area-integrated TKE normalized using $U$ and $D$ .

Table 2.

Comparison of the different statistical flow features of the near-body flow in the present DNS with experimental measurements and numerical results available in the literature. Here, LES stands for large-eddy simulation, $St=fD/U$ is the non-dimensional vortex shedding frequency, $\unicode[STIX]{x1D711}_{s}$ is the azimuthal separation angle, $C_{d}$ is the drag coefficient and $C_{pb}$ is the rearward pressure coefficient at $\unicode[STIX]{x1D711}=\unicode[STIX]{x03C0}$ .

To understand the remarkable regeneration of fluctuations in the near wake at low $Fr$ , contour plots of the azimuthal vorticity magnitude in the horizontal $(x_{1}{-}x_{2})$ and vertical $(x_{1}{-}x_{3})$ planes (figure 3) are examined. The near-wake dynamics changes qualitatively for cases with $Fr\leqslant 0.25$ , as elaborated below. The $Fr=1$ wake displays the anisotropy of a moderate- $Fr$ wake: a large spread in the horizontal plane (figure 3 a) and small-scale structures associated with the shear layer instability, while, in the vertical plane, the separated boundary layers (figure 3 b) contract, followed by an undulation of the wake. At $Fr=0.5$ (not shown here), the recirculation bubble is steady, the disintegration of the shear layer is suppressed in the horizontal plane and the separating shear layers dip to the centreline in the vertical plane. The shear layer formed by the separating boundary layer exhibits large steady waviness in the vertical plane, there is little unsteadiness in the near wake and, therefore, the TKE for $Fr=0.5$ is insignificant, as was shown in figure 2. A quasisteady recirculation bubble attached to the sphere is found in the horizontal plane (figure 3 c) for a larger stratification, $Fr=0.25$ . At the end of the recirculation zone, the wake undergoes an unsteady undulation with the shedding of vortices further downstream. The shear layer in the vertical direction (figure 3 d) manifests waviness (induced by lee waves), but the instability does not break down into turbulence. The flow between the upper and lower shear layers displays thin strips of enhanced vorticity symptomatic of vorticity layering.

Figure 3.

Instantaneous azimuthal vorticity magnitude on the horizontal $x_{1}{-}x_{2}$ centre plane ( $x_{3}=0$ ) and the vertical $x_{1}{-}x_{3}$ centre plane ( $x_{2}=0$ ). Snapshots compared among cases with different $Fr$ . The plotted vorticity is normalized using  $U$ and  $D$ .

The flow organization changes significantly with further decrease in $Fr$ to 0.125 and beyond. There is unsteady motion of the shear layers in the horizontal plane accompanied by patches of small-scale turbulence (figure 3 e) as compared with the steady recirculation bubble in the $Fr=0.25$ wake. This reappearance of small-scale fluctuations at $Fr=0.125$ occurs due to unsteady vortex shedding in the horizontal plane, which results in both flapping and destabilization of the shear layer. A similar vertical layering of vorticity to that at $Fr=0.25$ is also seen at $Fr=0.125$ but, in this case, the layers roll up intermittently to form Kelvin–Helmholtz (KH) billows (figure 3 f) which then break down into finer-scale fluctuations. A secondary instability of pancake vortices in the far wake to form KH rolls was noted in previous temporal simulations (Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011) for sufficiently high $Re$ . In the present near wake, the perturbations provided by the horizontal flapping motion and the value of the local $Re$ are sufficient to destabilize the vertically layered vorticity into KH billows. As $Fr$ approaches 0.025, the unsteady vortex shedding from the sphere in the horizontal plane becomes more noticeable. The TKE in the region $x/D<1$ which belongs to the very near wake is also the largest among all simulated cases, as was shown in figure 2. In the horizontal plane (figure 3 g), there are coherent vortices with interspersed threads of rolled-up vorticity. In the vertical plane (figure 3 h), layered vortical structures are seen but do not manifest KH billows. The fact that KH billows are absent in the $Fr=0.025$ case will be explained, based on the value of the buoyancy Reynolds number and the scaling analysis of Riley & deBruynKops (Reference Riley and deBruynKops2003) and Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), later in the paper. The vorticity pattern at $Fr=0.025$ appears to have less fine-scale activity relative to $Fr=0.125$ . Internal gravity waves at the body can be seen in the vertical plane (figure 3 d,f,h), but their discussion is deferred to future work.

Both the mean and the turbulent kinetic energy are increasingly dominated by horizontal motions as $Fr$ decreases to 0.25 and below. The evolution of the ratio of the area-integrated mean kinetic energies (MKEs) of the horizontal component $(\text{MKE}_{11}+\text{MKE}_{22})$ and the vertical component $(\text{MKE}_{33})$ is shown in figure 4(a). For $Fr=1$ , the horizontal MKE is larger near the sphere, but, beyond $x_{1}/D\approx 5$ , the MKE becomes similarly distributed among the horizontal and vertical components. The undulations after $x_{1}/D\approx 5$ signify the exchange of MKE between the horizontal and vertical components. The ratio $(\text{MKE}_{11}+\text{MKE}_{22})/\text{MKE}_{33}$ for $Fr=0.25$ and 0.125 characterizes the transition of the near wake into quasihorizontal motion. The case with $Fr=0.025$ exhibits the complete dominance of horizontal motion, present primarily in the form of layered coherent vortices that span a wide lateral ( $x_{2}$ ) extent. The streamwise variation of the components of the TKE for $Fr=1$ and 0.125 is presented in figure 4(b). The components of the TKE for $Fr=1$ evolve in a similar manner, whereas for $Fr=0.25$ (not shown here) the streamwise $(\text{TKE}_{11})$ and spanwise $(\text{TKE}_{22})$ components are larger relative to the vertical $(\text{TKE}_{33})$ component. A significant difference between the horizontal $(\text{TKE}_{11},\text{TKE}_{22})$ and vertical components is observed as $Fr$ is further decreased to 0.125 (shown here) and 0.025 (not shown here).

Figure 4.

Evolution of (a) the ratio of area-integrated horizontal and vertical MKEs and (b) the components of the integrated TKE, in the streamwise direction. The area integration is over the $x_{2},x_{3}$ plane normal to the streamwise direction.

Temporal spectra are examined to quantify buoyancy effects on the frequency content of the lateral velocity,  $v$ . Figure 5(a) shows that there is a significant decrease of energy at all frequencies when the stratification increases to change $Fr$ from 1 to 0.25. However, a further decrease of $Fr$ to 0.125 and 0.025 shows a re-energization of fluctuations at all frequencies. There is a strong low-frequency peak in these cases: (i) $St=\unicode[STIX]{x1D714}D/U=0.163$ for $Fr=0.125$ , (ii) $St=0.200$ for $Fr=0.025$ . Secondary peaks of $E_{vv}$ at harmonics of the low-frequency mode are also evident. There is substantial energy, much larger than at $Fr=0.25$ , at the intermediate frequencies as well. It should be noticed that for flow over a circular cylinder in an unstratified environment at $Re=3900$ , the shedding frequency is found to be ${\approx}0.2$ (Parnaudeau et al. Reference Parnaudeau, Carlier, Heitz and Lamballais2008). Therefore, with increasing stratification, the vortex shedding of a sphere shifts towards that of a circular cylinder. This is because the flow at depths larger than $O(U/N)$ with respect to the top of the sphere tends to divert around the sphere rather than over the sphere because of the potential energy barrier. We emphasize that the low- $Fr$ near wake, apart from the similarity of vortex shedding, is quite different from the unstratified cylinder wake, where the strong inhibition of vertical fluctuations by buoyancy is absent. For example, the vertical velocity spectra $E_{ww}$ (figure 5 b) at $Fr=0.125$ and $Fr=0.025$ have much smaller amplitudes relative to their corresponding horizontal counterparts, $E_{vv}$ , and also have smaller amplitudes with respect to $E_{ww}$ for the $Fr=1$ case.

Figure 5.

Energy spectra of (a) lateral $v$ fluctuations and (b) vertical $w$ fluctuations at a downstream point ( $x_{1}=1.6,x_{2}=0.51,x_{3}=0$ ) in the horizontal centre plane at various Froude numbers. Here, $E_{vv}$ , $E_{ww}$ and the Strouhal number, $St$ , are non-dimensional values based on  $U$ and  $D$ .

The mean velocity profiles change significantly with decreasing $Fr$ because of the preferential flow around the sphere rather than over it. Thus, the profile of the mean streamwise velocity (figure 6 a) along the lateral line $(x_{1}=x_{3}=0,x_{2}>0.5)$ shows enhanced horizontal shear in the vicinity of the sphere boundary at $x_{2}=0.5$ for the lower- $Fr$ cases in comparison with $Fr=1$ . At $x_{1}=1$ (figure 6 b), the shear is confined within a narrow band of $0.5<x_{2}<0.8$ for $Fr=1$ , whereas $Fr=0.25$ , 0.125 and 0.025 show progressively broader regions of shear. The lateral horizontal motion of the fluid near the sphere is also enhanced, as shown by the profile of the lateral velocity $U_{2,mean}(x_{2})$ on the line ( $x_{1}=x_{3}=0,x_{2}>0.5$ ) in figure 6(c). At $x_{1}=1$ , the variation of $U_{2,mean}$ as a function of $x_{2}$ (figure 6 d) is substantial for $Fr=0.25$ , 0.125 and 0.025 and has a complex shape because of the three-dimensional mean flow near the body.

Figure 6.

Streamwise ( $U_{1,mean}$ ) and lateral ( $U_{2,mean}$ ) mean velocity profiles are plotted as a function of the lateral coordinate $x_{2}$ at two streamwise locations ( $x_{1}/D=0,1$ ) in the horizontal central plane, $x_{3}=0$ . The plotted velocity has been normalized with $U$ and the $x_{2}$ coordinate with  $D$ .

Figure 7.

Shear production components for different $Fr$ cases, integrated over $x_{2}{-}x_{3}$ planes. The plotted production components have been normalized with $U$ and  $D$ .

The production of TKE is given by $P=-\overline{u_{i}^{\prime }u_{j}^{\prime }}\unicode[STIX]{x2202}\bar{U}_{i}/\unicode[STIX]{x2202}x_{j}$ , with the overbar denoting a mean value. The various components $P_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}}$ that comprise $P$ change in the near wake ( $x/D<5$ ) because of the buoyancy effect. Figure 7 shows the downstream evolution of the components $P_{\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}}$ integrated over the cross-stream $x_{2}{-}x_{3}$ plane. The integrated production for the $Fr=1$ wake is primarily dominated by the components $(P_{13},P_{31})$ involving vertical fluctuations $u_{3}^{\prime }$ , with some contributions from the components $(P_{12},P_{22})$ involving horizontal fluctuations $u_{2}^{\prime }$ , as shown in figure 7(a). This scenario changes when the stratification increases. As illustrated in figure 7(b) for $Fr=0.25$ , the components $P_{13}$ and especially $P_{31}$ are suppressed with respect to $Fr=1$ , and by $Fr=0.025$ (figure 7 d), both become negligible as the buoyancy effect strengthens to make $u_{3}^{\prime }$ negligible. However, $P_{12}$ and $P_{22}$ associated with horizontal fluctuations increase when $Fr$ is reduced to 0.125 from 0.25. The large lateral ( $x_{2}$ ) gradients of mean $U_{1}$ (figure 6 b) and mean $U_{2}$ (figure 6 d) enhance $P_{12}$ and $P_{22}$ respectively, making them the leading production terms for $Fr=0.125$ and 0.025.

The buoyancy Reynolds number, $Re_{b}=\unicode[STIX]{x1D700}/\unicode[STIX]{x1D708}N^{2}$ , where $\unicode[STIX]{x1D700}$ is the turbulent dissipation rate and $N$ is the background buoyancy frequency, is an often-used parameter to distinguish the turbulent nature of fluctuations in stratified flow. A similar parameter that distinguishes turbulence is $\mathscr{R}=ReFr_{h}^{2}$ , where $Fr_{h}=u/l_{h}N$ ( $l_{h}$ is the length scale and $u$ is the velocity scale of horizontal fluctuations) is the horizontal Froude number and $Re=ul_{h}/\unicode[STIX]{x1D708}$ . The choice of $l_{h}=u^{3}/\unicode[STIX]{x1D700}$ makes $\mathscr{R}$ identical to $Re_{b}$ . Riley & deBruynKops (Reference Riley and deBruynKops2003) estimated the Richardson number of layered motions in strongly stratified flow by $Ri\simeq 1/\mathscr{R}$ , and proposed that layer instability was possible if $Ri\lesssim 1$ or, equivalently, $\mathscr{R}\gtrsim 1$ . Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007) concluded that if $\mathscr{R}\gg 1$ , an energy cascade from large to small scales is possible, allowing an inertial range in horizontal energy spectra. In contrast, for $\mathscr{R}\ll 1$ , the dissipation $\unicode[STIX]{x1D700}$ is associated with quasi-two-dimensional scales. Arobone & Sarkar (Reference Arobone and Sarkar2010), in their DNS of a stratified fluid with horizontal shear, found a network of quasi-2D vortices with interspersed dislocations that were laminar for small $Re_{b}$ but exhibited secondary instability for larger  $Re_{b}$ .

We find that the values of $Re_{b}$ (figure 8) provide guidance to the observed differences in the state of fluctuating motion at different $Fr$ . The $Fr=1$ case has $Re_{b}$ values between 10–100 at $0.54<x_{1}/D<5.5$ , signifying broadband turbulence, as observed from the energy content at high frequencies in the horizontal and vertical energy spectra (figure 5 a,b). For the lower $Fr$ of 0.25, the streamwise locations $0.5<x_{1}/D<3$ have $0.1<Re_{b}<1$ . At these streamwise locations, the vortices are still attached, as shown in figure 3(c), and no small-scale features are present. Some of the small scales observed in the $Fr=0.25$ case (figure 3 c) at $x_{1}/D=4{-}5$ are consistent with $Re_{b}\gtrsim 1$ in this region. The small scales observed in figure 3(e) are consistent with the $O(1)$ values of $Re_{b}$ for $Fr=0.125$ at locations $1.14<x_{1}/D<2.75$ , where $Re_{b}<1$ and the flow transitions towards quasi-2D dissipation. For $Fr=0.025$ , $Re_{b}\ll 1$ at all $x_{1}/D$ locations. There is vertical shear between pancake eddies, as shown in figure 3(f,h), which is quasi-laminar for small $Re_{b}$ , consistent with Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007). Nevertheless, the flow is far from laminar. The horizontal motion is unsteady due to vortex shedding, there is broadband turbulence in the near wake, as shown by velocity spectra, and there are small scales, e.g. thin braid vortices between the vortices being shed from the sphere (figure 3 g) in the vorticity field.

Figure 8.

Variation of the buoyancy Reynolds number $Re_{b}=\unicode[STIX]{x1D700}/(\unicode[STIX]{x1D708}N^{2})$ for different $Fr$ at the centre line $x_{2}=0,x_{3}=0$ in the streamwise direction $x_{1}$ .

From figure 8, it can be seen that for $Fr=0.25$ and 0.125, the value of $Ri\approx 1/Re_{b}$ is ${\lesssim}1$ and, therefore, secondary KH instabilities are present in the vertical layers (figure 3 d,f). However, for $Fr=0.125$ at $x_{1}/D>5$ , the value of $Ri>1$ , and for $Fr=0.025$ , the value of $Ri\gg 1$ at all $x_{1}/D$ locations. Hence, secondary instability is absent in the vertical layers at the $x_{1}/D\approx 5$ location in figure 3(f) and at all locations in figure 3(h).