Vorticity dynamics and drag for flows over a sphere and a prolate spheroid

Abstract

The connection between the drag and vorticity dynamics for viscous flow over a bluff body is explored using the Josephson–Anderson (J–A) relation for classical fluids. The instantaneous rate of work on the fluid, associated with the drag force, is related to the vorticity flux across the streamlines of a background potential flow. The vorticity transport itself is examined by aid of the Huggins vorticity-flux tensor. The analysis is performed for three flows: flow over a sphere at Reynolds numbers $Re=200,3700$, and flow over a prolate spheroid at $Re=3000$ and $20^{\circ }$ incidence. In these flows, the vorticity transport shifts the flow away from and towards the ideal potential flow, with a net balance towards the former effect thus making an appreciable contribution to the drag. The J–A relation is first demonstrated for the flow over a sphere at $Re=200$. The drag power injection is related to the viscous flux of azimuthal vorticity from the wall into the fluid, and the advection of vorticity by the detached shear layer. In the wake, the azimuthal vorticity is advected towards the wake centreline and is annihilated by viscous effects, which contributes a reduction in drag. The analysis of the flow over a sphere at $Re=3700$ is reported for the impulsively started and stationary stages, with emphasis on the effects of unsteady two-dimensional separation and turbulent transport in the transitional wake. The turbulent flux in the wake enhances the transport of mean azimuthal vorticity towards the wake centreline, and is the main driver of the recovery of enthalpy between the rear point of the sphere and far downstream. The rate of work on the fluid by the drag force for a prolate spheroid is mostly due to the transport of vorticity along the separated boundary layers. Primary and secondary separation contribute oppositely to the power injection by the drag force, while the large-scale vortices only re-distribute vorticity without a net contribution. A mechanism for secondary separation is proposed based on the theory of vortex-induced separation.

1. Introduction

The relationship between drag and the transport of vorticity is established for some canonical flows; perhaps most known among them are the flows in channels and pipes (Taylor 1932). In these configurations, the transverse transport of vorticity leads to a drag force in the streamwise direction. Such interdependence provides a special perspective to interpret drag in viscous, vorticity-dominated and vortex-dominated flows. For flows over immersed bodies, this connection is given by the detailed Josephson–Anderson relation that equates the rate of work done by the drag force and the flux of vorticity against a background potential flow (Eyink 2021). In this study, we examine this relation for three-dimensional separated flows over spheres at moderate Reynolds number and over a spheroid at incidence.

The description of vorticity transport involves the notion of vorticity flux. In the conservative form of the vorticity equation, the Huggins flux tensor compactly captures the vorticity evolution by advection, tilting and stretching, and viscous diffusion (Huggins 1970, 1971). However, the vorticity equation dictates the form of the vorticity-flux tensor only up to a divergence-free contribution. A different form of the vorticity flux was introduced by Lighthill (1963) and Panton (1984) at the solid wall, and subsequently extended into the fluid interior (Kolár 2003). The physical interpretations of Huggins and Lighthill–Panton flux tensors were discussed and compared by Terrington, Hourigan & Thompson (2021). Both definitions similarly capture the spatial transport of vorticity, although the two forms contain different viscous contributions. The Huggins flux tensor measures the viscous transfer of circulation due to tangential acceleration in the fluid, while the viscous part of Lighthill–Panton flux includes only the terms that lead to the local change of vorticity. In the present study, the Huggins definition is adopted because it is uniquely related to the momentum and force balances. The turbulent part of the vorticity flux, which corresponds to Reynolds stress in the momentum equations, has been used to characterise the structure of near-wall turbulence (Klewicki 2013). Recently, the Huggins flux tensor has been extensively applied to examine the vorticity transport in flow over rotating and translating spheres (Terrington et al. 2021) and free-surface flows (Terrington, Hourigan & Thompson 2022a , b ). Here, we will quantify the vorticity motion for the flows over a sphere and a spheroid.

