2.1. Flow over a bluff body

We consider the flow over a bluff body, which is depicted schematically in figure 1. A solid isolated body occupying a spatial volume $B$ is held stationary. The space outside $B$ is filled with a viscous incompressible fluid with viscosity $\nu ^{*}$ , where the superscript star designates dimensional quantities. The non-dimensional form of the governing Navier–Stokes equations is

(2.1) \begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} & = 0, \end{align}

(2.2) \begin{align} \frac {\partial \boldsymbol{u} }{\partial t} & = \boldsymbol{u} \times \boldsymbol{\omega } - \boldsymbol{\nabla} h + \frac {1}{Re} \boldsymbol{\nabla} ^2 \boldsymbol{u}, \end{align}

where $\boldsymbol{u}=(u,v,w)$ is the velocity vector, and the generalised enthalpy $h \equiv ({p}/{\rho })+ ({1}/{2})|\boldsymbol{u}|^2$ is the sum of the pressure (divided by density that is normalised to unity) and kinetic energy per unit mass. The bulk Reynolds number is $Re=\left |\boldsymbol{U}^{*}\right |\,L^{*}/\nu ^{*}$ , where $\boldsymbol{U}^{*}$ is the free-stream velocity. The characteristic length $L^{*}$ is defined in terms of the size of the bluff body, using the diameter of the sphere and the length of the minor axis for the spheroid. Since the bluff body is stationary, the velocity field at $\partial B$ satisfies the no-slip boundary conditions. In the far field $|\boldsymbol{x}| \rightarrow \infty$ , the flow approaches the free-stream velocity. These boundary conditions are expressed as

(2.3) \begin{align} \left .\boldsymbol{u}\right |_{\partial B}=0, \quad \boldsymbol{u} \underset {|\boldsymbol{x}| \rightarrow \infty }{\sim } \boldsymbol{U}, \end{align}

where $\boldsymbol{U} = \boldsymbol{U}^{*}/\left |\boldsymbol{U}^{*}\right |$ is the non-dimensional free-stream velocity. The body exerts a drag force on the fluid, which is expressed as the surface integral of pressure and wall shear stress,

(2.4) \begin{align} \boldsymbol{F} = \int _{\partial B}\left (p\hat {\boldsymbol{n}}-2\mu \boldsymbol{S}\hat {\boldsymbol{n}}\right )\,\mathrm{d}S, \quad \boldsymbol{S}=\tfrac {1}{2}\left(\boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{\nabla} \boldsymbol{u}^{\top}\right), \end{align}

where $\hat {\boldsymbol{n}}$ is the unit normal vector on the surface of the body pointing into the fluid, $\mu$ is the dynamic viscosity, and $\boldsymbol{S}$ is the strain rate tensor. The definition of $\boldsymbol{F}$ is the force exerted by the body on the fluid, consistent with Eyink (Reference Eyink2021). We adopt this convention because the momentum and energy balances considered in the derivation of the J–A relation are for the fluid domain. The governing equations (2.1) and initial and boundary conditions provide a complete description of the velocity field and of the drag force.

Figure 1.

Schematic of the flow over a bluff body. A uniform flow with velocity $\boldsymbol{U}$ passes over a solid body $B$ (coloured in blue). The vortical structures are defined by the iso-surface of the Q-criteria $Q=0.5$ , and are coloured by enstrophy.

Taking the curl of the momentum equation (2.2) yields the Helmholtz equation for the vorticity $\boldsymbol{\omega }\equiv \boldsymbol{\nabla} \times \boldsymbol{u}$ :

(2.5) \begin{align} \frac {\partial \boldsymbol{\omega }}{\partial t} + \boldsymbol{\boldsymbol{\nabla} } \boldsymbol{\cdot} \boldsymbol{\Sigma } = 0, \end{align}

where

(2.6) \begin{align} \boldsymbol{\Sigma } = \boldsymbol{u}\boldsymbol{\omega }-\boldsymbol{\omega }\boldsymbol{u}-\nu \left(\boldsymbol{\nabla} \boldsymbol{\omega } - \boldsymbol{\nabla} \boldsymbol{\omega }^{\top}\right) \end{align}

is the Huggins flux tensor. The Helmholtz equation (2.5) is written as a conservation law by aid of the Huggins flux tensor $\boldsymbol{\Sigma }$ , which encodes the transport of vorticity by advection, stretching and tilting, and viscous diffusion. Specifically, the entry $\Sigma _{ij}$ represents the flux of the $j$ th vorticity component in the $i$ th coordinate direction. The Huggins tensor is anti-symmetric, and thus admits a dual representation by an axial vector $\boldsymbol{\eta }$ :

(2.7) \begin{align} \boldsymbol{\eta }=\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\nabla} \times \boldsymbol{\omega }, \quad \eta _i = \tfrac {1}{2}\varepsilon _{ijk}\Sigma _{jk},\quad \Sigma _{ij}=\varepsilon _{ijk}\eta _k. \end{align}

Note that the vorticity-flux tensor in (2.5) is not uniquely defined. Any tensor $\boldsymbol{\Sigma }' = \boldsymbol{\Sigma } + \boldsymbol{\Theta }$ with $\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\Theta }=0$ is a valid option that describes the same vorticity evolution as $\boldsymbol{\Sigma }$ . However, the Huggins flux tensor $\boldsymbol{\Sigma }$ is the unique choice that can be connected with momentum transport and pressure gradient. This connection is given by the rotational form of the momentum equation,

(2.8) \begin{align} \frac {\partial u_i}{\partial t}=\frac {1}{2}\epsilon _{ijk}\Sigma _{jk}-\frac {\partial h}{\partial x_i}. \end{align}

Equation (2.8) provides a direct relation between the vorticity flux and the total pressure gradient, which at stationary walls reduces to

(2.9) \begin{align} -\nu\, \boldsymbol{\nabla} \times \boldsymbol{\omega }=\boldsymbol{\nabla} p. \end{align}

An expression for the wall vorticity flux $\boldsymbol{\sigma }=\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\Sigma }$ follows by taking the cross-product between the wall-normal unit vector and (2.9):

(2.10) \begin{align} \boldsymbol{\sigma }=\boldsymbol{n}\boldsymbol{\cdot} \boldsymbol{\Sigma }=\boldsymbol{n}\times \left (\nu\, \boldsymbol{\nabla} \times \boldsymbol{\omega }\right )=-\boldsymbol{n}\times \boldsymbol{\nabla} p, \end{align}

where $\boldsymbol{n}= -\hat {\boldsymbol{n}}$ is the unit normal vector on the wall pointing towards the body. The physical meaning of $\boldsymbol{\sigma }$ is the rate of vorticity transport from the fluid out through the wall per unit surface area. Integrating the vorticity equation (2.5) in a finite domain enclosed by a surface, the time rate of change of vorticity in this domain is balanced by the surface integral of $\boldsymbol{\sigma }$ . According to (2.10), the viscous vorticity flux at the wall is closely related to the tangential pressure gradient, as discussed by Lighthill (Reference Lighthill1963) and Morton (Reference Morton1984). In this study, we numerically evaluate the Huggins vorticity flux (2.6), and use it to examine the transport of vorticity. We also explore the role of vorticity flux in the momentum balance using (2.8).

We end this subsection with a brief discussion of empirical observations related to boundary-layer separation and vorticity. The boundary layer develops along the surface of the body from the attachment line and evolves downstream, and the curvature of the surface induces favourable and adverse pressure gradients. Since the Reynolds numbers in the numerical studies considered in this work are subcritical – $Re_c=3\times 10^{5}$ for flow over a sphere (Achenbach Reference Achenbach1974), and $4.2\times 10^{5}$ for flow over a spheroid (Ahn Reference Ahn1992) – natural transition to turbulence does not take place prior to separation. Furthermore, depending on the geometry and inflow conditions, the boundary layer undergoes either two- or three-dimensional separation, as can be gleaned from the wall shear stress $\boldsymbol{\tau }_w=2\mu \boldsymbol{S}\hat {\boldsymbol{n}}$ (Tobak & Peake Reference Tobak and Peake1982). For instance, the steady axisymmetric flow over a sphere at $24\leqslant Re (:= ({U^{*}D^{*}}/{\nu ^{*}}) )\leqslant 200$ undergoes two-dimensional separation at a polar angle $0^{\circ }\lt \theta \lt 62^{\circ }$ where $\tau _w=0$ (Johnson & Patel Reference Johnson and Patel1999). For an example of three-dimensional separation, consider flow over a prolate spheroid with a non-zero angle of attack (Wang et al. Reference Wang, Zhou, Hu and Harrington1990). Depending on the angle of incidence, multiple separation and reattachment patterns can develop and result in a complex topology of wall-stress lines on the spheroid surface. Separation lines are identified as limiting friction lines connecting two singular points of the wall shear stress (where $\tau _w=0$ ), and have neighbouring wall shear stress lines converging towards them (Chapman & Yates Reference Chapman and Yates1991; Surana, Grunberg & Haller Reference Surana, Grunberg and Haller2006).