The connection between drag and vorticity fluxes has been the subject of influential works. Taylor (1932) introduced the correspondence between streamwise pressure gradient and the transverse transport of spanwise vorticity in shear flows. Lighthill (1963) quantified the generation of vorticity at solid walls, and established the balance between wall vorticity flux and the pressure gradient in channel and pipe flows. For the particular flows of interest in this work, namely flows over isolated bluff bodies, the drag force has been expressed as the integral of physical quantities within the flow interior. For example, Lighthill (1986) indirectly related the rate of work done by the drag force to the viscous dissipation of kinetic energy. Through a similar argument, Stone (1993) equated the rate of work done by the drag force over a steadily translating drop to the energy dissipation in the surrounding fluid and within the drop itself. The total dissipation is further decomposed into contributions from enstrophy and an interface vorticity term associated with surface curvature. Howe (1995) and Magnaudet (2011) expressed the drag force as an integral formula in terms of velocity and vorticity fields. In particular, Howe (1995) constructed potential flows using different boundary conditions, and expressed lift and torque using these potentials. The work by Eyink (2021), although derived through a different approach, yielded a similar expression that he termed the detailed Josephson–Anderson (J–A) relation. This expression relates the drag force to vorticity transport, and thus provides fluid dynamical insight. These ideas were demonstrated recently for the interpretation of drag in unsteady, separated laminar flow over a hill (Kumar & Eyink 2024). Here, our focus is on bluff-body flows. Using the J–A relation, we describe the rate of work done on the fluid by the drag force in terms of the vorticity flux crossing the potential-flow streamlines. This connection highlights important regions of the flow. In addition, the J–A relation can be interpreted in terms of an energy transfer between the vortical flow and the ideal, or potential, flow around the object. These ideas provide new perspectives on the role of vorticity and its dynamics in bluff-body flows and the generation of drag, which we explore in the present work for the flows around a sphere and a prolate spheroid.

The flow around a sphere exhibits a wealth of fluid dynamical phenomena, which have been the subject of experimental and numerical studies. The early works by Achenbach (1974) and Sakamoto & Haniu (1990) systematically documented the drag force and the flow states over a wide range of Reynolds numbers $400\leqslant Re\leqslant 5\times 10^{6}$ based on the sphere diameter. Direct numerical simulations (DNS) and large eddy simulations (LES) were performed in the subcritical range $Re\lt Re_c=3.7\times 10^{5}$ , where separation is laminar (Johnson & Patel 1999; Constantinescu & Squires 2004; Yun, Kim & Choi 2006; Rodriguez et al. 2011; Bazilevs et al. 2014). The work by Johnson & Patel (1999) documented several flow regimes, including steady axisymmetric flow ( $Re\leqslant 200$ ), steady non-axisymmetric flow ( $210\leqslant Re\leqslant 270$ ), and planar vortex shedding ( $Re=300$ ). Based on DNS at $Re=3700$ , Rodriguez et al. (2011) and Bazilevs et al. (2014) reported the flow statistics, morphology of the vortices, and the shedding mechanism. Vortices form and shed at random azimuthal angles due to the change in wall pressure, which results in a helical-shape wake. Yun et al. (2006) performed LES at $Re=3700, 10^{4}$ , and found that the higher $Re$ case exhibits a smaller recirculation region, and earlier transition and recovery of the wake. These studies establish a qualitative physical picture about the motion of vortices using visualisation of vortical structures, passive particle tracers, and frequency analysis of the shedding process. We will complement those efforts by providing a quantitative analysis, based on the Huggins flux tensor and directly linking the recovery of the wake to the vorticity transport.

Compared to the subcritical flow over a sphere, additional complexity is introduced when three-dimensional separation develops over the surface of a bluff body. An important example is the flow over a spheroid at incidence. In an early study, Wang et al. (1990) systematically investigated the friction lines and separation patterns on prolate spheroids with different aspect ratios and incidence angles, at subcritical Reynolds number. They discovered the phenomenon of open separation, and showed that multiple separations and reattachments can appear on the leeward side depending on the incidence angle. Ahn (1992) and Wetzel (1996) identified the critical Reynolds number $Re_c\sim 4.2\times 10^{5}$ based on the length of the minor axis of a 6 : 1 prolate spheroid, where transition occurs in the boundary layer prior to separation. They documented the flow statistics and friction patterns at both subcritical and supercritical conditions. Fu et al. (1994) studied the counter-rotating vortices and the crossflow separations above the leeward surface, and related the circulation in the large-scale vortices to the lift and lateral forces. Such a relation underscores the importance of vorticity dynamics in the near-body field. The DNS of flow over a spheroid have been limited to the subcritical flow regimes (El Khoury et al. 2010; Jiang et al. 2016), and have focused primarily on the vortical structures in the turbulent wake rather than of the separated boundary layer. Jiang et al. (2016) reported a non-symmetric helical wake behind a 6 : 1 prolate spheroid at $45^{\circ }$ incidence, and identified the axial separation line as the origin of the vortex generation. Ortiz-Tarin, Nidhan & Sarkar (2021) performed LES at $Re=10^{5}$ based on the length of the spheroid minor axis, and proposed a new decay rate for the velocity deficit based on the non-equilibrium dissipation scaling in the far wake. More recent LES by Plasseraud, Kumar & Mahesh (2023) demonstrated that tripping the boundary layer has a limited effect at $Re=7\times 10^{5}$ and $20^{\circ }$ incidence. For a theoretical treatment, Wu et al. (2000) and Surana et al. (2006) introduced criteria for the identification of three-dimensional separation, and discussed the relationship between vorticity and separation. Also relevant to the present effort are studies of vortex-induced separation (Peridier, Smith & Walker 1991; Doligalski, Smith & Walker 1994), since the dominant counter-rotating vortices on the leeward side of the spheroid can interact with the underlying boundary layer. Building on these previous efforts, we will directly evaluate the motion of axial vorticity using the flux tensor, and provide a vorticity-based mechanism for separation on the spheroid.

In § 2, we start by introducing the theoretical framework that is adopted for our study, including the set-up of the flow over a bluff body, the governing equations, and the J–A relation. The numerical simulations of the flows over a sphere and a prolate spheroid are presented in § 2.3. The J–A relation and the Huggins flux tensor are evaluated numerically. We first present the analysis of laminar flow over a sphere at $Re=200$ as a preliminary example in § 3.1. The more complex case of an impulsively started flow and the late-stage stationary state over the sphere at $Re=3700$ are analysed in §§ 3.2 and 3.3, with a focus on the two-dimensional unsteady separation and the turbulent wake dynamics. We then proceed to discuss the vorticity transport in the three-dimensional separation on the prolate spheroid in § 4. A summary and conclusion are provided in § 5.

2. Theoretical formulation and computational approach

In this section, we introduce the computational set-up for simulating the flow over bluff bodies, including the domain geometries, the governing equations and the boundary conditions. We discuss the vorticity dynamics, including the Helmholtz equation and the Huggins flux tensor. We then provide a brief derivation of the J–A relation, which is the theoretical approach that we adopt to interpret the power injection into the fluid in terms of vorticity fluxes.

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 (2021). 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 (1963) and Morton (1984). 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 1974), and $4.2\times 10^{5}$ for flow over a spheroid (Ahn 1992) – 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 1982). 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 1999). For an example of three-dimensional separation, consider flow over a prolate spheroid with a non-zero angle of attack (Wang et al. 1990). 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 1991; Surana, Grunberg & Haller 2006).

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. 2000). Second, the wall pressure gradient is related to the vorticity flux through (2.10) (Wu, Ma & Zhou 2007). 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 (2021). 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 2019; You & Zaki 2019). 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 2002). 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 1991). 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. Vorticity dynamics in flow over a sphere