In both two- and three-dimensional separation, the motion of vorticity is key to understanding the boundary-layer dynamics. First, the vorticity is directly related to the wall shear stress by $\boldsymbol{\tau }_w=\mu \boldsymbol{\omega }\times \hat {\boldsymbol{n}}$ , thus the characterisation of separation using the wall shear stress can be equivalently established in terms of the vorticity (Wu et al. Reference Wu, Tramel, Zhu and Yin2000). Second, the wall pressure gradient is related to the vorticity flux through (2.10) (Wu, Ma & Zhou Reference Wu, Ma and Zhou2007). Favourable pressure gradient creates vorticity with the same orientation as that within the boundary layer, which keeps the boundary layer attached to the wall; adverse pressure gradient creates vorticity with opposite orientation to the existing vorticity within the boundary layer, and thus promotes flow separation from the wall.

2.2. The detailed J–A relation

We start by stating the J–A relation for flow over a stationary isolated body, then proceed to summarise its derivation. For the interested reader, a detailed derivation is provided by Eyink (Reference Eyink2021). In the lab frame, the power injected into the fluids by the drag force is given by $-\boldsymbol{F}\boldsymbol{\cdot} \boldsymbol{U}$ . This quantity, in the body frame, can be expressed as

(2.11) \begin{eqnarray} -\boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{U} &=& -\rho \int _{\Omega }\boldsymbol{u}_{\phi }\boldsymbol{\cdot} \left (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )\,\mathrm{d}V \end{eqnarray}

(2.12) \begin{align} &= -\int \mathrm{d}J\int \left (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )\boldsymbol{\cdot} \mathrm{d}\boldsymbol{l} \end{align}

(2.13) \begin{align} & = -\tfrac {1}{2}\int \mathrm{d}J\int \varepsilon _{ijk}\Sigma _{ij}\,\mathrm{d}l_k. \end{align}

This power injection is equal to the integral over the fluid volume $\Omega$ of the vorticity flux against a background potential flow $\boldsymbol{u}_{\phi }$ , as shown by (2.11). The volume element is decomposed into $\mathrm{d}V = \mathrm{d}A\,|\mathrm{d}\boldsymbol{l}|$ , where $\mathrm{d}\boldsymbol{l}$ is a vector line element along $\boldsymbol{u}_{\phi }$ , and $\mathrm{d}A$ is an area element normal to $\boldsymbol{u}_{\phi }$ . By defining $\mathrm{d}J = \rho\, |\boldsymbol{u}_{\phi }|\,\mathrm{d}A$ as the potential mass current within a streamtube, the right-hand-side of the J–A relation is cast into (2.12). Using the identities in (2.7), the J–A relation is then written in terms of the Huggins flux tensor $\boldsymbol{\Sigma }$ in (2.13). We consider this right-hand side starting with $\mathrm{d}l_k$ , which is the line element aligned with the potential flow, and notice that $\varepsilon _{ijk}\Sigma _{ij}\,\mathrm{d}l_k$ vanishes if either $i$ or $j$ is equal to $k$ . The implication is that contributions to the integral arise only due to components in the Huggins tensor that correspond to the flux ( $i$ index) and the vorticity ( $j$ index) both being orthogonal to the potential flow, i.e. when $i \ne k$ and $j \ne k$ . Hence (2.11)–(2.13) represent the the rate of energy injection associated with the vorticity flux crossing the potential-flow streamlines, weighted by the potential flow speed. In other words, the rate of drag work performed by the immersed body on the fluid is equal to the amount of vorticity normal to the potential flow that crosses the potential mass current outwards. Conversely, inward vorticity flux across the potential-flow streamlines reduces drag power. Finally, vorticity flux along the potential-flow streamlines does not contribute to rate of work by drag.

The derivation of the J–A relation starts by defining $(\boldsymbol{u}_{\phi }, p_{\phi })$ as the potential-flow solution over the solid body with the same free-stream configuration, which can be obtained either analytically or numerically. The true velocity and pressure fields from the Navier–Stokes solution are then split into potential and vortical parts, with the vortical solution defined as

(2.14) \begin{align} \boldsymbol{u}_{\omega } := \boldsymbol{u}-\boldsymbol{u}_{\phi }, \quad p_{\omega } := p-p_{\phi }. \end{align}

The governing equations for the vortical flow $(\boldsymbol{u}_{\omega }, p_{\omega })$ are derived by subtracting the Euler and Navier–Stokes equations, which yields

(2.15) \begin{align} \partial _t \boldsymbol{u}_\omega =\boldsymbol{u} \times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} } \times \boldsymbol{\omega }-\boldsymbol{\boldsymbol{\nabla} } h_\omega , \end{align}

where $h_\omega =p_\omega +({1}/{2})\left |\boldsymbol{u}_\omega \right |^2+\boldsymbol{u}_\omega \boldsymbol{\cdot} \boldsymbol{u}_\phi$ is the total pressure for the vortical flow. A relation between the total vortical momentum $\boldsymbol{P}_{\omega }$ , the total vortical drag $\boldsymbol{F}_{\omega }$ , and the far-field vortical pressure is obtained by integrating (2.15) in space:

(2.16) \begin{align} \frac {\mathrm{d} \boldsymbol{P}_\omega }{\mathrm{d} t}=\boldsymbol{F}_\omega -\rho \lim _{R \rightarrow \infty } \int _{S_R} \hat {\boldsymbol{x}} h_\omega\, \mathrm{d} S, \quad \boldsymbol{P}_{\omega } =\int _{\Omega }\rho \boldsymbol{u}_{\omega }\,\mathrm{d}V ,\quad \boldsymbol{F}_{\omega } = \int _{\partial B}(p_{\omega }\hat {\boldsymbol{n}}-2\mu \boldsymbol{S}\hat {\boldsymbol{n}})\,\mathrm{d}S, \\[-16pt] \nonumber \end{align}

where $S_R$ is a sphere with radius $R$ , and $\hat {\boldsymbol{x}} = \boldsymbol{x}/|\boldsymbol{x}|$ is the radially outward unit vector. The total vortical momentum $\boldsymbol{P}_{\omega }$ represents the amount of momentum contained in the vortical flow. The total vortical force $\boldsymbol{F}_{\omega }$ is the same as the drag in (2.4). The first expression of (2.16) shows that the rate of change of the total vortical momentum is driven by the imbalance between the drag force $\boldsymbol{F}_{\omega }$ and the far-field vortical pressure. This time derivative ${{\rm d}\boldsymbol{P}_\omega }/{{\rm d}t}$ is generally not zero. In addition to the momentum equation (2.16), we also consider the energy equation. The following expression for the local kinetic energy is obtained from the potential–vortical decomposition of velocity:

(2.17) \begin{align} e = \frac {1}{2}\rho \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{u} = \frac {1}{2}\rho \boldsymbol{u}_{\phi }\boldsymbol{\cdot} \boldsymbol{u}_{\phi }+\rho \boldsymbol{u}_{\phi }\boldsymbol{\cdot} \boldsymbol{u}_{\omega } + \frac {1}{2}\rho \boldsymbol{u}_{\omega }\boldsymbol{\cdot} \boldsymbol{u}_{\omega }. \end{align}

An energy evolution equation can then be written for each of the three ingredients, namely the potential, interaction and vortical kinetic energies. The volume integral of the potential kinetic energy is infinite and conserved. We are specifically interested in the evolution equation of the interaction energy $\boldsymbol{u}_{\phi } \boldsymbol{\cdot} \boldsymbol{u}_{\omega }$ . Performing a dot product of the governing equation for $\boldsymbol{u}_{\omega }$ with $\boldsymbol{u}_{\phi }$ , we obtain

(2.18) \begin{align} \begin{aligned} \partial _{t}\left (\boldsymbol{u}_{\phi } \boldsymbol{\cdot} \boldsymbol{u}_{\omega }\right )+\boldsymbol{\boldsymbol{\nabla} } \boldsymbol{\cdot} \left [h_{\omega } \boldsymbol{u}_{\phi }+h_\phi \boldsymbol{u}_{\omega }\right ]=\boldsymbol{u}_{\phi } \boldsymbol{\cdot} (\boldsymbol{u} \times \omega -\nu\, \boldsymbol{\boldsymbol{\nabla} } \times \boldsymbol{\omega }), \end{aligned} \end{align}

where $h_\phi =p_\phi +({1}/{2})\left |\boldsymbol{u}_\phi \right |^2$ is the potential total pressure. The divergence term spatially transports the local interaction energy, and does not cause any net change in the total interaction energy. The right-hand side is the work done by the vortex force $(\boldsymbol{u} \times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} } \times \boldsymbol{\omega })$ along the potential flow $\boldsymbol{u}_{\phi }$ . Similar to momentum, (2.18) is integrated over the fluid domain:

(2.19) \begin{align} \frac {\mathrm{d}}{\mathrm{d} t}\int _{\Omega }\rho \boldsymbol{u}_{\phi } \boldsymbol{\cdot} \boldsymbol{u}_{\omega }\,\mathrm{d}V=\rho \int _{\Omega } \boldsymbol{u}_\phi \boldsymbol{\cdot} (\boldsymbol{u} \times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} } \times \boldsymbol{\omega })\, {\rm d} V-\rho \boldsymbol{U} \boldsymbol{\cdot} \lim _{R \rightarrow \infty } \int _{S_R} \hat {\boldsymbol{x}} h_\omega\, {\rm d}A. \\[-16pt] \nonumber \end{align}

Additionally, a multi-pole expansion of the vortical velocity field $\boldsymbol{u}_{\omega }$ provides the following expression for the total interaction energy:

(2.20) \begin{align} \int _{\Omega }\boldsymbol{u}_{\phi } \boldsymbol{\cdot} \boldsymbol{u}_{\omega }\,\mathrm{d}V = \int _{\Omega }\boldsymbol{u}_{\omega }\,\mathrm{d}V\boldsymbol{\cdot} \boldsymbol{U}. \end{align}

Combining (2.16), (2.19) and (2.20) yields the J–A relation (2.11), which is repeated below:

(2.21) \begin{align} -\boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{U} = \rho \int _{\Omega }\underbrace {\boldsymbol{u}_{\phi }\boldsymbol{\cdot} \left (-\boldsymbol{u}\times \boldsymbol{\omega }\right )}_{\Pi _a}+\underbrace {\boldsymbol{u}_{\phi }\boldsymbol{\cdot} \left (\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )}_{\Pi _\nu }\,\mathrm{d}V. \end{align}

This relation introduces a new perspective on the transport of vorticity, the transfer of energy, and the drag power. First, the fields $\Pi _{a}=-\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\boldsymbol{u}\times \boldsymbol{\omega } )$ and $\Pi _{\nu }=\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega } )$ are the energy fluxes that correspond to advective and viscous vorticity transport crossing the potential streamlines of $\boldsymbol{u}_{\phi }$ . Positive values represent advection and diffusion effects transporting negative vorticity outwards across potential streamlines, and the reverse for negative values. Second, $\Pi _{a}$ and $\Pi _{\nu }$ can be interpreted as spatial contributions to the power injection into the fluid by the drag force, with positive and negative values being the drag and anti-drag contributions. Finally, comparing the right-hand side of the J–A relation (2.11) to the interaction-energy equation (2.18) shows that drag is associated with a loss in the interaction energy. This loss appears as a source in the vortical energy equation:

(2.22) \begin{align} \begin{aligned} &\partial _{t}\left (\tfrac {1}{2}\left |\boldsymbol{u}_{\omega }\right |^{2}\right )+\boldsymbol{\boldsymbol{\nabla} } \boldsymbol{\cdot} \left [\left (p_{\omega }+\tfrac {1}{2}\left |\boldsymbol{u}_{\omega }\right |^{2}+\boldsymbol{u}_{\omega } \boldsymbol{\cdot} \boldsymbol{u}_{\phi }\right ) \boldsymbol{u}_{\omega }-\nu \boldsymbol{u} \times \boldsymbol{\omega} \right ] \\ &=-\boldsymbol{u}_{\phi } \boldsymbol{\cdot} (\boldsymbol{u} \times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} } \times \boldsymbol{\omega })-\nu\, |\boldsymbol{\omega }|^{2}. \end{aligned} \end{align}

In other words, $\Pi _a + \Pi _{\nu } = -\boldsymbol{u}_{\phi } \boldsymbol{\cdot} (\boldsymbol{u} \times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} } \times \boldsymbol{\omega })$ is exactly equal to the energy transfer from the interaction energy to the vortical energy. Positive values are transfers from interaction to vortical energy by advection and viscous effects, and the reverse for negative values. Ultimately, the vortical energy is dissipated into heat by the last term, $-\nu\, |\boldsymbol{\omega }|^2$ , in (2.21).

The physical interpretation of the J–A relation in a finite integration domain demands more careful discussion. A steady or statistically stationary state of the flow near the body can be reached at a sufficiently long time after an initial impulsive start. However, the flow far downstream remains non-stationary, for example in the region where the starting vortex is first observed. As such, the integral of the unsteady term in (2.18) does not have a contribution from the stationary flow near the body, and is solely due to the far wake. In other words, the rate of change of total interaction energy is primarily due to the far wake where the flow has not reached stationarity. The streamwise extent of this region of the wake is at least the product of the free-stream velocity and the transient time (from the initial condition to the stationary state), which is too large to include in DNS. Nevertheless, the J–A relation can still be approximately satisfied within the domains customarily adopted for DNS, because the right-hand side of (2.18) decays fast downstream – an argument that is numerically verified in this study.

Figure 2.

Computational domains and meshes. $(a)$ The flow set-up, the multi-block grid system, and a rotated view of the front block for the flow over the sphere. (b) The same visualisations for the flow over the prolate spheroid.

2.3. Numerical simulation of flow over a bluff body

In this subsection, we describe the numerical approach for simulating the flows over a sphere and a spheroid. The governing equations (2.1) are discretised and solved using a fractional-step approach with a local volume-flux formulation on a staggered curvilinear grid (Wang, Wang & Zaki Reference Wang, Wang and Zaki2019; You & Zaki Reference You and Zaki2019). The advection terms are discretised using the Adams–Bashforth scheme, and the Crank–Nicolson scheme is adopted for the diffusion terms. The pressure Poisson equation is solved using a bi-conjugate gradient stabilised method (BICGSTAB) with an algebraic multi-grid preconditioner provided by Hypre (Falgout & Yang Reference Falgout and Yang2002). These algorithms are first implemented and validated on a single-block domain. For the simulation of flow over a bluff body, the fluid domain is decomposed into six blocks, as shown in figure 2. Within each block, the advection and diffusion terms are discretised using the single-block description noted above, and on the block boundaries the diffusion terms are discretised using the Adams–Bashforth scheme. The pressure fields in all blocks are solved globally due to the ellipticity of the pressure Poisson equation.

The computational domains and grids for the flows over the sphere and spheroid are shown in figure 2. The fluid domain in the former case is formed by two concentric sphere surfaces with $R_1 = 0.5$ and $R_2 = 15$ , and divided into six blocks for the structured multi-block flow solver. No-slip and free-stream boundary conditions are imposed for the velocity field at the inner and outer sphere surfaces, respectively. Each block is discretised into a structured curvilinear mesh. The number of grid points on the interface between two blocks is the same on the two sides, which constrains the number of grid points within different blocks. On account of these constraints, there are four independent numbers of grid points $N_x, N_y, N_z, N_w$ that can be specified on this multi-block grid, which are shown in figure 1. The flow domain, mesh and boundary conditions for the flow over the prolate spheroid are configured similarly. The fluid domain is formed between two concentric spheroids with their axes aligned, decomposed into six blocks, and discretised on a structured curvilinear mesh similar to the sphere cases. The aspect ratio of the inner spheroid is $a/b=6:1$ , as shown in figure 1(b). The outer spheroid has aspect ratio close to unity, and radii equal to $16.44$ and $16.17$ , with its major axes aligned with the inner spheroid. The geometry, mesh parameters and Reynolds numbers are reported in table 1. The case designations ‘SL’, ‘ST’, ‘PT’ refer to laminar flow over the sphere, turbulent flow over the sphere, and turbulent flow over the prolate spheroid.