In this section, the vorticity dynamics and its connection to the rate of work exerted on the fluid by the drag force for flow over a sphere are studied quantitatively using the J–A relation. Although the geometry is simple, the flow exhibits rich physical phenomena at different Reynolds numbers. We choose $Re=200$ (case SL) as the first example because the flow at this Reynolds number is sufficient to build a physical picture of vorticity transport. We then proceed to study the vorticity transport in the unsteady, impulsively started, turbulent wake (case ST) $Re=3700$ .

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. 2021):

(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 1996), from laminar boundary-layer scaling. This $t^{-1/2}$ scaling can be related to the Basset–Boussinesq force (Basset 1888; Boussinesq 1903), 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 (1972) and by Wang (1969), and was similarly attributed by Kang & Leal (1988) 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 1982).

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. 2006; Rodriguez et al. 2011) 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. (2011), and () experimental data by Kim & Durbin (1988). 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.

4. Three-dimensional boundary-layer separation over a spheroid

The vorticity transport and the interpretation of drag power injection over a bluff body become more complex in three-dimensional separation. In this section, we consider flow over a prolate spheroid at incidence angle $\alpha = 20^{\circ }$ . A visualisation of the near-wall flow is shown in figure 17. The free-stream flow encounters the spheroid body at the stagnation point near the nose on the windward side. A thin boundary layer forms on the forebody and extends downstream. On the windward side, the attached boundary layer grows along the axial direction. On the two sides of the body, the boundary layers separate due to the combined effects of pressure gradient and curvature. The detached shear layers curve inwards towards the recirculation region, and form a pair of counter-rotating vortices. These dominant vortices extend to form the turbulent far wake. The large-scale vortices induce a complex pattern of wall shear stress, visualised by the friction lines on the body surface in figure 17(a). A primary separation line ( $S_1$ ) appears where the boundary layer detaches from the body. The reverse flow induced by the dominant vortices near $\varphi = 180^{\circ }$ separates again at the secondary separation line $S_2$ . A reattachment line $R_2$ appears between $S_1$ and $S_2$ . These separation and reattachment lines are well identified by the strong convergence and divergence of the wall-friction lines. We first investigate the balance given by the J–A relation, and analyse the regions that contribute most to drag work, then proceed to analyse the vorticity transport and its relation to separation.

Figure 17. $(a)$ Visualisation of friction lines on the surface of the spheroid and the velocity field on selected vertical planes, $x/a=0.16,0.5,0.84$ . $(b)$ Time history of the drag coefficient from the surface integral and the J–A relation. () Total drag work evaluated by surface integration of pressure and friction. () pressure work. () friction work. () total advective vorticity flux $\int _{\Omega }\Pi _a\,\mathrm{d}V$ . () total viscous vorticity flux $\int _{\Omega }\Pi _\nu\, \mathrm{d}V$ . () oscillation of J–A drag.

The time histories of the contributions to drag power injection were evaluated from the surface forces and from the J–A relation, and are reported in figure 17. Although the Reynolds number is similar to turbulent flow over a sphere (case ST), the form drag here is of magnitude similar to that of the friction drag due to the slender shape of the body. Turning to the terms in the J–A relation, the advective contribution oscillates due to the unsteadiness of the velocity and vorticity fields, caused by the breakdown of the primary vortices.

Figure 18. Contours of the instantaneous values of $(a)$ $\Pi _a+\Pi _\nu$ , $(b)$ $\Pi _a$ and $(c)$ $\Pi _\nu$ for the flow over the spheroid. Overlaid on the contours are potential-flow streamlines. The vertical sections in the three rows are taken at $x/a=0.16,0.5,0.84$ . The labels $B_1$ and $B_2$ identify the primary and secondary boundary layers, $S_1$ and $S_2$ mark the primary and secondary separations, and $R_2$ is the secondary reattachment. The primary reattachment is on the leeward plane of symmetry, and is not marked on the figure.