Table 1.

Geometries, Reynolds numbers, number of grid points, and resolutions for DNS. The resolution $\Delta y_{b}$ represents the wall-normal grid spacing at the solid wall, while $\unicode{x1D6E5} _w = (\Delta x_w\, \Delta y_w\, \Delta z_w)^{1/3}$ denotes the grid size at a point in the wake three units of length downstream of the trailing edge of the sphere and spheroid.

Figure 3.

Vortical structures visualised using iso-surface of the Q-criteria and coloured by the streamwise velocity, from simulations of the flows over the spheres and the spheroid. $(a)$ Flow over a sphere at $Re=200$ , $Q=0.5$ . $(b)$ Flow over a sphere at $Re=3700$ , $Q=0.1$ . $(c)$ Flow over a spheroid at $Re=3000$ , $Q=0.5$ .

Visualisations of the vortical structures in the three flows are shown in figure 3. The distinct characteristics of the vortical structures emphasise the different purposes of these three examples. The simplicity of case SL, reflected by figure 3(a), enables a concise demonstration of the theoretical elements discussed in §§ 2.1 and 2.2. Vortex shedding and a turbulent wake are introduced by considering case ST (figure 3 b), where a statistical description of the vorticity dynamics is required. The pair of vortices above the spheroid in case PT (figure 3 c) is closely related to the three-dimensional separation on the wall. The implication for vorticity transport by these vortices and the underlying separation patterns are the focus of case PT.

The J–A relation (2.11) and the Huggins flux tensor (2.6) are evaluated numerically for the three simulated flows described earlier. The vorticity is computed using a finite-volume scheme on a generalised curvilinear coordinate system (Rosenfeld, Kwak & Vinokur Reference Rosenfeld, Kwak and Vinokur1991). This vorticity is then used to evaluate the advection term in the J–A relation and the advective vorticity flux. The vorticity diffusion vector is computed based on the identity $\boldsymbol{\nabla} ^2\boldsymbol{u}=-\boldsymbol{\nabla} \times \boldsymbol{\omega }$ , where the Laplacian of the velocity field is obtained from the right-hand side of the momentum equation in the Navier–Stokes solver. This procedure ensures that the numerical evaluation of vorticity diffusion is consistent with the momentum balance enforced in the solver.

Both cylindrical and spherical coordinates will be utilised to present the numerical results, and they are shown schematically in figure 4. The origins of both coordinates are placed at the centre of the sphere for cases SL and ST, and the free-stream velocity lies along the $x$ -axis. For the flow over the spheroid, the centre of the spheroid is located at $(x,y,z)=({a}/{2},0,0)$ , and the major axis lies with the $x$ -axis. The free-stream velocity is $\boldsymbol{U}=(\cos (\alpha ),\sin (\alpha ),0)$ , where $\alpha$ is the incidence angle.

Figure 4.

Schematics of the $(a)$ cylindrical and $(b)$ spherical coordinate systems that are adopted in the analysis of vorticity fluxes. $(a)$ The $x$ -axis is aligned with the polar coordinate, the azimuthal angle is denoted by $\varphi$ , and the radial coordinate is $r$ . $(b)$ The polar angle $\theta$ is formed by the polar axis ( $x$ -axis) and the radial vector. The length of the radial vector is denoted by $\eta$ . The azimuthal angle $\varphi$ is formed with respect to the $y$ -direction.

3.1. Laminar flow over a sphere

The flow over a sphere at Reynolds number $Re=200$ is steady and axisymmetric. The forebody of the sphere is wrapped with an attached laminar boundary layer that separates at approximately $\theta = 62^{\circ }$ due to the adverse pressure gradient and curvature. A steady cylindrical shear layer that wraps the recirculation region forms behind the sphere. Farther downstream, the flow recovers from the defect profile towards the uniform free-stream velocity.

Figure 5.

Time history of drag coefficient for case SL from integration of the wall pressure and shear stress J–A relation. : total drag work evaluated by surface integration of pressure and friction. : pressure work. : friction work. : total drag work evaluated from the J–A relation. : total advective contribution $\int _{\Omega }\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (-\boldsymbol{u}\times \boldsymbol{\omega })\,\mathrm{d}V$ . : total viscous contribution $\int _{\Omega }\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega })\,\mathrm{d}V$ .

We compare two approaches to evaluating the power injection into the fluid $-\boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{U}$ . First, the drag force is evaluated by integrating the pressure and shear stress over the sphere surface using (2.4). In the second approach, we evaluate the J–A relation (2.11). In the latter case, the potential flow is given by the analytical expression

(3.1) \begin{align} \phi =U\left (\eta +\frac {R_1^{3}}{2 \eta ^{2}}\right ) \cos \theta , \qquad \psi =\frac {1}{2} U \eta ^{2} \sin ^{2} \theta \left (1-\frac {R_1^{3}}{\eta ^{3}}\right ), \end{align}

where $\eta$ is the radial distance in the spherical coordinate system. The comparison in figure 5 demonstrates that the two approaches to evaluating $-\boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{U}$ agree, both yielding equal time histories. The ingredients in terms of the wall pressure and wall shear stress, as well as from the vorticity fluxes by advection and diffusion, are also visualised in the same figure. The magnitudes of pressure and friction drag are similar and remain constant through time in this case. Instead of dividing the drag force into form and friction drag, the J–A relation expresses the rate of drag work in advective $\int _{\Omega }\Pi _{a}\,\mathrm{d}V$ and diffusive $\int _{\Omega }\Pi _{\nu }\,\mathrm{d}V$ terms. The advective part is most closely related to form drag since the advection of vorticity is accompanied with a pressure gradient in the transverse direction. The diffusive part, on the other hand, is most closely related to the friction drag. These four components of drag have similar magnitudes in figure 5 due to the low Reynolds number considered here. Therefore, it can be anticipated that the comparison between the two approaches to compute drag, and the differences between two contributions within each approach, will be more evident in the turbulent case ST.

Figure 6.

Schematics of potential flow and vorticity transport direction. (a) The translucent surfaces represent the iso-surfaces of $\psi =-0.1,-0.2,-0.4$ . The vectors $\boldsymbol{e}_{\phi }$ , $\boldsymbol{e}_{s}$ , $\boldsymbol{e}_{n}$ form a set of local orthogonal coordinates. (b) The solid black lines with arrows represent the azimuthal vortex rings. The red and blue arrows represent the outward and inward transport of vorticity crossing the iso-surface of $\psi$ .

For axisymmetric flows, the J–A relation (2.12) can be further rewritten using the flux of the azimuthal vorticity. We denote the unit vector in the azimuthal direction as $\boldsymbol{e}_{\varphi }$ , and that along the potential-flow streamlines as $\boldsymbol{e}_{s}$ . Therefore, the unit vector $\boldsymbol{e}_{n}$ normal to an iso-surface of the potential-flow streamfunction is given by $\boldsymbol{e}_{n} = \boldsymbol{e}_{\varphi } \times \boldsymbol{e}_{s}$ , as visualised in figure 6 $(a)$ . The J–A relation is then rewritten as

(3.2) \begin{align} -\boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{U} &= - \int \mathrm{d}J\int \left (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )\mathrm{d}\boldsymbol{l}\nonumber \\&= -\int \mathrm{d}J\int \left (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )\boldsymbol{\cdot} \boldsymbol{e}_{s}\, \mathrm{d}l\nonumber \\&= - \int \mathrm{d}J\int \left (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )\boldsymbol{\cdot} \left (\boldsymbol{e}_{n}\times \boldsymbol{e}_{\varphi }\right ) \mathrm{d}l\nonumber \\&= - \int \mathrm{d}J\int \boldsymbol{e}_{\varphi }\boldsymbol{\cdot} \left (\left (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega }\right )\times \boldsymbol{e}_{n}\right ) \mathrm{d}l \nonumber \\&= - \int \mathrm{d}J\int \boldsymbol{e}_{\varphi }\boldsymbol{\cdot} \boldsymbol{\Sigma }\boldsymbol{\cdot} \boldsymbol{e}_{n}\, \mathrm{d}l \nonumber \\&= - \int \mathrm{d}J\int \boldsymbol{e}_{n}\boldsymbol{\cdot} \boldsymbol{\Sigma }\boldsymbol{\cdot} \left (-\boldsymbol{e}_{\varphi }\right ) \mathrm{d}l, \end{align}

where $\mathrm{d}\boldsymbol{l}=\boldsymbol{e}_{s}\, \mathrm{d}l$ is the vector line element pointing along potential velocity, and $\mathrm{d}l$ is the length of this line element. Figure 6 $(b)$ is a schematic of vortex loops being generated from the sphere surface and transported into the wake. The red and blue arrows represent vortex loops crossing the potential streamline outwards and inwards, and correspond to the positive and negative rates of work by drag, respectively. The expression $\boldsymbol{e}_{n}\boldsymbol{\cdot} \boldsymbol{\Sigma }\boldsymbol{\cdot} (-\boldsymbol{e}_{\varphi } )$ represents the flux of negative azimuthal vorticity in the $\boldsymbol{e}_{n}$ -direction. In this form of the J–A relation, the rate of work done by the drag force is represented as a weighted integral of azimuthal vorticity flux crossing the iso-surface of the streamfunction. The vorticity-flux tensor $\boldsymbol{\Sigma }=\boldsymbol{\Sigma }_a+\boldsymbol{\Sigma }_{\nu }$ can be simplified into the following forms:

(3.3) \begin{align} \boldsymbol{\Sigma }_a = \begin{pmatrix} 0 & 0 & u_x\omega _{\varphi }\\ 0 & 0 & u_r\omega _{\varphi }\\ -u_x\omega _{\varphi } & -u_r\omega _{\varphi } & 0\\ \end{pmatrix}, \quad \boldsymbol{\Sigma }_{\nu } = \begin{pmatrix} 0 & 0 & -\nu \dfrac {\partial \omega _{\varphi }}{\partial x}\\[8pt] 0 & 0 & -\nu \dfrac {\partial \omega _{\varphi }}{\partial r}\\ \nu \dfrac {\partial \omega _{\varphi }}{\partial x} &\nu \dfrac {\partial \omega _{\varphi }}{\partial r} & 0 \end{pmatrix}. \end{align}

The flux of azimuthal vorticity can be obtained by dotting the above expressions with the azimuthal unit vector $\boldsymbol{e}_{\varphi }$ :

(3.4) \begin{align} \boldsymbol{\Sigma }_{a}\boldsymbol{\cdot} \boldsymbol{e}_{\varphi } = \begin{pmatrix} u_x \omega _{\varphi }\\ u_r \omega _{\varphi } \\ 0\\ \end{pmatrix},\quad \boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{\varphi } = \begin{pmatrix} -\nu \dfrac {\partial \omega _{\varphi }}{\partial x} \\[8pt] -\nu \dfrac {\partial \omega _{\varphi }}{\partial r} \\ 0\\ \end{pmatrix}. \end{align}

The explicit forms of $\boldsymbol{\Sigma }_{a}\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ and $\boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ obtained above can be clearly interpreted as the advection of $\omega _{\varphi }$ by the velocity field and the down-gradient diffusion of $\omega _{\varphi }$ by viscosity, respectively. The portion of these fluxes that contributes to the action of drag is the projection onto $\boldsymbol{e}_{n}$ , i.e. the flux normal to the potential flow, as shown in the last expression of (3.2).

Figure 7.

Two-dimensional contour of instantaneous value of vorticity fluxes for case SL overlapped with streamlines. Coloured contours show $(a)$ $\Pi _{a}$ and $(b)$ $\Pi _{\nu }$ . Panels $(a.ii)$ and $(b.ii)$ are zoomed-in views of boxed region $A$ , and panels $(a.iii)$ and $(b.iii)$ are zoomed-in views of boxed region $B$ .

The advective contribution $\Pi _a(\boldsymbol{x})= -\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\boldsymbol{u}\times \boldsymbol{\omega } ) =|\boldsymbol{u}_{\phi }|\,\boldsymbol{e}_{n} \boldsymbol{\cdot} \boldsymbol{\Sigma }_a\boldsymbol{\cdot} (-\boldsymbol{e}_{\varphi } )$ and viscous contribution $ \Pi _{\nu }(\boldsymbol{x}) =|\boldsymbol{u}_{\phi }|\,\boldsymbol{e}_{n} \boldsymbol{\cdot} \boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} (-\boldsymbol{e}_{\varphi } ) = \boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega } )$ are visualised in figure 7, overlapped with flow streamlines. In the near-wall region $A$ , the viscous contribution dominates inside the boundary layer, while the advective contribution is close to zero near the wall due to the no-slip condition. The positive viscous flux contributes to drag power injection in a region along the wall in the attached boundary layer, which is associated with the outward transport of the negative azimuthal vorticity that is introduced into the fluid domain by the favourable pressure gradient. The vortex force in the near-wall region $A$ drives the flow from potential to vortical, and results in a transfer from the interaction energy to vortical energy (interpretation from (2.18)). A weak negative viscous contribution appears near $\theta =90^{\circ }$ at the wall prior to separation, and extends downstream as shown in region $B$ in figure 7. This negative flux, which reduces the rate of drag work, is due to the change in pressure gradient from favourable to adverse and the associated wall vorticity flux $\sigma$ . After the negative azimuthal vorticity diffuses into the domain, it is advected along the detached boundary layer crossing the potential streamlines outwards. This flux drives the flow away from the ideal potential flow, and thus produces a positive contribution to drag power. In the laminar wake, the streamwise velocity recovers from a deficit profile to the free-stream potential flow. The negative contribution to drag power near the wake centreline indicates that the vortex force drives the flow from vortical to potential, and that the kinetic energy is transferred from vortical to interaction energy.

Figure 8.

Visualisation of the vorticity-flux vector field. The lines are tangent to the vector fields (a) $\boldsymbol{\Sigma }_a\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ and (b) $\boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ . The colour contours represent the streamfunction $\psi$ of the background potential flow.

The link between the vorticity flux and the rate of drag work underscores the importance of a quantitative description of vorticity transport, which is given by the Huggins flux tensor. The evolution equation for a general vorticity component $\omega _{\xi } = \boldsymbol{\omega }\boldsymbol{\cdot} \boldsymbol{e}_{\xi }$ can be derived by dotting a constant unit vector $\boldsymbol{e}_{\xi }$ with the vorticity equation (2.5):

(3.5) \begin{align} \frac {\partial \omega _{\xi }}{\partial t} + \boldsymbol{\boldsymbol{\nabla} } \boldsymbol{\cdot} \left (\boldsymbol{\Sigma }\boldsymbol{\cdot} \boldsymbol{e}_{\xi }\right ) = 0. \end{align}

The above relation is a planar conservation law for the transport of $\omega _{\xi }$ within a plane $S_{\boldsymbol{e}_{\xi }}$ that is normal to the unit vector $\boldsymbol{e}_{\xi }$ . Due to the anti-symmetry of the flux tensor $\boldsymbol{\Sigma }$ , the vorticity flux of $\omega _{\xi }$ in the direction $\boldsymbol{e}_{\xi }$ vanishes, since $\boldsymbol{e}_{\xi } \boldsymbol{\cdot} \boldsymbol{\Sigma }\boldsymbol{\cdot} \boldsymbol{e}_{\xi } = 0$ . In other words, in the absence of boundaries, the vorticity component $\omega _{\xi }$ is conserved within the plane $S_{\boldsymbol{e}_{\xi }}$ , as can be seen by integrating (3.5) on the plane $S_{\boldsymbol{e}_{\xi }}$ . The vector field $\boldsymbol{\Sigma }\boldsymbol{\cdot} \boldsymbol{e}_{\xi }$ is the flux of $\omega _{\xi }$ , and serves as an informative visualisation of its transport.