The instantaneous spatial contributions to the rate of drag work, $\Pi _a(\boldsymbol{x})= -\boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\boldsymbol{u}\times \boldsymbol{\omega } )$ and $\Pi _{\nu }(\boldsymbol{x}) = \boldsymbol{u}_{\phi }\boldsymbol{\cdot} (\nu\, \boldsymbol{\boldsymbol{\nabla} }\times \boldsymbol{\omega } )$ , and their sum are visualised in figure 18 in vertical planes at the three locations marked in figure 17 $(a)$ , $x/a=0.16,0.5,0.84$ . The contours are overlapped with the potential-flow streamlines. On each vertical plane, the boundary layer starts from the stagnation point at $\varphi =0$ on the windward side, and separates at different azimuthal angles. The separated layer, marked $B_1$ in the figure, crosses potential streamlines carrying vorticity into the free stream. This advection of vorticity shifts the flow away from being ideal and thus contributes to drag power injection through $\Pi _a$ . The secondary boundary layer, marked $B_2$ in the figure, is induced on the leeward side by the large-scale vortices, and has a favourable influence of reducing drag due to the reversed sign of $\omega _x$ . Further away from the body, the pair of primary vortices advects vorticity across the potential-flow streamlines, transporting vorticity inwards at the top and outwards at the bottom, thus simultaneously contributing to and detracting from drag in these two regions, respectively. The viscous flux $\Pi _{\nu }(\boldsymbol{x})$ is mostly concentrated inside the attached primary and secondary boundary layers, where azimuthal vorticity is generated by streamwise pressure gradient and diffuses outwards.

Figure 19. Contours of the streamwise vorticity, overlaid by the vorticity-flux vectors. The lines are tangent to the vector fields due to the ( $a$ ) advective flux $\boldsymbol{\Sigma }_a\boldsymbol{\cdot} \boldsymbol{e}_{x}$ and ( $b$ ) viscous flux $\boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{x}$ . (i–iii) The three panels correspond to the vertical sections at $x/a=0.16, 0.5, 0.84$ .

For a detailed description of the transport of vorticity, we examine the Huggins flux tensor, specifically as it appears in the planar vorticity conservation equation (3.5). We stress that this description is Eulerian, thus it is focused on how the vorticity is transported locally. We consider $\omega _{x} = \boldsymbol{\omega }\boldsymbol{\cdot} \boldsymbol{e}_{x}$ since the streamwise vorticity is most relevant to the three-dimensional separation in the present flow. The vector $\boldsymbol{\Sigma }\boldsymbol{\cdot} \boldsymbol{e}_{x}$ represents the two-dimensional transport of $x$ vorticity inside the vertical plane. The advection $\boldsymbol{\Sigma }_a\boldsymbol{\cdot} \boldsymbol{e}_{x}$ and viscous $\boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{x}$ components are visualised in figure 19 at the same three axial locations examined in figures 17 and 18. Near the wall, the advection flux is parallel to the surface, and the viscous flux lines start from the wall. The streamwise vorticity is generated at the wall at low azimuthal angles by the azimuthal pressure gradient, and diffuses into the fluid. Within the detached primary boundary layer, vorticity is primarily advected. The dominant vortices are regions of accumulation of streamwise vorticity in the recirculation region, as shown in the planes $x/a=0.5,0.84$ in figure 19. A zoomed-in view in figure 20 shows one of the two counter-rotating vortices, which induces a secondary boundary layer (and similarly for the other vortex across the vertical symmetry plane). The induced boundary layer contains $\omega _x$ of the opposite sign to the vorticity in the large-scale vortex. The secondary boundary layer detaches again from the leeward side, which is known as the secondary separation (Wang et al. 1990; Wetzel 1996). The streamwise vorticity $\omega _x$ is advected along this separated layer, and diffuses towards the large-scale vortex and the surface. Note that a secondary vortex appears near the end of this secondary layer, shown by the streamlines in figure 20(c).

The wall-friction lines and vorticity-flux contours are plotted in figure 21. The locations of the separation and reattachment lines can be identified easily by the converging/diverging shear stress lines. The viscous flux is positive at lower azimuthal locations, and becomes negative upstream of the primary separation. Flows migrating from the upstream region carry excess vorticity that needs to be absorbed into the wall before separation, thus the wall vorticity flux changes sign prior to the primary separation line.

Figure 20. Visualisation of streamwise vorticity and velocity contours, overlaid by the vorticity-flux lines and the velocity vectors. The shown region is the boxed area in figure 19. The lines are tangent to $(a)$ $\boldsymbol{\Sigma }_{a}\boldsymbol{\cdot} \boldsymbol{e}_{x}$ , $(b)$ $\boldsymbol{\Sigma }_{\nu }\boldsymbol{\cdot} \boldsymbol{e}_{x}$ and $(c)$ $\boldsymbol{u}$ .

Figure 21. Visualisation of wall vorticity flux and friction lines. Colour indicates viscous flux $\Pi _{\nu }$ . Lines are friction lines that are tangent to wall shear stress $\boldsymbol{\tau }_w$ .

Figure 22. Schematic of vortex-induced separation, following the lecture notes on turbulence theory by G. Eyink. Positive and negative symbols represent the signs of the local streamwise vorticity, and colour indicates the magnitude. Positive values point into the page.

In the rest of this section, we discuss the mechanism of vortex-induced secondary separation. The interpretation is motivated by the local flow structures, and invokes ideas from two-dimensional separation. We stress, however, that the secondary separation on the spheroid is three-dimensional, and that the present discussion is intended to build intuition regarding the local flow near the secondary separation. Following Doligalski et al. (1994), consider an inviscid vortex with negative circulation above a solid wall, as shown schematically in figure 22(a). Since the wall is stationary, a vortex sheet is introduced at the surface, with strength $\gamma _w=\boldsymbol{n}\times \boldsymbol{w}(z)$ , where $\boldsymbol{w}(z)$ is the fluid velocity. The velocity and pressure satisfy the inviscid balance

(4.1) \begin{align} -\frac {\mathrm{d}p}{\mathrm{d}z}=w\frac {\mathrm{d}w}{\mathrm{d}z}. \end{align}

Since the singular vortex sheet is proportional to the velocity magnitude, this equation can be regarded as a balance between the destruction of the vortex sheet by pressure gradient and the advection of the sheet. If a small amount of viscosity is introduced into the flow, the singular vortex sheet becomes finite in thickness, as shown in figure 22(b). The vorticity magnitude reduces, thus the advection of vorticity becomes smaller than its destruction by the pressure gradient. A thin layer of negative vorticity thus appears near the wall. This region of reverse flow ejects the boundary layer into the free stream, and the detached boundary layer could break up and form a secondary vortex, as shown in figures 22(c,d).

We return to the flow field and the vorticity flux above the spheroid at $x/a=0.83$ , as visualised in figure 20. The dominant vortex with negative streamwise vorticity induces the secondary boundary layer on the leeward surface with positive vorticity (figure 20 $(b)$ ). This boundary layer separates due to an adverse pressure gradient imposed by the dominant vortex. Streamwise vorticity is advected along the secondary boundary layer and ejected to form a secondary vortex. Viscous diffusion of positive $\omega _x$ points from the secondary boundary layer towards the wall, where the outflux is related to the azimuthal pressure gradient. These observations are consistent with the vortex-induced separation mechanism described in figure 22.

5. Conclusions