In the present case of laminar separation (SL), the vorticity field is comprised of an azimuthal component only. We therefore consider the azimuthal unit vector $\boldsymbol{e}_{\varphi }$ , and visualise the vector field $\boldsymbol{\Sigma }\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ . The streamtraces parallel to the vector fields $\boldsymbol{\Sigma }_a\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ and $\boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ are visualised in figure 8. The advective flux $\boldsymbol{\Sigma }_a\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ does not start or end at the wall, which implies that this flux only redistributes vorticity within the volume, and thus changes the vorticity magnitude locally without altering the volume integral of vorticity. In contrast, the viscous flux is mainly responsible for the outward diffusion of vorticity from the wall into the fluid and its annihilation at the wake centreline. This interpretation is consistent with the classical Eulerian picture of vorticity transport:

Negative azimuthal vorticity is generated at the wall by pressure gradient, and vorticity diffuses into the boundary layer. Within the cylindrical shear layer, vorticity is mostly advected. The azimuthal vorticities from different azimuthal angles cancel each other at the wake centreline. The inward transport of vorticity observed in figure 8 corresponds to the anti-drag region appearing in figure 7 near the wake centreline. The rate of vorticity annihilation, defined by the vorticity flux $ (\boldsymbol{\Sigma }_\nu )_{r\varphi }$ , can be related to streamwise total pressure gradient by simplifying (2.8) at the wake centreline (Terrington et al. Reference Terrington, Hourigan and Thompson2021):

(3.6) \begin{align} \left (\boldsymbol{\Sigma }_\nu \right )_{r\varphi }=-\nu \frac {\partial \omega _{\varphi }}{\partial r}=\frac {\partial h}{\partial x}. \end{align}

The recovery of total pressure is accompanied by inward transport of azimuthal vorticity towards the wake centreline. This inward flux only contains a viscous contribution at the centreline, meaning that viscous effect is solely responsible for the annihilation of azimuthal vorticity.

3.2. Impulsively started flow over a sphere: multiple two-dimensional separations

The flow over a sphere at $Re=3700$ is simulated and is designated ST (see table 1). The simulation starts from an initial condition with uniform free-stream velocity in the fluid domain everywhere. The time history of the terms in the J–A relation, as well as the rates of form and friction drag work, are visualised in figure 9. The total drag force drops from an unbounded value at $t=0$ to a finite value at later times due to the $t^{-1/2}$ singularity of the drag force at the initial time (Mei & Lawrence Reference Mei and Lawrence1996), from laminar boundary-layer scaling. This $t^{-1/2}$ scaling can be related to the Basset–Boussinesq force (Basset Reference Basset1888; Boussinesq Reference Boussinesq1903), which accounts for the history effect in unsteady Stokes flow around a sphere. The relative acceleration between the body and the free stream influences the development of the boundary layer, which subsequently affects the drag force. This history effect is captured in the Basset–Boussinesq force,

(3.7) \begin{align} F_B(t) = 6 a^2 \sqrt {\pi \nu } \int _{-\infty }^t \frac {\ddot {X}(t^{\prime })}{\sqrt {t-t^{\prime }}}\mathrm{\,d} t^{\prime }, \end{align}

where $\ddot {X} (t^{\prime } )$ is the acceleration of the body at $t'$ . In the present case of an impulsively started flow over a sphere, the acceleration is a delta function at $t=0$ . Therefore, for $t\gt 0$ , the Basset–Boussinesq force is

(3.8) \begin{align} F_B(t) =6 a^2U\sqrt {\frac { \pi \nu }{t}}. \end{align}

This force decays proportionally to $t^{-1/2}$ , which directly corresponds to the initial decay observed in the total drag. Figure 9(c) shows the initial decay of the drag force for the impulsively started flow over a sphere, and the corresponding Basset–Boussinesq force $F_B$ . The $t^{-{1/2}}$ decay of the drag force agrees well with the decay of $F_B$ . Note that although at $t\to 0^{+}$ the thickness of the boundary layer over the sphere approaches zero, and most of the flow field is identical to potential flow, the total pressure drag does not approach zero. The divergence of the Lamb vector $\boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u}\times \boldsymbol{\omega })$ is the source term in the Poisson equation of the enthalpy function $h = ({p}/{\rho })+ ({1}/{2}) \vert \boldsymbol{u}\vert ^2$ . The thin singular vortex sheet on the sphere wall exerts a non-local effect on the pressure field that induces a non-zero form drag even at the initial time. The asymmetry of the pressure profile with respect to the polar angle at the initial time was verified by Dennis & Walker (Reference Dennis and Walker1972) and by Wang (Reference Wang1969), and was similarly attributed by Kang & Leal (Reference Kang and Leal1988) to the viscous generation of a vortex sheet on the surface for flow over a spherical bubble. The total drag from the J–A relation agrees well with the combination of form and friction drags. A distinction arises between the conventional division of drag into friction/form parts, and the J–A division into advection/diffusion parts at $t\to 0^+$ . The pressure and friction drags are both non-zero at the initial time; meanwhile, the J–A relation attributes all of the initial drag power injection to the viscous flux through the wall. The advection contribution is zero at $t \to 0$ , since the Lamb vector $\boldsymbol{u}\times \boldsymbol{\omega }$ is orthogonal to the potential-flow direction in the vortex sheet on the wall. At this initial time, the J–A interpretation of the drag force in terms of advective and diffusive fluxes of vorticity is more instructive compared to the conventional decomposition into shear stress and pressure. Specifically, diffusion of vorticity is the sole relevant actor at the initial time, while the shear stress and pressure are by-products of the action of the instantaneous diffusive flux. The viscous contribution drops rapidly during the initial development stage since the viscous flux reduces as the boundary layer grows. The advective contribution increases gradually from zero because of the detachment of the boundary layer and vortex formation on the back of the sphere. At later times, $t\gt 1$ , the viscous contribution in the J–A relation reduces, while the advection part increases and dominates the total rate of work by drag due to the detachment of vorticity from the near-wall region.

Figure 9.

( $a$ ) Time history of drag coefficient for case ST from integration of wall pressure and wall shear stress, and from the J–A relation. ( $b$ ) A zoomed-in view of $(a)$ during $0\leqslant t \leqslant 2$ .: total drag work evaluated by surface integration of pressure and friction.: pressure work.: friction work.: total drag work evaluated from the J–A relation. : total advective contribution $\int _{\Omega }\Pi _a\,\mathrm{d}V$ . : total viscous contribution $\int _{\Omega }\Pi _\nu\, \mathrm{d}V$ . $(c)$ A further zoomed-in view of $(b)$ in $0.08\leqslant t\leqslant 0.5$ in log scale, and including () the Basset–Boussinesq force $F_B$ .

Figure 10.

Contours of instantaneous vorticity fluxes for case ST at (i–iii) $t=0.3, 0.9, 1.5$ . The region enclosed in the black box is visualised in figure 11.

Figure 11.

Two-dimensional contours of instantaneous value of drag contribution and pressure coefficients for case ST. The top and bottom rows correspond to $t=0.9,1.5$ , respectively.

The spatial contributions to drag power injection $\Pi _a$ and $\Pi _\nu$ are visualised in figure 10 at three different times, $t=0.3,0.9,1.5$ . At the first instant, all of the vorticity in this flow is confined to the thin boundary layer on the sphere surface. Most of the drag originates from the viscous effect within the surface layer, where the vorticity transport pushes the flow away from ideal to vortical and produces drag force. A thin reverse-flow layer is already forming at $\varphi \gt 90^{\circ }$ due to the adverse pressure gradient. The locations $S_1$ and $R_1$ represent the primary separation and reattachment, which are lines of zero friction along the azimuthal direction. At $t=0.9$ , the total advection and viscous contributions become comparable. The advection part $\Pi _a$ becomes stronger due to the detachment of the primary boundary layer from the geometry and the formation of an axisymmetric primary vortex, as also shown in figure 11. The separation location $S_1$ moves upstream due to the adverse tangential pressure gradient. At $t=1.5$ , the primary vortex starts to detach from the solid wall, with a low-pressure region forming at the vortex core. The secondary boundary layer beneath the vortex starts to separate due to a pressure gradient that is adverse for the secondary layer but favourable for the primary flow. A secondary vortex forms beneath the primary one between $S_2$ and $R_2$ , which appear in a pair due to the topological constraints on singular points on the surface (Tobak & Peake Reference Tobak and Peake1982).

3.3. Turbulent flow over a sphere

After the starting stage discussed in § 3.2, the wake continues to extend downstream. The detached cylindrical shear layer breaks down to turbulence, and a statistically stationary state is established within the simulation domain after a sufficiently long time. The transient time from the initial condition to the stationary stage was approximately $150$ advective time units. Flow statistics were evaluated by averaging over $100$ time units and in the azimuthal direction. At this Reynolds number, although the boundary layer on the sphere is still laminar, the detached shear layer develops instability at approximately $x=2$ (Yun et al. Reference Yun, Kim and Choi2006; Rodriguez et al. Reference Rodriguez, Borell, Lehmkuhl, Segarra and Oliva2011) where localised azimuthal vortex roll-up starts to appear. At the end of the separation bubble, corrugated structures along the azimuthal direction develop, and ultimately, the wake breaks down to turbulence, which is observed in the instantaneous visualisation of streamwise velocity in figure 12 $(a)$ . The mean streamwise velocity field, as well as its comparison with previous DNS and experimental data, is shown in figures 12(b,c). We focus on the physical interpretation of the J–A relation and the vorticity dynamics in the presence of the turbulent wake.

Figure 12.

The streamwise velocity for case ST. (a) A snapshot of instantaneous $u$ velocity around the sphere during the statistically stationary stage. (b) Mean streamwise velocity $\overline {u}$ . (c) Comparison of the mean-velocity profiles from () the present simulations, with () previous DNS by Rodriguez et al. (Reference Rodriguez, Borell, Lehmkuhl, Segarra and Oliva2011), and () experimental data by Kim & Durbin (Reference Kim and Durbin1988). The $\overline {u}$ profiles are plotted along the radial direction in a cylindrical coordinate at streamwise locations  $x=0.2,1.6,3.0$ .

Figure 13.

(a) Time history of the drag coefficient for case ST from integration of the wall pressure and shear stress and from the J–A relation: () Total drag work evaluated by surface integration of pressure and friction. () time-averaged drag work evaluated by the J–A relation. () pressure work. () friction work. () total advective vorticity flux $\int _{\Omega }\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (-\boldsymbol{u}\times \boldsymbol{\omega })\,\mathrm{d}V$ . () total viscous vorticity flux $\int _{\Omega }\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\nu \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega })\,\mathrm{d}V$ . () oscillation of J–A drag. $(b)$ Symbols represent the time averaged J–A drag evaluated over a domain of radius $r$ . () time-averaged drag force from the summation of form and friction drag.

The contributions to the drag coefficient are presented in figure 13. The form drag is significantly higher than friction drag due to the higher Reynolds number considered relative to case SL. In addition, small temporal variations in the form and total drag can be discerned, which are due to the unsteadiness of the flow field. The similarity of values of form drag and advective flux contribution, as well as friction drag and viscous flux contribution, are observed in figure 13 as expected (the mathematical relations are provided in Appendix A). The advective term from the J–A relation exhibits high-frequency oscillations. The finite size of the simulation domain can accommodate only a portion of the full streamwise extent of the wake, here $15$ diameters. Vortices that leave the domain cause oscillations in the total advective vorticity flux and the rate of drag work from the J–A relation. The grey dotted line in figure 13(a) represents the spatial integration of the J–A relation outside the domain, which is computed by subtracting the J–A drag from the surface integral drag. The exchange of circulation between the computational domain and the outside flow causes the oscillation of J–A drag inside the computational domain. To examine the influence of the domain size on the accuracy of the J–A relation, we calculated the J–A integrals on domains $\Omega _r$ with different radii:

(3.9) \begin{align} \boldsymbol{F}_r \boldsymbol{\cdot} \boldsymbol{U}= \rho \int _{\Omega _r} \boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\boldsymbol{u}\times \boldsymbol{\omega }-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega })\,\mathrm{d}V. \end{align}

The mean value

(3.10) \begin{align} \overline {\boldsymbol{F}}_r \boldsymbol{\cdot} \boldsymbol{U}= \frac {1}{T_2-T_1}\int _{T_1}^{T_2} \boldsymbol{F}_r \boldsymbol{\cdot} \boldsymbol{U}\, \mathrm{d}t \end{align}

is visualised in figure 13(b), and converges rapidly with the increase of domain radius $r$ . The results indicate that the domain size is appropriate.

Figure 14.

Two-dimensional view of the instantaneous vorticity fluxes for case ST. Colour contours represent $(a)$ $\Pi _a+\Pi _\nu$ , ( $b.i$ , $c.i$ ) $\Pi _a$ , and ( $b.ii$ , $c.ii$ ) $\Pi _\nu$ . Lines correspond to the potential-flow streamlines with $\psi = 0, -0.1, -0.2, -0.3$ .

The instantaneous fields of $\Pi _a$ , $\Pi _{\nu }$ and their sum are visualised along with the potential-flow streamlines in figure 14. The advective term $\Pi _a$ dominates the total flux in most of the flow, while the viscous term $\Pi _{\nu }$ is important only near the wall. The advective flux $\Pi _a$ in this case exhibits behaviour similar to that observed in case SL in figure 7, where the circular detached shear layer near the sphere carries vorticity outwards crossing the potential streamlines, and produces drag. Farther downstream, the shear layer breaks down, carrying vorticity inwards towards the centreline. Crossing the potential streamline in this direction amounts to a force aligned with the potential flow, and thus reduces the rate of drag work. In this example, a direct analysis of the instantaneous vorticity transport is difficult due to the turbulence in the wake. We therefore proceed to analyse the mean vorticity transport, and exploit its symmetry.

Figure 15.

The mean vorticity fluxes cross the potential streamline for case ST. The solid black lines in the background represent potential function iso-surfaces $\psi =-0.1,-0.2,-0.3,-0.4$ .

The contributions to the mean flux tensor are simplified to the following form due to the axisymmetry of the flow statistics:

(3.11) \begin{align} \overline {\boldsymbol{\Sigma }}_{\overline {a}} = \begin{pmatrix} 0 & 0 & \overline {u}_x\,\overline {\omega }_{\varphi }\\[4pt] 0 & 0 & \overline {u}_r\,\overline {\omega }_{\varphi }\\[4pt] -\overline {u}_x\,\overline {\omega }_{\varphi } & -\overline {u}_r\,\overline {\omega }_{\varphi } & 0\\ \end{pmatrix},\quad \overline {\boldsymbol{\Sigma }}_{\nu } = \begin{pmatrix} 0 & 0 & -\nu \dfrac {\partial \overline {\omega }_{\varphi }}{\partial x}\\[8pt] 0 & 0 & -\nu \dfrac {\partial \overline {\omega }_{\varphi }}{\partial r}\\[4pt] \nu \dfrac {\partial \overline {\omega }_{\varphi }}{\partial x} & \nu \dfrac {\partial \overline {\omega }_{\varphi }}{\partial r} & 0 \end{pmatrix}, \end{align}