In this work, we performed a numerical investigation of the vorticity dynamics and its relation to the drag force for the flow over a bluff body. We first invoked the Josephson–Anderson (J–A) relation that expresses the rate of work done by the drag force as the spatial integration of advective and viscous vorticity fluxed across the potential-flow streamlines. To understand the contribution to drag power injection by different flow regions, we investigated the Huggins vorticity-flux tensor. Three numerical simulations of flow over a bluff body were performed: flows over a sphere at $Re = 200, 3700$ , and flow over a prolate spheroid at $Re = 3000$ and $20^{\circ }$ incidence angle. For each of these flows, the balance of the J–A relation, the contributions to drag from the spatially distributed vorticity fluxes, and the flux vectors were computed and visualised. Our analysis addressed the most relevant features of each flow, including steady and unsteady two-dimensional separations, turbulent wake and three-dimensional separation.

For the flow over a sphere at $Re=200$ , the power injection of drag force calculated from the J–A relation and the integration of the surface stresses agree well. At the sphere surface, the azimuthal vorticity is diffused into the fluid interior by viscosity, which contributes to drag through the viscous term in the J–A relation. The vorticity in the boundary layer is first advected outwards crossing the potential-flow streamlines, thus contributing to drag, and is then advected inwards towards the wake centreline, thus contributing anti-drag. At the higher Reynolds number $Re=3700$ , we first find that the drag from the J–A relation agrees well with the integration of the surface stresses during the impulsively starting stage. The J–A relation attributes all of the drag at $t=0^+$ to the viscous vorticity flux across the wall, while the friction and pressure drag are both non-zero at that time. The development of multiple unsteady separations on the sphere wall is accompanied by the ejection of the primary and secondary vortices into the fluid interior, which contributes to the advection part of drag power injection. During the statistically stationary stage, the turbulent transport of vorticity in the near wake becomes important. The annihilation of vorticity along the wake centreline is balanced by the gradient of the total pressure, and is dominated by the turbulent vorticity flux. For the flow over a prolate spheroid at $Re=3000$ , the large-scale axial vortices induce multiple three-dimensional separation and reattachment lines over the leeward surface. The axial vorticity is generated in the windward boundary layer and advected into the vortices in the near wake. The primary and secondary separated boundary layers contribute opposite effects, namely an increase and a reduction, to the drag force. In addition, the separation of the secondary boundary layer was interpreted in terms of the theory of vortex-induced separation.

There are several possible avenues of future research that involve the J–A relation and vorticity dynamics, in the context of the flow over a bluff body. First, transitional boundary-layer flows can be studied (Zaki et al. 2010; Wang, Eyink & Zaki 2022). Across laminar-to-turbulence transition, the wall vorticity magnitude increases significantly, and an up-gradient turbulent flux of spanwise vorticity appears inside the boundary layer (Lighthill 1986). Both phenomena introduce more complex mechanisms of vorticity transport in the near-wall region that are of both theoretical and practical interest. Second, the J–A relation provides a new framework to examine drag-reduction strategies, from a vorticity dynamics point of view, including riblets (García-Mayoral & Jiménez 2011; Choi, Moin & Kim 1993) and superhydrophobic surfaces (Daniello, Waterhouse & Rothstein 2009; Jelly, Jung & Zaki 2014). In addition, instead of attempting to minimise the full drag term, drag-reduction strategies can be sought that target the reduction of the outward flux of vorticity. Third, other physical effects can be incorporated into the J–A relation, such as a spatiotemporal body force. For example, in viscoelastic flows (Terrapon et al. 2004; Li & Graham 2007; Hameduddin et al. 2018), the influence of the polymer stress can be included in the momentum equation as a distributed body force. These effects can be taken into consideration in the J–A relation by accounting for the body forces in the derivation. Finally, while the present study adopted an Eulerian perspective to examine the vorticity transport, future efforts should consider the Lagrangian evolution of vorticity (Xiang, Eyink & Zaki 2025) in these bluff-body flows.