(3.12) \begin{align} \overline {\boldsymbol{\Sigma }}_{a'} = \begin{pmatrix} 0 & \overline {u_x'\omega _{r}'}-\overline {u_{r}'\omega _{x}'} & \overline {u_x'\omega _{\varphi }'}-\overline {u_{\varphi }'\omega _{x}'}\\[5pt] \overline {u_r'\omega _{x}'}-\overline {u_{x}'\omega _{r}'} & 0 & \overline {u_{r}'\omega _{\varphi }'}-\overline {u_{\varphi }'\omega _{r}'}\\[5pt] \overline {u_{\varphi }'\omega _{x}'}-\overline {u_{x}'\omega _{\varphi }'} & \overline {u_{\varphi }'\omega _{r}'}-\overline {u_{r}'\omega _{\varphi }'} & 0 \end{pmatrix}, \end{align}

where the overbars indicate mean quantities, and primes denote fluctuations. The flux of azimuthal vorticity is obtained by dotting the above flux tensors with the azimuthal unit vector $\boldsymbol{e}_{\varphi }$ :

(3.13) \begin{align} \overline {\boldsymbol{\Sigma }}_{\overline {a}} \boldsymbol{\cdot} \boldsymbol{e}_{\varphi } = \begin{pmatrix} \overline {u}_x \overline {\omega }_{\varphi }\\[4pt] \overline {u}_r \overline {\omega }_{\varphi } \\[4pt] 0\\ \end{pmatrix}, \quad \overline {\boldsymbol{\Sigma }}_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{\varphi } = \begin{pmatrix} -\nu \dfrac {\partial \overline {\omega }_{\varphi }}{\partial x} \\[8pt] -\nu \dfrac {\partial \overline {\omega }_{\varphi }}{\partial r} \\[4pt] 0\\ \end{pmatrix}, \quad \overline {\boldsymbol{\Sigma }}_{a'}\boldsymbol{\cdot} \boldsymbol{e}_{\varphi } = \begin{pmatrix} \overline {u_{x}'\omega _{\varphi }'}-\overline {u_{\varphi }'\omega _{x}'} \\[4pt] \overline {u_{x}'\omega _{\varphi }'}-\overline {u_{\varphi }'\omega _{r}'} \\[4pt] 0\\ \end{pmatrix}. \end{align}

The physical interpretation of the mean advection and diffusion fluxes of azimuthal vorticity is similar to the laminar case (§ 3.1). The turbulent flux $\overline {\boldsymbol{\Sigma }}_{a'}\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ is a combination of instantaneous advection and tilting effects. For example, the flux of azimuthal vorticity in the radial direction $\overline {u_r'\omega _{\varphi }'}-\overline {u_{\varphi }'\omega _{r}'}$ is comprised of (i) the transport of vorticity fluctuation $\omega _{\varphi }'$ by the fluctuating radial velocity $u_{r}'$ , and (ii) the tilting of fluctuating radial vorticity $\omega _{r}'$ into the azimuthal direction by the azimuthal fluctuation velocity $u_{\varphi }'$ .

With the detailed expressions for the vorticity fluxes, the mean J–A relation can be derived by ensemble averaging (2.12), which yields

(3.14) \begin{align} -\overline {\boldsymbol{F}} \boldsymbol{\cdot} \boldsymbol{U} &= -\rho \int _{\Omega }\boldsymbol{u}_{\phi }\boldsymbol{\cdot} \left (\overline {\boldsymbol{u}}\times \overline {\boldsymbol{\omega }}+\overline {\boldsymbol{u}'\times \boldsymbol{\omega }'}-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \overline {\boldsymbol{\omega }}\right )\mathrm{d}V \end{align}

(3.15) \begin{align}& = -\rho \int \mathrm{d}J\int \boldsymbol{e}_{n}\boldsymbol{\cdot} \left (\overline {\boldsymbol{\Sigma }}_{\overline {a}}+\overline {\boldsymbol{\Sigma }}_{a'}+\overline {\boldsymbol{\Sigma }}_{\nu }\right )\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }\, \mathrm{d}l. \end{align}

In this equation, the rate of work done by the mean drag is related to the spatial integration of mean advection $\overline {\Pi }_{\overline {a}} = \boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\overline {\boldsymbol{u}}\times \overline {\boldsymbol{\omega }} )=|\boldsymbol{u}_{\phi }|\,\boldsymbol{e}_{n} \boldsymbol{\cdot} \overline {\boldsymbol{\Sigma }}_{\overline {a}}\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ , turbulent advection $\overline {\Pi }_{a'} = \boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\overline {\boldsymbol{u}'\times \boldsymbol{\omega }'} )=|\boldsymbol{u}_{\phi }|\,\boldsymbol{e}_{n} \boldsymbol{\cdot} \overline {\boldsymbol{\Sigma }}_{a'}\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ , and diffusion $\overline {\Pi }_{\nu }= \boldsymbol{u}_{\phi }\boldsymbol{\cdot} (-\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \overline {\boldsymbol{\omega }} )=|\boldsymbol{u}_{\phi }|\,\boldsymbol{e}_{n} \boldsymbol{\cdot} \overline {\boldsymbol{\Sigma }}_\nu \boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ . The spatial distributions of these terms are shown in figure 15. The terms $\overline {\Pi }_{\overline {a}}$ , $\overline {\Pi }_{a'}$ and $\overline {\Pi }_{\nu }$ can also be interpreted as fluxes of mean azimuthal vorticity crossing iso-surfaces of the potential-flow streamfunction, weighted by the potential flow speed. The mean term $\overline {\Pi }_{\overline {a}}$ dominates the advective flux $\overline {\Pi }_{a}=\overline {\Pi }_{\overline {a}}+\overline {\Pi }_{a'}$ . Effectively, at $x\lt 1.5$ , the mean flow carries the mean azimuthal vorticity in the detached shear layer outwards, crossing the potential flow, which pushes the flow away from the ideal flow, and generates drag. Further downstream, the mean flow advects the mean azimuthal vorticity inwards towards the wake centreline, which shifts the flow back towards the potential flow, and corresponds to anti-drag. The turbulent flux $\overline {\Pi }_{a'}$ dominates near the wake centreline, at $x\gt 2$ . In this region, the turbulent flux accelerates the advection towards the wake centreline where the vorticity is annihilated. Note that the turbulent flux includes not only advection by the fluctuating velocity field, but also tilting of the fluctuating vorticity by the turbulent velocity.

Figure 16.

Balance between pressure gradient and vorticity fluxes. (a.i,a.ii) The flux lines tangent to $(\overline {\boldsymbol{\Sigma }}_a+\overline {\boldsymbol{\Sigma }}_\nu )\boldsymbol{\cdot} \boldsymbol{e}_{\varphi }$ (black lines with arrow) and iso-contours of the streamfunction $\psi$ (line and coloured contours) for cases SL and ST. (b.i,b.ii) The corresponding total pressure recovery, with lines representing total pressure difference: () $h(x)-h(R_1)$ ; () $T_{\nu }$ ; () $T_{a'}$ ; () $\pm (h_\infty -h_{R_1})$ .

The mean momentum equation relates the annihilation rate of mean azimuthal vorticity to the total pressure gradient, according to

(3.16) \begin{align} \left (\overline {\boldsymbol{\Sigma }}_{a'}\right )_{r\varphi }+\left (\overline {\boldsymbol{\Sigma }}_\nu \right )_{r\varphi }=\overline {u_r'\omega _{\varphi }'} - \overline {u_{\varphi }'\omega _{r}'}-\nu \frac {\partial \overline {\omega }_\varphi }{\partial r}=\frac {\partial \overline {h}}{\partial x}. \end{align}

Integrating (3.16) along the wake centreline from point $x'=R_1$ to a downstream location $x'=x$ yields a relation between the total pressure recovery and the vorticity flux along this line segment:

(3.17) \begin{align} \overline { T}_{a'} + \overline {T}_{\nu } = \int _{R_1}^{x}\left (\overline {\boldsymbol{\Sigma }}_{a'}\right )_{r\varphi }\mathrm{d}x'+\int _{R_1}^{x}\left (\overline {\boldsymbol{\Sigma }}_\nu \right )_{r\varphi }\mathrm{d}x'=\overline {h}(x)-\overline {h}(R_1). \end{align}

The flux lines of azimuthal vorticity and the recovery of total pressure along the wake centreline are both visualised in figure 16, for both SL and ST cases. Overall, the direction of the flux lines is similar in the laminar and turbulent cases. The azimuthal vorticity is generated at the wall, and annihilated along the wake centreline. However, the major mechanism of vorticity annihilation is different for the two cases. For case SL, the recovery of total pressure is due mainly to the viscous flux $ (\overline {\boldsymbol{\Sigma }}_{\nu } )_{r\varphi }$ , which is active at all $x\gt R_1$ and contributes a gradual recovery towards the free-stream condition. For case ST, the mean viscous flux is no longer important. Instead, the turbulent flux balances the mean total pressure gradient. Within the recirculation bubble, $R_1\lt x\lt 2$ , the turbulence is nearly absent, thus the total pressure remains nearly constant. Beyond this region, the detached shear layer starts to break down. The turbulence advection and tilting contribute to the flux of azimuthal vorticity towards the wake centreline. This physical interpretation of the vorticity transport reinforces our understanding of the mechanism of drag, by aid of the J–A relation and analysis of the Huggins flux tensor.