Back-in-time analysis of vorticity in viscous separated flows over immersed bodies

Abstract

We present a back-in-time analysis for the origin of vorticity in viscous separated flows over immersed bodies, using the adjoint-vorticity framework recently introduced by Xiang et al. (2025 J. Fluid Mech. vol. 1011, A33. The solution of the adjoint-vorticity equations yields the volume density of mean deformation, which captures the stretching and tilting of the earlier vorticity that leads to the terminal value. The analysis also takes into account the boundary contributions of vorticity and its flux. Three examples are considered. Steady, axisymmetric separation in the flow over a sphere at Reynolds number $Re=200$ is shown to be established due to wall flux from both upstream and downstream of separation, the latter contribution being absent from the classical description by Lighthill. For unsteady separation at higher $Re=300$, the streamwise vorticity within the wake hairpin vortex is traced back, quantitatively, to the azimuthal vorticity on the sphere. The third configuration is the flow over a prolate spheroid at $Re=3000$. The null vorticity at three-dimensional separation originates from the cancellation of opposite interior contributions adjacent to the separation surface. The contribution from the downstream side migrates across the separation surface into the upstream region due to a tilting effect – a fundamental distinction between two- and three-dimensional separation. We also examine the detached vortical structures. The streamwise vorticity in the primary vortex originates from tilting of near-wall azimuthal vorticity, differing from Lighthill’s conjecture that the origin is streamwise near-wall vorticity that arises due to the reduced Coriolis force. Finally, a necklace vortex in the turbulent wake is traced back in time, and is shown to have contributions from the spheroid trailing-edge shed shear layer and the large-scale counter-rotating primary vortices.

1. Introduction

Vorticity is a striking feature of bluff body flows, and its transport is relevant to understanding the dynamics in the boundary layer around the body, separation, vortex shedding and the wake. In this work, we investigate the back-in-time origin of vorticity in flows over a sphere and a prolate spheroid. The origin of vorticity at the onset of the two-dimensional separation in a flow over a sphere is evaluated and contrasted with the classical theory by Lighthill (1963). Furthermore, the development of the three-dimensional separation in the flow over a prolate spheroid is analysed from a vorticity dynamics perspective and contrasted to two-dimensional separation.

1.1. Flow over an immersed body

Incompressible flows over immersed bodies involve a vast range of phenomena that vary depending on the Reynolds numbers, the body shape and the free-stream conditions. A comprehensive review is therefore not attempted, and instead we focus this section on key elements relevant to this work. We first review separation of the boundary layer around the body, and then discuss important physical observations of the dominant vortical structures in its wake.

Separation is classified into two-dimensional and three-dimensional due to their distinct definitions and features (Tobak & Peake 1981; Chapman & Yates 1991; Délery 2001). The rigorous definition of two-dimensional separation was first given by Prandtl (1905), stating that separation occurs where the wall shear stress vanishes, $\tau _w=0$ , and its tangential derivative is negative. Haller (2004) developed a first-principle framework for two-dimensional unsteady separation with an asymptotic mean. In that work, separation profiles are characterised as non-hyperbolic unstable manifolds in a Lagrangian frame. By considering the evolution of material line curvature, precursors of separation involving material spike formation and material backbone are identified and examined (Surana & Haller 2008; Serra et al. 2018, 2020; Klose, Jacobs & Serra 2020; Santhosh et al. 2023). In particular, Klose et al. (2020) identify asymptotic separation manifolds using backward-time Lyapunov exponents, with which the present work shares the back-in-time perspective. However, rather than being concerned with material lines, we focus on the back-in-time Lagrangian dynamics of vorticity. Steady separated flows around a two-dimensional wing or past a sphere fall into the category of two-dimensional separation. In these flows, the generation and advection of vorticity from the wall are key aspects of the separation dynamics (Lighthill 1963; Wu et al. 2000). An Eulerian description of vorticity transport is commonly adopted in the existing literature on boundary-layer separation. The relevance of the vorticity dynamics is twofold. Firstly, the wall shear-stress $\boldsymbol{\tau }_w$ is related to the wall vorticity $\boldsymbol{\omega }_w$ by the relation, $\boldsymbol{\tau }_w=\mu \boldsymbol{\omega }_w\times \hat {\boldsymbol{n}}$ , where $\mu$ is the dynamic viscosity, and $\hat {\boldsymbol{n}}$ is the wall-normal unit vector pointing towards the fluid. The zero-stress condition by Prandtl (1905) is therefore equivalent to a zero-vorticity condition. Secondly, separation is often induced by an adverse pressure gradient along the streamwise direction, which is equivalent to a flux of negative vorticity out of the fluid. According to Lighthill (1963), the wall-vorticity flux by favourable pressure gradient supplies the boundary-layer vorticity. Adverse pressure gradient extracts that vorticity from the boundary layer and ultimately causes separation. This theory by Lighthill (1963) is examined and verified in later works by Wu & Wu (1998) and Terrington, Hourigan & Thompson (2021) in flows over a wing and a sphere, respectively. However, Lighthill’s theory is qualitative since it does not deliver precise information about the detailed cancellation of vorticity sources and sinks that produce the zero vorticity at the separation location. Furthermore, the extension of Lighthill’s theory to three-dimensional separation is difficult, which underscores the need for a rigorous and quantitative framework to study the Lagrangian dynamics of vorticity in such flows.

The two-dimensional separated flows discussed above are degenerate since they contain only one component of vorticity, and the separation line is formed by zero-vorticity points aligned along the azimuthal or spanwise direction. In three-dimensional separation, such degeneracy no longer exists, and the topology of the wall shear-stress field serves as an indicator of type of separation (Legendre 1956; Lighthill 1963; Hunt et al. 1978). Lighthill (1963) recognised the convergence of neighbouring friction lines towards a limiting line as separation. By Poincaré–Bendixson theorem (Hunt et al. 1978), such friction lines must start and end at a zero-stress point, a periodic orbit or a connected set of zero-stress points. Different characterisations of three-dimensional separations were developed since then by Surana, Grunberg & Haller (2006) using dynamical system theory, and by Wu et al. (2000) using indicators constructed from vorticity. The notion of separation lines linking saddle points to stable spirals or nodes of the wall-friction vector field (Lighthill 1963) is challenged by the observation of open separation in flows over prolate spheroids (Wang et al. 1990). Open separation lines start and end at locations with a non-zero friction vector, which seems to contradict the argument by Lighthill (1963). Open separation appears with multiple friction lines suddenly converging into one, and the starting point of this one separation line becomes un-identifiable due to the non-uniqueness of the backward integration of friction lines. (Han & Patel 1979; Yates & Chapman 1992; Wetzel, Simpson & Chesnakas 1998; Surana et al. 2006). Although three-dimensional separation is tightly related to pressure gradient and the wall-vorticity flux, a description of the Lagrangian vorticity dynamics around separation has not been put forward. In the present work, we trace the origin for the null component of vorticity along the three-dimensional separation line back in time for flow over a spheroid, and we quantify the cancellations among volume and boundary contributions that result in this null component.

The separated boundary layers detach from the immersed bodies and form various vortical structures in their wake. The wake behind a sphere exhibits rich phenomenology depending on the Reynolds number, including axisymmetric steady separation ( $20\lt Re\lt 210$ ), planar steady symmetric wakes ( $210\lt Re\lt 270$ ), planar shedding of hairpin vortices ( $270\lt Re\lt 500$ ) and three-dimensional shedding ( $Re\gt 500$ ) (Sakamoto & Haniu 1990, 1995; Johnson & Patel 1999; Tomboulides & Orszag 2000). Although the vorticity generated along a sphere surface is predominantly azimuthal, strong streamwise vorticity appears in the legs of the shed hairpin vortices (Tomboulides & Orszag 2000). Johnson & Patel (1999) related vortex shedding in the laminar planar shedding regime to the pressure gradient along the axis of the toroidal vortices and the motion of velocity focus points. However, a physical explanation of the appearance of streamwise vorticity in the shed vortices is absent. For the flow over a prolate spheroid, the dynamics of vortical structures around the body depends on the incidence angle. Experimental measurements of vortical structures, wall shear stresses and velocity statistics were documented by Fu et al. (1994) and Chesnakas & Simpson (1996) for flow over a prolate spheroid at $\alpha =20^{\circ }$ incidence. A pair of counter-rotating vortices is observed above the spheroid leeward surface, which contains most of the circulation in the flow. The counter-rotating vortical structures are strongly related to the lift force on the body (Fu et al. 1994), however, the specific origin of the streamwise vorticity contained in these structures has not been determined quantitatively. In this work, we quantify the back-in-time contributions to the streamwise vorticity in the shed primary vortices in the flow over a sphere and a prolate spheroid.

1.2. Lagrangian vorticity dynamics

The transport of vorticity in inviscid flows has powerful geometric and Lagrangian interpretations. In the absence of viscosity, a vorticity vector translates and deforms according to the same dynamical equation as an infinitesimal material line segment (Cauchy 1815). A Cauchy formula expresses vorticity at any point in space and time as the vorticity in the same material element at an earlier time, multiplied by the deformation gradient along the Lagrangian trajectory. The origin of vorticity can be easily recovered in inviscid flows by following the Lagrangian trajectories backward in time. However, such back-in-time Lagrangian analysis cannot be applied directly to viscous flows, since the inviscid Cauchy formula does not hold in the presence of viscosity. The mathematical works by Constantin & Iyer (2008) and Constantin & Iyer (2011) generalise the ideal Cauchy formula for viscous flows by introducing stochasticity into the Lagrangian trajectories. This viscous extension of the Cauchy formula represents the vorticity at a certain spatio-temporal location as the expectation of contributions from different backward stochastic Lagrangian trajectories. As such, an important distinction is established between the vorticity deformation tensor and material deformation. The latter is no longer directly linked to stretching and tilting of vorticity since material deformation is deterministic while the Lagrangian dynamics of vorticity involves stochastic Brownian motion due to viscosity.

The presence of walls, which can act as sources of vorticity, requires careful treatment in the stochastic Cauchy formulation of the origin of vorticity. Eyink et al. (2020a ) developed a Dirichlet boundary condition, where particles stick to the wall when their back-in-time stochastic trajectories reach the surface and their contribution becomes frozen. Using this boundary condition, Eyink et al. (2020b ) analysed the history of vorticity in near-wall sweep and eject events in a turbulent channel flow. Wang et al. (2022a ) introduced a Neumann boundary condition, where particles are effectively reflected at the wall, and applied the stochastic Lagrangian approach to study the origin of the wall stress increase in a transitional boundary layer.

While the stochastic Cauchy formula incorporates viscous diffusion into the Lagrangian dynamics of vorticity exactly, the implementation requires Monte Carlo simulations for sampling of the stochastic trajectories, which is computationally demanding. The convergence of statistics for these trajectories is slow with respect to the number of samples. In the recent work by Xiang, Eyink & Zaki (2025), the contribution to the vorticity averaged over the stochastic trajectories was shown to be governed by a partial differential equation referred to as the adjoint-vorticity equation. Using this new approach, they directly traced back near-wall vorticity in turbulent channel flow to its origin, and compared the relative importance of the stretching and tilting of interior vorticity with the wall contribution. They also examined the origin of high wall shear-stress events. In the present work, we use this new approach to investigate the origin of zero vorticity for two- and three-dimensional separation, and quantify the contribution to complex vortical structures in flows over a sphere and a spheroid.

We present the methodology in § 2, including the flow set-up and governing equations, the background of the stochastic Lagrangian analysis and the adjoint-vorticity equations. In § 3, we describe the numerical method and the details of the simulations of the flows over a sphere and a prolate spheroid. The back-in-time origins of vorticity at the onset of two-dimensional separation and inside the shed vortices for the flow over a sphere are discussed in § 4. The adjoint-vorticity analysis of three-dimensional separation and of the dominant vortical structures in the flow over a prolate spheroid are discussed in § 5. Conclusions are presented in § 6.

Figure 1. Flow over a fixed body. The free-stream velocity is $\boldsymbol{U}$ . Iso-surfaces of the Q-criterion, $Q=0.1$ , show the abundance of vortical structures. Colour shows the magnitude of vorticity.

2. Methodology

2.1. Flows over immersed bodies

We consider flow over an immersed body, as shown by the example in figure 1. A solid body occupying space $B\subset \mathbb{R}^3$ is held fixed in a infinitely large domain filled with a viscous fluid. The free stream flows at a constant velocity $\boldsymbol{U}^{*}$ . The motion of the fluid is governed by the incompressible Navier–Stokes equations

(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} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} &= - \boldsymbol{\nabla }p + \frac {1}{Re} {\nabla} ^2 \boldsymbol{u}, \end{align}

which are non-dimensionalised using the free-stream speed and a characteristic length $L^{*}$ of the solid body. The latter will herein be either the diameter of a sphere or the length of the minor axis of a spheroid. The velocity vector is $\boldsymbol{u}=(u,v,w)$ , and $p$ is the pressure. The Reynolds number is defined by $Re = |\boldsymbol{U}^{*} |L^{*}/\nu ^{*}$ , where $\nu^{*}$ is the kinematic viscosity. The (2.1) can be solved after prescribing suitable initial and boundary conditions. The flow velocity is equal to zero at the wall and $\boldsymbol{U}$ in the free stream

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

2.2. The stochastic Cauchy invariant

The vorticity equation is obtained by taking the curl of the viscous momentum (2.2)

(2.4) \begin{equation} \frac {\partial \boldsymbol{\omega }}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega } = \boldsymbol{\omega } \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} + \frac {1}{Re} {\nabla} ^2 \boldsymbol{\omega }, \end{equation}

where $\boldsymbol{\omega }=\boldsymbol{\nabla }\times \boldsymbol{u}$ . Vorticity in the fluid domain can be introduced at the free-stream boundary or through flux at the solid walls. The latter is given by

(2.5) \begin{equation} \boldsymbol{\sigma } := -\nu \hat {\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega }=-\hat {\boldsymbol{n}}\times \boldsymbol{\nabla }p-\nu \boldsymbol{\nabla }\boldsymbol{\omega }\boldsymbol{\cdot }\hat {\boldsymbol{n}}, \end{equation}

where the wall-normal unit vector $\hat {\boldsymbol{n}}$ points into the fluid domain, and the last equality is obtained by evaluating the momentum equation at the wall. The inviscid version of (2.4) is time reversible, and can therefore be marched forward or backward in time. As such, any observation of vorticity can be evolved back in time to identify its origin. This relation is shown schematically in figure 2(a). In contrast, the viscous (2.4) is not time reversible due to the viscous term. In order to trace back the origin of an observed vorticity, a stochastic Lagrangian approach was introduced by Constantin & Iyer (2008, 2011).

Figure 2. Schematics of inviscid and viscous Lagrangian dynamics of vorticity. Panel $(a)$ : vorticity transport along deterministic Lagrangian trajectory in inviscid flow. Solid blue line represents the deterministic Lagrangian trajectory. Panel $(b)$ : vorticity transport along stochastic trajectories. Different backward stochastic trajectories (coloured with red, green and blue) sample vorticity at time $s$ at different spatial locations. The black arrow $\boldsymbol{\omega }(\boldsymbol{x}_{\kern-1pt f}, T)=\mathbb{E}[\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})\boldsymbol{\cdot }\boldsymbol{\omega }(\tilde {\boldsymbol{A}}_T^{s }(\boldsymbol{x}_{\kern-1pt f}), s)]$ is the average of contributions from trajectories. Panel $(c)$ : geometric interpretation of mean deformation gradient $\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})$ . The vorticity at $(\boldsymbol{x},s)$ is transported along different stochastic trajectories to $(\boldsymbol{x}_{\kern-1pt f},T)$ , and their average is $\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})\boldsymbol{\omega }(\boldsymbol{x},s)$ .

The stochastic Lagrangian approach is shown schematically in figure 2(b). We consider an observation of vorticity at the terminal time $T$ and position $\boldsymbol{x}_{\kern-1pt f}$ . Particles are released at $(\boldsymbol{x}_{\kern-1pt f},T)$ , and their back-in-time trajectories are stochastic. In effect, by computing the backward trajectories we are efficiently identifying the original labels of particles that arrive at $(\boldsymbol{x}_{\kern-1pt f},T)$ . The back-to-label map $\tilde {\boldsymbol{A}}_t^s$ satisfies the equation

(2.6) \begin{equation} \hat {\textrm{d}} \tilde {\boldsymbol{A}}_T^s(\boldsymbol{x}_{\kern-1pt f})=\boldsymbol{u}\Big (\tilde {\boldsymbol{A}}_T^s(\boldsymbol{x}_{\kern-1pt f}), s\Big ) \textrm{d} s+\sqrt {2 v} \,\hat {\textrm{d}} \tilde {\boldsymbol{W}}(s), \quad s\lt T ; \quad \tilde {\boldsymbol{A}}_T^T(\boldsymbol{x}_{\kern-1pt f})=\boldsymbol{x}_{\kern-1pt f}. \end{equation}

The symbol $\hat {\textrm{d}}$ represents the backward Itô differential. The motion of these stochastic particles in (2.6) is driven by not only the background velocity field $\boldsymbol{u}$ , but also a vector backward Brownian motion $\tilde {\boldsymbol{W}}$ , where each realisation of $\tilde {\boldsymbol{W}}$ produces an individual backward Lagrangian trajectory. We denote $\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})$ as the deformation gradient along the stochastic trajectory. The vorticity $\boldsymbol{\omega }(\boldsymbol{x}_{\kern-1pt f}, T)$ is then given by the expectation of contributions from all sample trajectories

(2.7) \begin{equation} \boldsymbol{\omega }(\boldsymbol{x}_{\kern-1pt f}, T)=\mathbb{E}\Big [\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})\boldsymbol{\cdot }\boldsymbol{\omega }\Big (\tilde {\boldsymbol{A}}_T^{s }(\boldsymbol{x}_{\kern-1pt f}), s\Big )\Big ], \,\,\,\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})=\Big (\boldsymbol{\nabla }_x \tilde {\boldsymbol{A}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})\Big )^{-\top },\quad s\lt T . \end{equation}

Equation (2.7) was first derived by Constantin & Iyer (2008) and is explained with aid of figure 2 $(b)$ . The stochastic particles sample the vorticity at earlier time $s$ at different spatial locations $\tilde {\boldsymbol{A}}_T^{s}$ . The sampled vorticity $\boldsymbol{\omega } (\tilde {\boldsymbol{A}}_T^{s }(\boldsymbol{x}_{\kern-1pt f}), s )$ is tilted and stretched by the deformation matrix, $\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})\boldsymbol{\cdot }\boldsymbol{\omega } (\tilde {\boldsymbol{A}}_T^{s }(\boldsymbol{x}_{\kern-1pt f}), s )$ , along each stochastic trajectory in the same manner as an infinitesimal material line element. The terminal vorticity $\boldsymbol{\omega }(\boldsymbol{x}_{\kern-1pt f}, T)$ is then the expectation of the deformed vorticity over all possible trajectories. The stochastic formulation (2.6)–(2.7) is the viscous extension of the inviscid theory.

Note that, at an early time $s\lt T$ , the location of stochastic particle $\tilde {\boldsymbol{A}}_T^{s}$ could hit the wall. The interaction between the backward stochastic particles and the wall is specified by either a stick (Dirichlet) or reflection (Neumann) boundary condition. In the former, the stochastic particles adhere to the wall and in the latter they bounce back. The details of the boundary conditions are provided in Eyink et al. (2020a ) and Wang et al. (2022a ).

The stochastic Lagrangian analysis requires sampling of a large number of particle trajectories, which is computationally demanding. The sample mean converges to the true expectation in (2.7) with an error proportional to $\sigma /\sqrt {N}$ , where $N$ is the number sampled stochastic trajectories and $\sigma$ is the variance of $\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})\boldsymbol{\cdot }\boldsymbol{\omega } (\tilde {\boldsymbol{A}}_T^{s }(\boldsymbol{x}_{\kern-1pt f}), s )$ . The latter grows exponentially in backward time, and therefore it becomes prohibitively expensive to evaluate the expectation using the Monte Carlo approach. Xiang et al. (2025) recently introduced an adjoint-based approach to evaluate the back-in-time origin of vorticity, by solving the dual to the vorticity equation. The adjoint formulation is equivalent to the stochastic Lagrangian analysis and is computationally more efficient.

2.3. Adjoint-vorticity equation

Adjoint variational methods are known for their role in data assimilation (for a recent review, see Zaki 2025). The versatility of these methods has been exploited to answer fundamental questions regarding the accuracy of reconstructing turbulence from under-resolved observations (Wang & Zaki 2021) and to characterise measurements in wall-bounded turbulence (Wang et al. 2022b ). Rather than start from the Navier–Stokes equations, the present analysis focuses on the vorticity equation. The derivation follows the work by Xiang et al. (2025), and is the formulation that will be adopted throughout this work to analyse vorticity and separation in flows over immersed bodies.

Consider $\boldsymbol{f}$ and $\boldsymbol{g}$ as two time-dependent vector fields. An inner product $ ( \boldsymbol{f}, \boldsymbol{g} )$ is defined as

(2.8) \begin{equation} \left ( \boldsymbol{f}, \boldsymbol{g}\right ) := \int _s^T \int _{\mathcal{D}} \boldsymbol{f}^{\top } \boldsymbol{g} \textrm{d}V \textrm{d}t, \end{equation}

where the temporal and spatial integrals are performed over the interval $[s,T]$ and the fluid domain $\mathcal{D}=\mathbb{R}^3\setminus B$ . Taking the inner product between an adjoint vector field $\boldsymbol{\varOmega }$ and the forward vorticity (2.4) yields

(2.9) \begin{align} & \left (\boldsymbol{\varOmega }, \frac {\partial \boldsymbol{\omega }}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega }-\boldsymbol{\omega } \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}-\nu {\nabla} ^2 \boldsymbol{\omega }\right ) =\left (\boldsymbol{\omega }, \frac {\partial \boldsymbol{\varOmega }}{\partial \tau }-\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\varOmega }-\boldsymbol{\nabla } \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\varOmega }-\nu {\nabla} ^2 \boldsymbol{\varOmega }\right ) \nonumber \\& \qquad + \int _{\mathcal{D}} \boldsymbol{\varOmega }\left (\boldsymbol{x},T\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},T\right ) \textrm{d}V -\int _{\mathcal{D}} \boldsymbol{\varOmega }\left (\boldsymbol{x},s\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V \nonumber\\& \qquad + \int _s^T \oint _{\partial \mathcal{D}} \left (\boldsymbol{\varOmega } \boldsymbol{\cdot }\boldsymbol{\omega }\right ) \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n} + \nu \left (\left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\varOmega }\right ) \boldsymbol{\cdot }\boldsymbol{\omega }- \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right ) \boldsymbol{\cdot }\boldsymbol{\varOmega } \right ) \textrm{d}S \textrm{d}t, \end{align}

where $\tau =T-t$ is the backward time and the unit-normal vector $\boldsymbol{n}$ points outward from the fluid. The above expression introduces an important assumption. Unlike the four-dimensional variational approach to the Navier–Stokes equations (e.g. Zaki & Wang 2021), we do not seek an adjoint to the forward velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ which is instead held frozen. This assumption is essential for the derivation of the particular adjoint-vorticity equation that is equivalent to the stochastic Lagrangian formulation, as we will show in § 2.4. Since the forward vorticity equation is satisfied pointwise in the flow field, (2.9) is equal to zero. The first term on the right-hand side vanishes if the adjoint vorticity $\boldsymbol{\varOmega }$ is chosen to satisfy

(2.10) \begin{equation} \frac {\partial \boldsymbol{\varOmega }}{\partial \tau } - \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varOmega } - \boldsymbol{\nabla }\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\varOmega } -\nu {\nabla} ^2 \boldsymbol{\varOmega } = 0, \end{equation}

or in index notation

(2.11) \begin{equation} \frac {\partial \varOmega _i}{\partial \tau } - u_j\frac {\partial \varOmega _i}{\partial x_j} - \frac {\partial u_j}{\partial x_i} \varOmega _j -\nu \frac {\partial \varOmega _i}{\partial x_j \partial x_j} = 0. \end{equation}

Equation (2.10) can be solved backward in time once supplemented with suitable boundary and backward initial conditions. The boundary term $ (\boldsymbol{\varOmega } \boldsymbol{\cdot }\boldsymbol{\omega } ) \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}$ in (2.9) vanishes in absence of free-stream vorticity, at no-slip surfaces and at periodic boundaries – the first two conditions being relevant to the flows that we consider. The remaining terms in (2.9) define a forward-adjoint duality relation

(2.12) \begin{align} \int _{\mathcal{D}} \boldsymbol{\varOmega }\left (\boldsymbol{x},T\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},T\right ) \textrm{d}V & = \int _{\mathcal{D}} \boldsymbol{\varOmega }\left (\boldsymbol{x},s\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V \nonumber \\&\quad + \int _s^T \oint _{\partial \mathcal{D}} \nu \left ( - \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\varOmega }\right ) \boldsymbol{\cdot }\boldsymbol{\omega }+ \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right ) \boldsymbol{\cdot }\boldsymbol{\varOmega } \right )\,\textrm{d}S \textrm{d}t. \end{align}

The interpretation of (2.12) depends on the initial condition for the adjoint field $\boldsymbol{\varOmega } (\boldsymbol{x},T )$ . Consider, for example, a vector impulse located at $\boldsymbol{x}_{\kern-1pt f}$ pointing in the $k^{ {th}}$ coordinate direction $\boldsymbol{e}_k$

(2.13) \begin{equation} \boldsymbol{\varOmega }^k\left (\boldsymbol{x},T\right ) = \boldsymbol{e}_k \delta (\boldsymbol{x}-\boldsymbol{x}_{\kern-1pt f}). \end{equation}

The corresponding forward-adjoint duality becomes

(2.14) \begin{align} \omega _k\left (\boldsymbol{x}_{\kern-1pt f},T \right ) & = \int _{\mathcal{D}}\boldsymbol{\varOmega }^k\left (\boldsymbol{x},s\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V \nonumber \\& \quad - \int _s^T \oint _{\partial \mathcal{D}} \nu \big ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\varOmega }^k\big ) \boldsymbol{\cdot }\boldsymbol{\omega } \,\textrm{d}S \textrm{d}t + \int _s^T \oint _{\partial \mathcal{D}} \nu \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right )\boldsymbol{\cdot }\boldsymbol{\varOmega }^k \,\textrm{d}S \textrm{d}t \nonumber \\& = \mathcal{I}(s) +\mathcal{B}_D(s) + \mathcal{B}_N(s) . \end{align}

The superscript in $\boldsymbol{\varOmega }^k$ denotes the vector direction of the adjoint initial condition. Equation (2.14) expresses the target vorticity component $\omega _k (\boldsymbol{x}_{\kern-1pt f},T )$ in terms of an interior contribution $\mathcal{I}$ at earlier time $s$ plus two cumulative wall contributions, $\mathcal{B}_D(s)$ and $\mathcal{B}_N(s)$ . The interior contribution $\mathcal{I}$ involves the spatial integration of an inner product between the adjoint and forward vorticity fields at time $s$ . This inner product $\boldsymbol{\varOmega }^k (\boldsymbol{x},s ) \boldsymbol{\cdot }\boldsymbol{\omega } (\boldsymbol{x},s )$ thus accounts for the stretching and tilting of the vorticity at a certain spatio-temporal location $ (\boldsymbol{x},s )$ which contributes to the target vorticity component $\omega _k (\boldsymbol{x}_{\kern-1pt f},T )$ . The boundary terms, $\mathcal{B}_D(s)$ and $\mathcal{B}_N(s)$ , are assigned subscripts ‘D’ and ‘N’ since they involve the wall-vorticity values and fluxes, respectively. These terms are integrations over the solid boundaries and the time interval $[s,T]$ . The integrands $B_D=\nu ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\varOmega }^k ) \boldsymbol{\cdot }\boldsymbol{\omega }$ and $B_N = \nu \boldsymbol{\varOmega }^k \boldsymbol{\cdot }( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } )$ represent the instantaneous distributed contribution from solid walls to the target vorticity. A homogeneous Dirichlet boundary condition $\boldsymbol{\varOmega }^k=0$ fixes $\mathcal{B}_N(s)=0$ such that only the Dirichlet boundary term $\mathcal{B}_D(s)$ is sampled during the backward evolution. A Neumann boundary condition $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varOmega }^k=0$ , on the other hand, fixes $\mathcal{B}_D(s)=0$ and evaluates the contribution of the Lighthill wall-vorticity flux $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega }$ through $\mathcal{B}_N(s)$ . Based on (2.14), the full vorticity vector $\boldsymbol{\omega } (\boldsymbol{x}_{\kern-1pt f},T )$ can be expressed as

(2.15) \begin{align} \boldsymbol{\omega }\left (\boldsymbol{x}_{\kern-1pt f},T \right ) & = \int _{\mathcal{D}}{\unicode{x1D76E}}\left (\boldsymbol{x},s\right ) \boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V \nonumber \\& \quad - \int _s^T \oint _{\partial \mathcal{D}} \nu \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } {\unicode{x1D76E}}\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\textrm{d}S \textrm{d}t + \int _s^T \oint _{\partial \mathcal{D}} \nu \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right )\boldsymbol{\cdot }{\unicode{x1D76E}}\kern1pt \textrm{d}S \textrm{d}t, \end{align}

where $\unicode{x1D76E}=[{\boldsymbol{\varOmega }^{1}, \boldsymbol{\varOmega }^{2}, \boldsymbol{\varOmega }^{3}}]^{\top }$ . The operator $\unicode{x1D76E}$ is formed by solving the adjoint vorticity (2.10) using three initial impulses each in different coordinate direction, then combining these adjoint solution vectors into a matrix. In (2.15), the vorticity vector $\boldsymbol{\omega } (\boldsymbol{x}_{\kern-1pt f},T )$ is expressed in terms of the deformation of earlier vorticity in the domain interior and a time integral of wall contributions. The two boundary terms in (2.15) are again the wall vorticity sampled by the adjoint-vorticity flux, and the Lighthill wall-vorticity flux sampled by the adjoint.

Figure 3. Schematics of $(a)$ stretching and $(b)$ tilting of vorticity and adjoint vorticity. In inviscid flows, vorticity (blue vector) satisfies the equation of an infinitesimal material line element, and the adjoint vorticity (red area and area vector) satisfies the equation of an area element. $(a)$ Forward stretching of $\omega _y$ and compression of $\varOmega _y$ , and the back-in-time opposite behaviour. $(b)$ Forward tilting of $\omega _y$ and $\varOmega _y$ , and their back-in-time counterparts. The faint colouring indicates that the forward-time evolution of the adjoint vorticity and the backward-time evolution of vorticity are absent in viscous flows.

An intuitive demonstration of the kinematics of $\boldsymbol{\varOmega }$ can be provided in the inviscid limit, and complements the well-established description of the kinematics of $\boldsymbol{\omega }$ . In the inviscid limit, time reversibility permits both forward and backward evolutions of vorticity and of the adjoint. Figure 3 shows these evolutions in ( $a$ ) simple straining and ( $b$ ) rotational flows. The forward evolution of vorticity is governed by the same equation as an infinitesimal material line element, $\pm {\textrm{d}\boldsymbol{\omega }}/{\textrm{d}t} = \pm \boldsymbol{\omega } \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}$ , where $\textrm{d}\boldsymbol{\omega }/{\textrm{d}t}$ is the time derivative along a Lagrangian trajectory. Stretching in forward time becomes compression in backward time (figure 3 $a$ ). The forward evolution of the adjoint-vorticity variable satisfies the same equation for an area element, $\pm {\textrm{d}\boldsymbol{\varOmega }}/{\textrm{d}t} = \mp \boldsymbol{\nabla }\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\varOmega }$ . The straining flow reduces the area in forward time, and expands it in backward time. In this description, the analogue to incompressibility is the conservation of $\boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{\varOmega }$ , which is the inviscid form of the duality relation (2.12). In the rotational flow (figure 3 b), forward and adjoint vorticity initially aligned with the $y$ -direction are tilted into the $x$ -direction by $\omega _y ( {\partial u}/{\partial y})$ and $-\varOmega _y ( {\partial v}/{\partial x})$ , respectively. In backward time, these actions reverse sign. In viscous flows, only the forward evolution of vorticity and the backward evolution of the adjoint variable are physically meaningful, due to the time irreversibility of diffusion. The faint colouring in the figure thus indicates that, in the presence of viscosity, the adjoint variable cannot be evolved forward in time, just as the vorticity cannot be evolved backward in time.

2.4. Equivalence between the adjoint and stochastic representations

The adjoint vorticity (2.15) is equivalent to the stochastic representation (2.7). The solution $\tilde {\boldsymbol{A}}_t^s(\boldsymbol{x}_{\kern-1pt f})$ of the stochastic differential (2.6) is a stochastic process, and its probability density function $\varrho (\boldsymbol{x},\tau )$ is the solution of a backward Fokker–Planck equation with an impulse initial condition at $\boldsymbol{x}_{\kern-1pt f}$

(2.16) \begin{equation} \frac {\partial \varrho }{\partial \tau } - \boldsymbol{\nabla }\boldsymbol{\cdot }\left (\varrho \boldsymbol{u}\right ) = \nu {\nabla} ^2 \varrho , \qquad \varrho (\boldsymbol{x},T)=\delta (\boldsymbol{x}-\boldsymbol{x}_{\kern-1pt f}). \end{equation}

With the probability density $\varrho$ available, the stochastic representation (2.7) is re-written as

(2.17) \begin{eqnarray} \boldsymbol{\omega }(\boldsymbol{x}_{\kern-1pt f}, T)&=&\mathbb{E}\left [\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f}) \boldsymbol{\cdot }\boldsymbol{\omega }\Big (\tilde {\boldsymbol{A}}_{T}^{s}(\boldsymbol{x}_{\kern-1pt f}), s\Big )\right ] \nonumber \\&=& \mathbb{E}\left [\mathbb{E}\Big [\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f}) \left .\right \vert \tilde {\boldsymbol{A}}_{T}^{s}(\boldsymbol{x}_{\kern-1pt f})\Big ]\boldsymbol{\cdot }\boldsymbol{\omega }\Big (\tilde {\boldsymbol{A}}_{T}^{s}(\boldsymbol{x}_{\kern-1pt f}), s\Big )\right ]\nonumber \\&=&\int _{\mathcal{D}}\mathbb{E}\Big [\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f}) \left .\right \vert \tilde {\boldsymbol{A}}_{T}^{s}(\boldsymbol{x}_{\kern-1pt f})=\boldsymbol{x}\Big ]\boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x}, s\right )\varrho (\boldsymbol{x},s)\textrm{d}\boldsymbol{x}. \end{eqnarray}

On the second line of (2.17), the expectation over all possible stochastic trajectories is decomposed into two nested expectations by the tower property of conditional expectations. The inner expectation is taken conditioned on a fixed backward terminal location of particle $\tilde {\boldsymbol{A}}_{t}^{s}(\boldsymbol{x})$ , and the outer expectation is taken over all possible $\tilde {\boldsymbol{A}}_{t}^{s}(\boldsymbol{x})$ . The third line of (2.17) re-writes the outer expectation as a spatial integral, and is in an identical form as the interior term in (2.15) if we take the adjoint-vorticity operator $\unicode{x1D76E}$ to be

(2.18) \begin{equation} \unicode{x1D76E}\kern1pt(\boldsymbol{x},s)=\varrho (\boldsymbol{x},s)\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x}) \quad \text{with} \qquad \mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})=\mathbb{E}\Big [\tilde {\boldsymbol{D}}_T^{s}(\boldsymbol{x}_{\kern-1pt f})\left .\right \vert \tilde {\boldsymbol{A}}_{T}^{s}(\boldsymbol{x}_{\kern-1pt f})=\boldsymbol{x}\Big ]. \end{equation}

The mean conditional deformation matrix $\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})$ describes how vorticity is rotated and stretched on average from $(\boldsymbol{x},s)$ to $(\boldsymbol{x}_{\kern-1pt f},T)$ . There are infinitely many possible stochastic trajectories connecting $(\boldsymbol{x},s)$ to $(\boldsymbol{x}_{\kern-1pt f},T)$ , and $\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})\boldsymbol{\omega }(\boldsymbol{x},s)$ represents the contribution of point $(\boldsymbol{x},s)$ to the final vorticity at $(\boldsymbol{x}_{\kern-1pt f},T)$ averaged over all possible trajectories (see figure 2 c). This term alone does not, however, take into account the probability of those trajectories being realised from different $(\boldsymbol{x},s)$ locations. This point requires multiplication by $\varrho (\boldsymbol{x},s)$ . The total vorticity $\boldsymbol{\omega }(\boldsymbol{x}_{\kern-1pt f}, T)$ is therefore the probability weighted superposition of $\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})\boldsymbol{\omega }(\boldsymbol{x},s)$ from all spatial locations $\boldsymbol{x}$ at time $s$ , as given by (2.17). The decomposition $\unicode{x1D76E}\kern1pt(\boldsymbol{x},s)=\varrho (\boldsymbol{x},s)\mathcal{S}_T^s(\boldsymbol{x}_{\kern-1pt f},\boldsymbol{x})$ is a new contribution that was not discussed in the paper by Xiang et al. (2025).

The above equivalence of the stochastic Lagrangian and the adjoint vorticity approaches did not consider the presence of a solid wall. A detailed proof for the equivalence with wall effects is presented in Xiang et al. (2025). Here, we simply remark that the Dirichlet and Neumann boundary conditions in the adjoint formulation are mathematically equivalent to the ‘stick’ and ‘reflection’ boundary conditions in the stochastic Lagrangian approach. Specifically, the Dirichlet condition corresponds to particles in the stochastic formulation adhering to the wall upon contact, and having their contribution frozen. The Neumann condition corresponds to the stochastic particles reflecting from the wall, in search of the primordial origin of the vorticity.

3. Numerical simulations of flows over immersed bodies

The adjoint-vorticity framework introduced in § 2 will be adopted to analyse the origin of vorticity at points of interest in the flow around an immersed body. The numerical solution of the adjoint (2.10) in backward time requires access to spatially and temporally resolved forward velocity fields throughout the time horizon. The forward flow fields are obtained by direct numerical simulations (DNS) of the incompressible Navier–Stokes equations (2.1)–(2.2) over the bodies. The equations are discretised using a fractional-step approach with a local finite volume flux formulation on a staggered curvilinear grid (Rosenfeld, Kwak & Vinokur 1991; Wang, Wang & Zaki 2019). The advection and diffusion terms are discretised in time using Adam–Bashforth and Crank–Nicolson schemes, respectively. To enforce incompressibility, the pressure Poisson equation is solved using stabilised bi-conjugate gradient method with an algebraic multigrid preconditioner from Hypre (Falgout & Yang 2002). The flow domain is decomposed into six structured blocks on which the above single-block numerical algorithms are adopted (see figure 4). An overlapping grid method is applied for communication of the velocity data between adjacent blocks. The pressure field is solved globally across all six blocks, due to the ellipticity of pressure Poisson equation.

Figure 4. Schematics of flow set-up, domain and grid for flow over sphere (top row) and prolate spheroid (bottom row). Left column: domain set-up and inflow. Middle column: grid number marked on a multi-block domain. Right column: the blue-coloured block in middle column rotated vertically. The sizes of the bodies in the left column are exaggerated for the purpose of a clear visualisation.

Two types of immersed bodies are considered, namely a sphere and a prolate spheroid with six-to-one aspect ratio (see table 1). The flow set-ups and grid systems are visualised in figure 4. For the flow over a sphere, the free-stream velocity is $\boldsymbol{U} = (1,0,0)$ . The flow equations are solved between two concentric spheres with radii $R_1=0.5$ and $R_2=15$ . This domain is decomposed into six blocks, and each block is discretised using a structured curvilinear mesh. The surface meshes on the interfaces between adjacent blocks are identical. The numbers of grid points on the block surfaces is $(N_x, N_y, N_z) = (129,129,129)$ and in the wall-normal direction is $N_w = 193$ . The computational domain and grid set-up are qualitatively similar for the flow over the prolate spheroid. The spheroid major-to-minor-axis ratio is $a\,:\,b = 6\,:\,1$ , and the free-stream velocity is $\boldsymbol{U} = (\cos (\alpha ),\sin (\alpha ),0 )$ , where $\alpha =20^{\circ }$ is the incidence angle. The flow domain comprises two concentric spheroids with their axes aligned. The inner spheroid has axis lengths $a=6$ and $b=1$ , while the outer spheroid has aspect ratio close to unity and axis lengths equal to $32.88$ and $32.34$ . The flow domain is decomposed into six blocks and discretised on a structured curvilinear mesh. When we study separation and reattachment, the grid sizes are $(N_x, N_y, N_z,N_w) = (481,241,241,385)$ , which clusters points in on the spheroid. When we study the wake region, the grid sizes are $(N_x, N_y, N_z,N_w)=(241,361,361,769)$ which clusters more grid points in the wall-normal direction to ensure that the wake structures are resolved. The grid quality is verified by how well we achieve forward-adjoint duality for vorticity (e.g. 2.14), which is a stringent test because it requires accurate forward and adjoint simulations. As our results will show, duality is always accurately satisfied.

Table 1. Flow set-up and computational parameters. For flow over a sphere, $R_1$ and $R_2$ are the radii of the inner and outer domain boundaries. For flow over a spheroid, $a_1$ and $b_1$ are the axis lengths of the spheroid, and $a_2$ and $b_2$ are those of the outer spheroidal domain boundary. The resolution $\Delta y_{b,\textit{min}}$ is the wall-normal grid size at the solid wall, $\varDelta _w=(\Delta x_w\Delta y_w\Delta z_w)^{1/3}$ is the grid size at a location three diameters (length of minor axis) downstream of the trailing edge of the sphere and spheroid. For the last case of the flow over the prolate spheroid with impulse in the wake, a smaller grid size $\varDelta _w=0.016$ is adopted in the wake.

The instantaneous velocity fields are stored and used in the adjoint simulations. A continuous adjoint formulation is adopted for the discretisation of the adjoint (2.10), using schemes similar to those for the incompressible Navier–Stokes equations. The advection term $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varOmega }$ and adjoint tilting term $ \boldsymbol{\nabla }\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\varOmega }$ in (2.10) are discretised in time using Adam–Bashforth, and the diffusion term is treated using Crank–Nicolson. We use localised impulses as initial conditions for the adjoint simulations. The regions of interest where we place these impulses are listed in table 1, and the detailed information is presented in the corresponding sections.

Where beneficial, we will use either cylindrical or spherical coordinates for the presentation of numerical results in §§ 4 and 5. The coordinate systems are shown in figure 5. The origins of both systems lie at the centre of sphere or spheroid, and the $x$ -axis is along the free-stream velocity for the sphere and the major axis of the spheroid. The forward-adjoint duality (2.14), which we repeat below, can be expressed in these coordinates

(3.1) \begin{align} \omega _k\left (\boldsymbol{x}_{\kern-1pt f},T \right ) & = \int _{\mathcal{D}}\boldsymbol{\varOmega }^k\left (\boldsymbol{x},s\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V \,\,\, \nonumber \\& \quad - \int _s^T \oint _{\partial \mathcal{D}} \nu \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\varOmega }^k\right ) \boldsymbol{\cdot }\boldsymbol{\omega } \,\textrm{d}S \textrm{d}t + \int _s^T \oint _{\partial \mathcal{D}} \nu \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right )\boldsymbol{\cdot }\boldsymbol{\varOmega }^k \,\textrm{d}S \textrm{d}t \nonumber \\ & = \mathcal{I}(s) +\mathcal{B}_D(s) + \mathcal{B}_N(s). \end{align}

Figure 5. Schematics of the $(a)$ cylindrical and $(b)$ spherical coordinate systems adopted in the presentation of the adjoint-vorticity results. $(a)$ The axial direction is along the $x$ -axis, and the azimuthal angle $\varphi$ is formed by the $-y$ axis and the radial vector. The length of the radial vector is $r$ . $(b)$ The polar axis is along the $x$ -axis. The polar angle $\theta$ is formed between $x$ -axis and the radial vector with length $\eta$ . The definition of the azimuthal angle $\varphi$ is the same as in $(a)$ .

We decompose the interior contribution in the cylindrical coordinates

(3.2) \begin{equation} \mathcal{I}^{k} = \int _{\mathcal{D}} {\boldsymbol \varOmega }^{k}\boldsymbol{\cdot }\boldsymbol{\omega } \,\textrm{d}V, \quad \mathcal{I}_x^{k} = \int _{\mathcal{D}} \varOmega _x^{k}\omega _x \,\textrm{d}V, \quad \mathcal{I}_r^{k} = \int _{\mathcal{D}} \varOmega _r^{k}\omega _r \,\textrm{d}V, \quad \mathcal{I}_\varphi ^{k} = \int _{\mathcal{D}} \varOmega _\varphi ^{k}\omega _\varphi \,\textrm{d}V, \end{equation}

and the time histories of these components reflect the transformation of vorticity among different directions. For flow over a sphere, while the interior contribution is best represented in cylindrical coordinates, the wall contribution is naturally expressed in spherical coordinates. Note that both coordinates share the same azimuthal angle. We adopt the notation

(3.3) \begin{align} B^{k}(\theta ,\varphi ,\tau ) &= \nu \big ( \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right ) \boldsymbol{\cdot }{\boldsymbol \varOmega }^{k} - \big ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } {\boldsymbol \varOmega }^{k}\big ) \boldsymbol{\omega } \big ) , \nonumber\\\overline {B^k}^{\varphi }(\theta ,\tau ) &= \int _0^{2\pi } \nu \big ( \left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right ) \boldsymbol{\cdot }{\boldsymbol \varOmega }^{k} - \big ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } {\boldsymbol \varOmega }^{k}\big ) \boldsymbol{\omega } \big )\textrm{d}\varphi , \nonumber\\\overline {B^k}^{\varphi \tau }(\theta ,\tau ) &= \int _0^{\tau } \overline {B^k}^{\varphi }(\theta ,\xi )\textrm{d}\xi . \end{align}

The boundary term $B^{k}(\theta ,\varphi ,\tau )$ represents the instantaneous wall contribution to the target vorticity at $(\boldsymbol{x}_{\kern-1pt f},\tau =0)$ , from a specific wall location $(\theta ,\varphi )$ at backward time $\tau$ . The azimuthally integrated boundary term $\overline {B^k}^{\varphi }(\theta ,\tau )$ is the instantaneous wall contribution to the target vorticity at $(\boldsymbol{x}_{\kern-1pt f},\tau =0)$ from all azimuthal locations with the same polar angle $\theta$ . The accumulated wall contribution over the backward time interval $(0,\tau )$ is denoted as $\overline {B^k}^{\varphi \tau }(\theta ,\tau )$ .

4. Origin of vorticity in laminar flow over a sphere

In this section we examine the origin of vorticity in laminar flow over a sphere, as an illustrative example of the adjoint approach. Two Reynolds numbers $Re\colon = |\boldsymbol{U}^{*}|D^{*}/\nu ^{*}=\{200,300\}$ are considered. At the lower value, we demonstrate the different earlier-in-time contributions to the vorticity at a point within the steady axisymmetric boundary layer and at separation onset. At the higher Reynolds number, we consider a point in the shed vortices within the wake.

4.1. An illustrative example: axisymmetric separated flow

The flow over a sphere at $Re=200$ is axisymmetric, steady and only contains azimuthal vorticity $\omega _\varphi$ . Figure 6 visualises the streamlines, azimuthal vorticity and the pressure coefficient. An attached boundary layer forms on the forebody, where negative azimuthal vorticity is generated at the wall. The boundary layer separates at $\theta _s = 68^{\circ }$ and carries negative azimuthal vorticity into the wake. A toroidal vortex forms behind the sphere, which is visualised by the closed streamlines in figure 6 $(b)$ . Two points of interest are marked on figure 6, which will be examined using the adjoint approach: the first point ‘B’ is located in the adverse pressure gradient part of the boundary layer, at a wall-normal distance equals to the leading-edge boundary-layer thickness. The second point ‘S’ is at the wall, at the onset of separation. For both ‘B’ and ‘S’, the adjoint impulses will be directed in the azimuthal direction, and by (2.14) we represent the azimuthal vorticity as the summation of interior and boundary contributions. We first consider point ‘B’, which is within the fluid domain. As such, we can evaluate the adjoint evolution with either Neumann or Dirichlet boundary conditions, and compare the two interpretations. When point ‘S’ is considered, only the Neumann condition is compatible with the adjoint pulse being introduced at the wall.

Figure 6. Laminar flow over a sphere and locations of points of interest. $(a)$ Colour contours of $\omega _\varphi$ and streamlines. $(b)$ Contours of the pressure coefficient $C_p=(p-p_{\infty })/(({1}/{2})\rho |\boldsymbol{U}|^2)$ and flow streamlines. $(c)$ Zoomed-in view of $(b)$ , showing the locations of the two points of interest () within the boundary layer and () at the onset of separation.

With the definitions of the terms in (3.1), (3.3), we present the first set of results for case ‘B’. As remarked in § 3, the interior terms in the duality relation are evaluated in cylindrical coordinates (3.1) and the wall terms are reported in spherical coordinates (3.3). Recall that both systems share the same azimuthal direction. The duality relation (2.14) is visualised in figure 7 $(a{,}b)$ for Dirichlet and Neumann boundary conditions. Recall that the former condition evaluates the contribution by the wall vorticity itself, while the latter is in terms of the wall-vorticity flux. The sum of the interior and boundary terms, $\mathcal{I} + \mathcal{B}$ , makes up the target vorticity at B and remains constant through the back-in-time tracking. This conservation verifies the numerical accuracy of the adjoint simulation. Going back in time, the interior contribution $\mathcal{I}$ drops to zero within the considered backward-time horizon, and the boundary terms grow to the value of the target vorticity whether Dirichlet or Neumann boundary conditions are adopted in the adjoint simulation. These results thus track the vorticity at B to its origin at the wall, quantitatively. Figure 7 $(\textit{ii},\textit{iii})$ visualises the instantaneous $\overline {B^{\varphi }}^{\varphi }(\theta ,\tau )$ and integrated $\overline {B^{\varphi }}^{\varphi \tau }(\theta ,\tau )$ boundary contributions versus time. The results show which portion of the wall contributes to the vorticity at B, and how much, as a function of backward time. Such detailed, quantitative and accurate view was not possible before and has never been reported. The figure also shows how the relevant region on the wall shifts upstream with backward time, at a rate set by the adjoint evolution, until the boundary terms sample the vorticity (for Dirichlet) or its flux (for Neumann) at the leading edge.

Figure 7. Duality history and boundary terms for case ‘B’ with $(a)$ Dirichlet and $(b)$ Neumann adjoint boundary conditions. $(i)$ The lines represent () $\mathcal{I}+\mathcal{B}$ , () $\mathcal{I}$ , () $\mathcal{B}$ , () $\mathcal{I}_x$ , () $\mathcal{I}_r$ , () $\mathcal{I}_{\varphi }$ . The boundary contributions are $(a.i)$ $\mathcal{B}=\mathcal{B}_D$ and $(b.i)$ $\mathcal{B}=\mathcal{B}_N$ . $(\textit{ii})$ The instantaneous boundary terms $(a.\textit{ii})$ $\overline {B^{\varphi }_{D}}^{\varphi }(\theta ,\tau )$ and $(b.\textit{ii})$ $\overline {B^{\varphi }_{N}}^{\varphi }(\theta ,\tau )$ . $(\textit{iii})$ The time-integrated boundary terms $(a.\textit{iii})$ $\overline {B^{\varphi }_{D}}^{\varphi \tau }(\theta ,\tau )$ and $(b.\textit{iii})$ $\overline {B^{\varphi }_{N}}^{\varphi \tau }(\theta ,\tau )$ . The polar pressure coefficient profile from DNS is visualised twice identically in $(c,d)$ . The symbols are reference data from Johnson & Patel (1999). The angles $\theta _m$ and $\theta _s$ correspond to the pressure minimum and separation.

The differences between Dirichlet and Neumann boundary conditions are also evident. The Dirichlet boundary terms $\overline {B^{\varphi }_{D}}^{\varphi }$ and $\overline {B^{\varphi }_{D}}^{\varphi \tau }$ reach their final value by $\tau =1.5$ , while the Neumann counterparts require a longer backward time $\tau =3$ . In the former case, once the adjoint field reaches the wall, it is annihilated due to the condition $\boldsymbol{\varOmega }^{\varphi }|_{{wall}}=0$ , and the sampled value of vorticity $(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varOmega }^{\varphi })\boldsymbol{\cdot }\boldsymbol{\omega }|_{{wall}}$ is considered the origin of the wall contribution. In the language of the stochastic representation, once particles reach the wall, they are stopped and their contribution is frozen. In reality, the sampled wall vorticity has its own history, but it is not accessible using the Dirichlet condition. In the Neumann case, $\boldsymbol{\varOmega }^{\varphi }|_{{wall}}\neq 0$ , and we sample the wall-vorticity flux $(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega }^{\varphi })\boldsymbol{\cdot }\boldsymbol{\varOmega }|_{{wall}}$ . When the adjoint reaches the wall, it is effectively retained within the domain due to the boundary condition $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varOmega }^{\varphi }|_{{wall}}=0$ . In the stochastic representation, particles are reflected into the domain, and we continue to sample the vorticity back to its primordial origin. The Neumann condition can thus be regarded as providing a more complete history of the observed vorticity.

Another notable difference is in the signs of the boundary terms. While $\mathcal{B}_{D}$ and $\overline {B^{\varphi }_{D}}^{\varphi }(\theta ,\tau )$ remain negative through all of time, $\mathcal{B}_{N}$ and $\overline {B^{\varphi }_{N}}^{\varphi }(\theta ,\tau )$ are positive for $0\lt \tau \lt 0.7$ . The former terms being always negative is due to wall vorticity remaining negative upstream of separation $\theta _s$ , and therefore $(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varOmega }^{\varphi })\boldsymbol{\cdot }\boldsymbol{\omega }\lt 0$ . In contrast, in the Neumann case, $(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\omega })\boldsymbol{\cdot }\boldsymbol{\varOmega }^{\varphi }$ changes sign because the vorticity flux reverses direction with the pressure gradient changing from favourable to adverse. The wall pressure coefficient is plotted versus the polar angle in figure 7 $(c{,}d)$ , and the location of the pressure minimum $\theta _m$ is marked on the curves. The positive value of $\overline {B^{\varphi }_{N}}^{\varphi }(\theta ,\tau )$ slightly prior to separation is more compatible with classical theory (Lighthill 1963; Morton 1984). It is also important to remark that only the Neumann adjoint boundary condition is applicable if we wish to examine the origin of wall vorticity, i.e. starting from a point on the wall. We will therefore limit our analysis to the Neumann case for the remainder of this work.

To provide a more intuitive understanding of the interior contributions $\mathcal{I}$  (3.1), we visualise the spatial contours of the inner product $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ in figure 8(a). We also plot the forward and adjoint azimuthal vorticities, $\omega _\varphi$ and $\varOmega ^{\varphi }_\varphi$ , in figure 8(b). Initially the inner product $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ is concentrated near the point B. At the larger adjoint times $\tau =\{1,2\}$ , the interior contribution to the target vorticity diffuses and migrates upstream along the boundary. Ultimately, $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ decreases to zero, as the vorticity at the target location has its primordial origin at the wall. On the visualisation plane $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega } = \varOmega _\varphi ^{\varphi }\omega _\varphi$ , thus the individual fields $\varOmega _\varphi$ and $\omega _\varphi$ are informative (figure 8 b). During the back-in-time advection and diffusion of $\varOmega ^{\varphi }_\varphi$ , only regions where $\varOmega _\varphi ^{\varphi }$ and $\omega _\varphi$ are finite contribute to the terminal vorticity. Since $\omega _\varphi$ is appreciable inside the boundary layer and vanishes in the free stream, only the segment of the boundary layer that is sampled by $\varOmega _\varphi ^{\varphi }$ produces large value of $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ . At larger backward time $\tau$ , the adjoint field $\varOmega _\varphi ^{\varphi }$ leaves the near-wall region where $\omega _\varphi$ is finite, and hence the inner product $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ gradually vanishes ( $\tau =3$ ). At this time all of the target vorticity has been traced back to the wall.

Figure 8. Visualisation of the inner product of the forward and adjoint vorticity and their values, on the $x$ - $y$ plane, for case B. The three rows, from bottom to top, correspond to backward times $\tau =\{1,2,3\}$ . $(a)$ Thin solid lines are velocity streamlines. The thick solid line is the separation surface. Colour contours are the interior vorticity contribution $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ . $(b)$ The () positive and () negative values of $\omega _{\varphi }$ vary from $-19$ to $3$ with an increment of $2$ . The zero $\omega _{\varphi }$ contour is highlighted (). Colour contours are the adjoint vorticity $\varOmega _{\varphi }^{\varphi }$ .

We can further interpret the adjoint-vorticity operator $\unicode{x1D76E}$ by evaluating the probability density $\varrho$ and mean deformation matrix $\mathcal{S}$  (2.18), which are examined in figure 9. The evolution of the density $\varrho$ is obtained by solving the Fokker–Planck (2.16) in backward time with Neumann boundary conditions. The matrix $\mathcal{S}$ is calculated by dividing the adjoint vorticity by the density, $\mathcal{S}=\unicode{x1D76E}/\varrho$ . Starting from an initial impulse at B, the probability density function diffuses into a larger region, then migrates upstream due to the back-in-time advection. The probability density in the free stream samples zero vorticity, while only the density within the boundary layer samples finite vorticity. The opaque iso-surface in figure 9 is $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ , coloured by the logarithm of the deformation projected on the directions of the target $ (\boldsymbol{e}_{(\boldsymbol{x}_{\kern-1pt f},\tau =0)} )$ and sampled $ (\boldsymbol{e}_{(\boldsymbol{x},\tau )} )$ vorticity. Note that the structure of $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ is flat and wide in the tangential direction of the wall. Although all points on the iso-surface of $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ contribute equally to the target vorticity, the edge points of the iso-surface experience larger mean deformation since the local density is small.

Figure 9. Visualisation of the probability density function $\varrho =6$ (translucent iso-surfaces) and interior contribution $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }=250$ (opaque iso-surfaces) coloured by $\log (\boldsymbol{e}_{(\boldsymbol{x},\tau )}^{\top }\mathcal{S}\boldsymbol{e}_{(\boldsymbol{x}_{\kern-1pt f},\tau =0)})$ , at $\tau = \{1,2,3\}$ from bottom to top.

We now turn to examining the origin of the zero vorticity at the onset of separation at point ‘S’. This zero vorticity must arise from the exact cancellation of positive and negative vorticity contributions from earlier times, and the herein adopted approach enables a quantitative analysis of how separation develops. The classical description (figure 10) is that a favourable pressure gradient along a surface leads to an influx of negative vorticity that forms the boundary layer; the subsequent adverse pressure gradient is an outflux of negative vorticity (or influx of positive vorticity), which ultimately leads to a vanishing wall vorticity and, therefore, the onset of separation. Implicit in this description is a parabolic view that proceeds from the leading edge downstream to separation onset. Using our herein introduced approach, we will provide a quantitative account of the onset of separation and, for the first time, demonstrate and quantify the role of the pressure gradient downstream of the separation point.

Figure 10. Schematic of Lighthill’s theory for two-dimensional separation. The free-stream flow attaches to the surface at point $A_1$ . The boundary layer is bounded by the dashed line and the solid body. The negative vorticity is denoted by the ‘ $-$ ’ sign, and is oriented along the span into the page. The blue and red arrows represent the flux of negative vorticity into and out of the fluid, respectively. The faint flux arrows beyond the separation surface were not part of the classical description, and their importance is demonstrated herein.

The adjoint vector is initialised as an impulse in the azimuthal direction, along the separation line on the sphere. The total interior and boundary terms in the duality relation (2.14) are shown in figure 11 $(a)$ , and figure 11 $(b{,}c)$ shows the polar distribution of the instantaneous boundary term $\overline {B^{\varphi }_{N}}^{\varphi }(\theta ,\tau )$ and its cumulative contribution in backward time $\overline {B^{\varphi }_{N}}^{\varphi \tau }(\theta ,\tau )$ . The internal $\mathcal{I}$ and Neumann boundary $\mathcal{B}_N$ terms in figure 11 $(a)$ cancel each other, resulting in a conserved zero total vorticity throughout the time horizon. At the end of the considered horizon, the interior contribution nearly vanishes and so does that from the boundary. However, the latter is more interesting since it involves the history of the vorticity flux along the surface. By $\tau =6$ , the azimuthal distribution of the cumulate surface flux, $\overline {B^{\varphi }_N}^{\varphi \tau }$ in figure 11 $(c)$ , shows a cancellation of negative (favourable pressure gradient) and positive (adverse pressure gradient) contributions, the latter originating near the separation point $\theta =\theta _s$ at $\tau \in [0,1]$ . Positive values of $\overline {B^{\varphi }_{N}}^{\varphi }(\theta ,\tau )$ appear on both sides of the separation line. Recall that these fluxes, along with the vanishingly small interior contribution, lead to zero vorticity at separation. The present analysis shows, for the first time, the quantitative contribution to separation onset by the vorticity flux downstream of $\theta =\theta _s$ . Specifically $\mathcal{B}_N(\theta \lt \theta _s)/\mathcal{B}_N(\theta \gt \theta _m)=0.2$ , compared with the adverse pressure gradient contribution upstream of separation $\mathcal{B}_N(\theta _m\gt \theta \gt \theta _s)/\mathcal{B}_N(\theta \gt \theta _m)=0.8$ . These results more clearly delineate the importance of the vorticity flux downstream of separation, which was not highlighted in earlier studies and which should be weaved into our intuition regarding separation. The present analysis is the first to quantify this effect.

Figure 11. $(a)$ Time histories of the terms in the duality relation, including () $\mathcal{I}+\mathcal{B}_N$ , () $\mathcal{I}$ , () $\mathcal{B}_N$ , () $\mathcal{I}_x$ , () $\mathcal{I}_r$ , () $\mathcal{I}_{\varphi }$ . Panels $(b)$ and $(c)$ are instantaneous and time-integrated boundary terms, $\overline {B^{\varphi }_N}^{\varphi }(\theta ,\tau )$ and $\overline {B^{\varphi }_N}^{\varphi \tau }(\theta ,\tau )$ . The angles $\theta _m$ and $\theta _s$ correspond to the pressure minimum and separation.

A closely related phenomenon is material spiking, which has been identified as a precursor to separation (Serra et al. 2018, 2020). Material spiking arises from excessive deceleration of the boundary layer due to downstream resistance, which is typically caused by an adverse pressure gradient. In axisymmetric flow, this deceleration corresponds to large negative values of $\partial u_{\theta }/\partial \theta$ and, by the continuity equation, large wall-normal velocity gradients $\partial u_r/\partial r$ . Using the adjoint-vorticity formulation, we quantitatively relate the two-dimensional separation criterion $\omega _{\varphi }=0$ to the wall-vorticity flux and pressure gradient. Notably, our analysis identifies a strong adverse pressure gradient contribution downstream of the separation location, which reinforces the physical picture of material spiking and provides an additional mechanism contributing to its onset.

Figure 12. $(a$ , $b)$ Coloured contours of the inner product $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ , overlaid with () streamlines, () the separation surface and () the iso-contour $\omega _{\varphi }=0$ . $(c)$ Colour contours of adjoint vorticity and line contours of () positive and () negative forward vorticity $-19 \leqslant \omega _{\varphi } \leqslant 3$ with an increment of $2$ . Also marked is the iso-contour () $\omega _{\varphi }=0$ . From bottom to top, the rows correspond to visualisations at reverse times $\tau =\{0.2, 1.2, 2.2, 3.2\}$ .

The interior contribution $\mathcal{I}$ is only appreciable during $\tau \in [0,4]$ , which we examine in more detail. This term is visualised in figure 12. While the total contribution by this term is negative, at the earliest time $\tau =0.2$ there is a region of positive $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ within the domain. At this early time, the adjoint solution $\varOmega ^{\varphi }_{\varphi }$ has diffused from the impulse at the separation point into a finite region, and samples the forward vorticity $\omega _{\varphi }$ on both sides of the separation surface. Since $\varOmega ^{\varphi }_{\varphi }$ is positive, the sign of the inner product $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ depends on the forward vorticity itself. This inner product is negative upstream of separation and also within a small region downstream of it. However, this term changes sign where $\omega _{\varphi }$ is positive. The bisecting iso-contour $\omega _\varphi = 0$ intersects with the sphere at the separation point. At larger reverse times, the adjoint vorticity $\boldsymbol{\varOmega }^{\varphi }$ upstream of separation is further diffused and advected, sampling the vorticity of the upstream boundary layer and continuing to make a negative contribution. At $\tau =3.2$ , the adjoint solution leaves the sphere surface, and mostly samples the free stream which does not contain any vorticity. This explains why the interior contribution drops to zero a large reverse times. In contrast, the positive interior contribution $\boldsymbol{\varOmega }^{\varphi }\boldsymbol{\cdot }\boldsymbol{\omega }$ remains downstream of the iso-surface of $\omega _\varphi =0$ , and weakens with reverse time.

As remarked above, the positive wall contribution downstream of the separation point (figure 11) and the existence of the positive interior contribution from downstream of the iso-surface $\omega _\varphi =0$ (figure 12) were not discussed in the original literature (e.g. Lighthill 1963). The present analysis exposes these effects. More importantly, the adjoint-vorticity formulation provides a theoretical and quantitative framework to study them.

4.2. Flow over a sphere at $Re=300$

Unsteadiness and vortex shedding start to appear for flow over a sphere at Reynolds number $Re=300$ . In this regime, hairpin vortices are periodically shed behind the sphere in a symmetry plane, as shown in figure 13. The shedding corresponds to a Strouhal number $St := f^{*}D^{*}/ |\boldsymbol{U}^{*} |=0.133$ , or non-dimensional period $1/St = 7.5$ . The symmetry plane has a random azimuthal angle but is fixed. For the presentation of our results, the coordinate system is oriented such that the $x$ - $y$ plane coincides with the plane of symmetry.

Figure 13. $(a)$ Side and top views of the vortical structures for flow over a sphere at $Re=300$ . The iso-surface is $Q = 0.2$ , and is coloured by $\omega _x$ . $(b)$ Visualisation of $\omega _x$ on the plane marked by the dashed line in panel $(a)$ , and () the location of the adjoint initial impulse.

As shown in figure 13, the legs of the shed hairpin vortices contain strong streamwise vorticity. For this instance within the shedding cycle, the detached vortex immediately downstream of the sphere features an extended streamwise segment with appreciable $\omega _x$ . The generation of this streamwise component of vorticity from the earlier azimuthal vorticity on the sphere will be examined. For this purpose, the initial adjoint-vorticity impulse is placed at the marked location in figure 13( $b$ ), and is oriented in the streamwise $x$ direction in order to study the history of $\omega _x$ at this location. Unlike the forward flow, the adjoint solution is transient and does not exhibit a periodic behaviour. The adjoint field advects upstream back in time, samples the vorticity of the forward flow and its wall fluxes and ultimately advects into the free stream where the vorticity is zero. As such, in the limit of long backward time, the full contribution to the observed vorticity is due to the wall fluxes. Our interest is not, however, the long backward-time limit. Instead, we examine the early history and trace back the streamwise vorticity in the separated vortical structures to the azimuthal vorticity within the boundary layer on the sphere. In this manner, we exercise the tilting term that was not part of the previous discussion of axisymmetric separation.

Figure 14. $(a)$ Time history of the terms in the duality relation: () $\mathcal{I}+\mathcal{B}$ , () $\mathcal{I}$ , () $\mathcal{B}$ , () $\mathcal{I}_x$ , () $\mathcal{I}_r$ , () $\mathcal{I}_{\varphi }$ . $(b)$ Vortical structures visualised using translucent iso-surfaces of $Q=0.2$ . The opaque surfaces correspond to the iso-level of the interior contribution, $\boldsymbol{\varOmega }^{x}\boldsymbol{\cdot }\boldsymbol{\omega }=20$ , and are coloured by $\omega _x$ . From bottom to top, the panels correspond to the reverse times $\tau =\{0.6, 2.4, 4.2, 6.0\}$ .

Figure 14(a) shows the history of the terms in the duality relation (2.14). Only the interior term is appreciable, and the wall contribution is small within the time window considered here. Initially, the interior contribution $\mathcal{I}$ consists of only $\mathcal{I}_x$ , while $\mathcal{I}_r=\mathcal{I}_\varphi =0$ . This initial condition is consistent with the set-up of the initial impulse. At larger reverse times, the interior $\mathcal{I}_x$ term gradually reduces to zero. In figure 14 $(b)$ the transparent iso-surfaces show the $Q$ -criterion of the forward field, to provide a sense of the vortical structures present in the flow. The opaque surface is an iso-level of the inner product $\boldsymbol{\varOmega }^{x}\boldsymbol{\cdot }\boldsymbol{\omega }$ and is coloured by $\omega _x$ . This spatial region, which contributes significantly to the target streamwise vorticity, contains progressively less $\omega _x$ going back in time. These visualisations are shown at $\tau =\{0.6,2.4,4.2,6.0\}$ , or at phases $\tau \,St=\{0.08,0.32,0.56,0.80\}$ of the shedding cycle. The adjoint vorticity samples part of the hairpin vortex which shrinks in streamwise extent and tilts upwards at $\tau \sim 2.4$ . The $\omega _x$ contained in the legs of the hairpin vortex is tilted into a large positive $\omega _r$ and negative $\omega _\varphi$ . During this phase of the back-in-time evolution, $\mathcal{I}_r$ and $\mathcal{I}_\varphi$ are both amplifying, although largely opposing one another. By $\tau =4.2$ , the adjoint field expands and encompasses part of the wake and the boundary layer on the sphere surface. The visualised inner product shows two regions where $\varOmega ^{x} \boldsymbol{\cdot }\omega$ is large. At larger time $\tau =6.0$ , the hairpin vortex becomes nearly toroidal and aligns in the azimuthal direction. The duality history (figure 14 a) shows that the target vorticity $\omega _x$ inside the shed hairpin vortex is ultimately traced back to the azimuthal vorticity within the boundary layer. We can thus quantify the interior tilting and stretching of the earlier azimuthal vorticity that led to the streamwise vorticity in the shed vortices. The present analysis thus dispenses with the appeal to qualitative potentialities, for example that the near-wall azimuthal pressure gradient can generate the streamwise vorticity.

5. Origin of vorticity in flow over a prolate spheroid

In this section, we consider the flow over a prolate spheroid, which undergoes three-dimensional open separation. We do not attempt to establish a new theory for open separation, but rather we focus on the Lagrangian dynamics of vorticity at locations along the separation line and within vortical structures that shed from the body.

The Reynolds number based on the minor axis is $Re\equiv Ub/\nu =3000$ , or $18\,000$ based on the major axis, and the incidence angle is $\alpha =20^{\circ }$ (see figure 4). For these parameters, the flow is in the transitional regime. Figure 15 provides a general view of the flow using visualisation of the vortical structures and wall-friction lines. The free stream attaches on the body at the windward centreline. A boundary layer forms from the forebody and extends downstream. Due to the adverse pressure gradient in the azimuthal direction and the surface curvature, the boundary layer separates. The primary separation is identified in the surface stress lines and is marked $S_1$ . The detached shear layers from $S_1$ on the right and left sides curve inward towards the rear point in the cross-flow plane and form a pair of primary counter-rotating vortices which is marked as ‘PV’ in figure 15. Near the centreline on the leeward side, the counter-rotating vortices induce a pair of reverse-flow boundary layers containing negative $v$ velocity, which we refer to as the secondary boundary layers. These secondary boundary layers undergo secondary separation along the line marked $S_2$ . Note that, in the friction patterns depicted in figure 15, separation lines $S_1$ and $S_2$ cannot be traced back to friction zeros, but rather emerge from the sudden convergence of isolated friction lines. The lack of a uniquely identifiable separation line at the beginning of $S_1$ is a characteristic of open separation according to Han & Patel (1979) and Surana et al. (2006). The two counter-rotating primary vortices extrude downstream and break down to turbulence due to shear-layer instability.

Figure 15. $(a)$ Vortical structures visualised using translucent iso-surfaces of $Q=1$ , coloured by $\omega _x$ . $(b)$ Friction lines on the spheroid surface, and contours of the axial velocity $u$ on three $y$ - $z$ planes at $x=\{-2, 0, 2\}$ . $(c)$ Visualisation of in-plane streamlines at five $x$ -locations. Symbols mark the ( $S_1$ , $S_2$ ) separation and ( $R_2$ ) reattachment locations.

Figure 15 $(c)$ shows the $z$ - $y$ in-plane velocity streamlines at five axial positions along the spheroid. The primary and secondary separation points $S_1$ and $S_2$ are marked, and the expansion of the primary vortex is clearly captured. Beneath each primary vortex, two secondary streamwise vortices appear. Between these two secondary vortices, a reattachment line is induced at the wall and is denoted as $R_2$ . Although figure 15 $(b)$ may suggest that $S_2$ and $R_2$ are already present at $x=0$ , closer inspection of the stress lines shows that $S_2$ is not established until farther downstream and the start of $R_2$ is less certain. Considering the velocity fields visualised in figure 15 $(c)$ as a family of two-dimensional dynamical systems that are continuously parameterised by the $x$ location, the emergence of $S_2$ and $R_2$ between $x=0$ and $x=1$ is evidence of a saddle-node bifurcation in this region (Thompson & Stewart 2002).

On each side of the plane of symmetry, the near-wall flow is divided into four regions demarcated by the separation and reattachment lines, $S_1$ , $S_2$ and $R_2$ . Four sample friction lines, $F_1$ $F_4$ , in these regions are visualised in figure 16( $a$ , $b$ ). Starting from the windward centreline, the friction line $F_1$ is representative of the near-wall flow prior to the primary separation $S_1$ . Since the near-wall flow is essentially steady at $x\lt 1.5$ , the friction line approximates the path of near-wall fluid. The line $F_1$ converges to $S_1$ at approximately $x=-0.5$ where the near-wall flow starts to leave the surface (see figure 15 $c$ ). The lines $F_2$ to $F_4$ all start from the leeward side centreline where the free-stream fluid approaches the spheroid due to the streamwise favourable pressure gradient and the downward entrainment induced by the large primary vortices. Using the visualisations of these friction lines, we can relate $R_2$ and $S_2$ on the wall to the vortical structures above the surface. The lines $F_2$ and $F_3$ , essentially coincide before $x=0$ where the secondary boundary layer is attached, and diverge in the region of $R_2$ . The three-dimensional structure corresponding to $R_2$ is the downward entrainment between the counter-rotating near-wall secondary vortices (see figure 15 $c$ , $x\geqslant 1$ ). The entrained fluid is advected towards $S_1$ and $S_2$ , away from $R_2$ , which causes the initially overlapped $F_2$ and $F_3$ to diverge at $R_2$ . The last sample $F_4$ lies along the secondary boundary layer and ultimately converges onto $S_2$ . The flow structure associated with $S_2$ is the ejection of the secondary boundary layer. These separation lines are key features of the wall stress, and are the footprint of the complex vortical structures above the wall. We will examine the formation of these separation and reattachment lines from a vorticity dynamics point of view using the adjoint approach, which unwinds the history of the wall vorticity.

Figure 16. $(a)$ Sample friction lines on the spheroid surface. The grey markers indicate the streamwise locations of the five visualisation planes shown in figure 15 $(c)$ . Colour contours show the pressure coefficient. $(b)$ Friction pattern and sample friction lines mapped onto the $x$ - $\varphi$ coordinates.

5.1. The origin of the null vorticity at three-dimensional separation

On a solid surface, the vorticity component aligned with a friction line vanishes. This fact also holds for three-dimensional separation lines which are regarded as limiting friction lines that attract the neighbouring ones. The null vorticity tangential to the separation line can be traced back in time using the methodology presented in § 2.3 in order to determine the detailed cancellations between stretching/tilting and boundary fluxes in the preceding dynamics.

In the previous sections, the adjoint field was initialised as an impulse $\boldsymbol{\varOmega }^{k} = \boldsymbol{e}_k\delta (\boldsymbol{x}-\boldsymbol{x}_{\kern-1pt f} )$ in order to identify the origin of the terminal vorticity at a point, $\omega _k(\boldsymbol{x}_{\kern-1pt f},T)$ . Here, we make use of another form of the adjoint-vorticity initial condition aimed at examining separation and reattachment

(5.1) \begin{equation} \begin{aligned} \boldsymbol{\varOmega }^{\perp }\left (\boldsymbol{x},T\right ) = \sum _{m=1}^{M}\boldsymbol{\eta }_{(x_m,T)}\delta (\boldsymbol{x}-\boldsymbol{x}_m). \end{aligned} \end{equation}

The above expression corresponds to initial adjoint vector fields that are orthogonal to the vorticity vector, which we indicate by the superscript $\perp$ . These initial fields are therefore aligned along the wall stress vector $\boldsymbol{\tau }_w(x_m,T)$ , with unit vectors $\boldsymbol{\eta }_{(x_m,T)}$ . The index $m=\{1,\ldots ,M\}$ corresponds to $M$ -impulses, and $\boldsymbol{x}_m$ is the spatial location of each impulse at the wall. With this initial condition and the homogeneous Neumann boundary condition, the duality relation (2.12) is re-written as

(5.2) \begin{equation} \begin{aligned} 0=& \int _{\mathcal{D}} \boldsymbol{\varOmega }^{\perp }\left (\boldsymbol{x},s\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V + \int _s^T \oint _{\partial \mathcal{D}} \nu \big ( \boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right ) \big ) \,\textrm{d}S \textrm{d}t. \end{aligned} \end{equation}

The meaning of $\boldsymbol{\varOmega }^{\perp }$ can be readily interpreted from (5.2). The adjoint field $\boldsymbol{\varOmega }^{\perp }$ shows how prior events can superpose to cancel each other and lead to zero vorticity in the $\boldsymbol{\eta }_{(x_m,T)}$ direction, along separation and reattachment lines.

To study the primary separation, a series of vector impulses are placed along $S_1$ using the expression for $\boldsymbol{\varOmega }^{\perp }$ in (5.1). The set-up of the initial impulses, and the time history of contributions to vorticity (5.2), are reported in figure 17(a). Starting from $\tau =0$ , the interior and wall contributions develop negative and positive values, which cancel exactly throughout the back-in-time evolution. The total interior contribution $\mathcal{I}$ appears nearly equal to $\mathcal{I}_x$ , because the other two larger terms, $\mathcal{I}_\varphi$ and $\mathcal{I}_r$ , nearly cancel. This interpretation is, however, potentially misleading as we will demonstrate below. Here we simply caution that the correct interpretation requires visualisation of the fields that make up these integral values. Notwithstanding the details of how $\mathcal{I}=\mathcal{I}_x+\mathcal{I}_\varphi +\mathcal{I}_r$ is established, this term is then perfectly opposed by the wall-vorticity flux to lead to zero vorticity at the target locations and times. Iso-surfaces of the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ are plotted in figure 17 $(b{,}c)$ . The surfaces diffuse around the separation line $S_1$ and travel upstream due to advection. The zoomed-in view in figure 17 $(c)$ reveals a significant difference in the three-dimensional separation relative to the earlier two-dimensional example in § 4.1. A pair of positive and negative values of $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ forms on the two sides of the separation line. The approximate location of the separation surface is identified by the two-dimensional streamlines on a plane orthogonal to $S_1$ in figure 17( $c$ ). The positive isosurface travels, from downstream, across the separation surface and migrates upstream. Such migration across the separation surface was not observed in the two-dimensional case (§ 4.1), and is a feature of the present three-dimensional separation. We will explain this behaviour by examining the details of the terms that make up the interior contribution, $\mathcal{I}$ .

Figure 17. $(a)$ Time histories of the terms in the duality relation: () $\mathcal{I}+\mathcal{B}$ , () $\mathcal{I}$ , () $\mathcal{B}$ , () $\mathcal{I}_x$ , () $\mathcal{I}_r$ , () $\mathcal{I}_{\varphi }$ . $(b)$ Initial impulses at $\tau =0$ (marked by ) and iso-surfaces of the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ at reverse times $\tau = \{0.5, 1.0, 1.5, 2.0\}$ . () Primary separation $S_1$ and () nearby friction lines are visualised on the surface. $(c)$ Zoomed-in view of the region marked by the black square in panel ( $b$ ). The iso-surfaces show the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ . Two-dimensional streamlines of the local forward flow field are shown on a plane orthogonal to $S_1$ , marked by the green lines in panel $(b)$ .

Figure 18. Component-wise interior contributions to the terminal vorticity, at reverse time $\tau =0.5$ . The red and blue iso-surfaces are ( $a$ ) $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ , ( $b$ ) $\varOmega _x^{\perp }\omega _x$ , ( $c$ ) $\varOmega _r^{\perp }\omega _r$ , ( $d$ ) $\varOmega _\varphi ^{\perp }\omega _\varphi$ , with values equal to $\pm 3000$ in the four panels. The separation line $S_1$ , in-plane streamlines normal to $S_1$ , and the separation surface stemming from $S_1$ are also visualised.

Figure 19. Visualisation of the near-wall contributions to the terminal vorticity, at reverse time $\tau =0.5$ . $(a)$ A global view showing the location of the visualisation window relative to the spheroid surface. $(b)$ Zoomed-in view of $(a)$ , displaying the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ on the wall, and the in-plane streamlines (black lines) normal to $S_1$ . Two rectangular regions are marked in $(b)$ , both of which are tangent to the wall at $S_1$ , one at the wall and the second at height $h=0.036$ . Both regions are projected onto the $x$ - $y$ plane and visualised in panels $(c)$ and $(d)$ . ( $c$ – $d$ ) Colour contours show the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ , and yellow arrows are the in-plane adjoint-vorticity vectors $\boldsymbol{\varOmega }^{\perp }$ .

Figure 18 shows the iso-surface of the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ and its components $\varOmega _x^{\perp }\omega _x$ , $\varOmega _r^{\perp }\omega _r$ and $\varOmega _\varphi ^{\perp }\omega _\varphi$ , at $\tau =0.5$ (same as in figure 17 $c$ ). The iso-surface of the azimuthal part $\varOmega _\varphi ^{\perp }\omega _\varphi$ resembles that for the total inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ . Although the other two components, $\varOmega _x^{\perp }\omega _x$ and $\varOmega _r^{\perp }\omega _r$ , contribute to the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ , they do not alter its overall structure. This observation addresses the potential misconception from figure 17, that $\mathcal{I}$ is due to $\mathcal{I}_x$ and that $\mathcal{I}_r$ and $\mathcal{I}_\varphi$ nearly cancel. The present figure 18 shows the importance of $\mathcal{I}_\varphi$ . We therefore direct our attention on the $\varphi$ adjoint equation. We consider a local coordinate system $(x', r', \varphi ')$ as shown in figure 19, where $x'$ is in the direction along the separation line, $r'$ is orthogonal to the wall and $\varphi '$ is orthogonal to the separation line and tangent to the wall. For the $\varphi '$ -adjoint-vorticity equation, the advective terms vanish at the wall due to the no-slip condition on the velocity field. Additionally, since the adjoint field is initialised with only an $x'$ -component (so $\varOmega _{r'}^{\perp }=\varOmega _{\varphi '}^{\perp }=0$ at $t=0$ ), the diffusion term and two of the tilting terms are also zero at that instant, leaving the initial-time equation

(5.3) \begin{equation} \frac {\partial \varOmega ^{\perp }_{\varphi '}}{\partial \tau }= \varOmega ^{\perp }_{x'}\frac {\partial u_{x'}}{\partial \varphi '}. \end{equation}

Equation (5.3) demonstrates that the generation of non-zero $\varOmega ^{\perp }_{\varphi '}$ is due to tilting of the adjoint $\varOmega _{x'}^{\perp }$ . Initially the adjoint field $\varOmega ^{\perp }_{x'}$ is positive, and diffuses from a Dirac delta located on the separation line to both sides of the line. The sign of $\varOmega ^{\perp }_{\varphi '}$ at later times depends on the sign of $\varOmega ^{\perp }_{x'}( {\partial u_{x'}}/{\partial \varphi '})$ . Since the flow satisfies the no-slip boundary condition, the velocity $u_{x'}$ at separation line is a local minimum in the $\varphi '$ coordinate compared with points on the tangent plane. We immediately see that, for locations where $\varphi '\gt 0$ , $ ({\partial u_{x'}}/{\partial \varphi '}) \gt 0$ and, as a result, $\varOmega ^{\perp }_{\varphi '}\gt 0$ by (5.3). The same analysis carries over to the other half of the tangent plane with $\varphi '\lt 0$ . We ultimately observe a pattern of the $\boldsymbol{\varOmega }^{\perp }$ vector that is predominantly parallel to, but slightly diverging from, the separation line (see figure 19 $c$ ). Since the vorticity at the vicinity of $S_1$ points into the $\varphi '$ direction, the diverging pattern of $\boldsymbol{\varOmega }^{\perp }$ results in a pair of $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ regions with opposite signs. The adjoint field $\boldsymbol{\varOmega }^{\perp }$ then diffuses across the separation surface. The visualisation in figure 19 $(d)$ shows an off-wall plane, upstream of the separation surface based on the azimuthal component of velocity. The positive value of $\varOmega _{x'}^{\perp }$ in this plane arises from diffusion of positive $\varOmega _{x'}^{\perp }$ across the separation surface from the downstream side, while negative values are directly generated from the tilting effect in the upstream region. This pair of positive and negative $\varOmega _{x'}^{\perp }$ samples the local $\omega _{\varphi }$ with opposite contributions for the target null vorticity.

In summary, the migration of the positive contribution to vorticity across the separation surface, from the downstream to the upstream side, is due to a combination of tilting and viscous diffusion. In contrast, the tilting effect was not active in two-dimensional separation (§ 4.1), and this migration across the separation surface did not exist. The secondary separation $S_2$ has previously been explained from an Eulerian perspective using the theory of vortex-induced separations (Du & Zaki 2024). For the Lagrangian interpretation, We traced the null vorticity along $S_2$ back in time. Similarly to $S_1$ , a pair of positive and negative interior contributions emerges on the two sides of $S_2$ . The negative interior contribution on the downstream side of $S_2$ crosses the separation surface. The results and conclusions from the back-in-time analysis of $S_2$ are qualitatively consistent with those for $S_1$ , and are therefore not included.

Figure 20. $(a)$ Locations () of the initial adjoint-vorticity impulses along $R_2$ . The surface friction lines are shown in grey, and the in-plane streamlines orthogonal to the reattachment line $R_2$ are shown in brown. $(b)$ Visualisation of the marked black rectangle in panel ( $a$ ), showing iso-surfaces of $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }=\{-5000,5000\}$ at reverse times $\tau =\{1.0, 1.5, 2.0\}$ . $(c)$ Time histories of the duality relation. () $\mathcal{I}+\mathcal{B}_N$ ; () $\mathcal{I}$ ; () $\mathcal{B}_N$ ; () $\mathcal{I}_x$ ; () $\mathcal{I}_r$ ; () $\mathcal{I}_{\varphi }$ .

Figure 21. (a) Contours of axial vorticity, $\boldsymbol{\omega _x}$ . The vectors are tangent to the in-plane vorticity. (b) Contours of the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ . $(c,d,e)$ Contours of the component-wise products of the forward and adjoint vorticity. The marker $R_2$ shows the location of secondary separation in this plane.

The set-up of initial adjoint impulses along the separation line is replicated along the reattachment line $R_2$ . This reattachment is established between two near-wall small counter-rotating vortices, and is induced by the downward mass flux towards the wall (see figure 15 $c$ ). The initial impulses and the friction lines are visualised in figure 20 $(a)$ along with the streamlines in a cross-sectional plane. The initial condition consists of impulses aligned tangential to the reattachment line, since we seek to quantify how the null vorticity along this line is established. Iso-surfaces of the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ are visualised in figure 20 $(b)$ at different backward times. Similar to the analysis of separation, a pair of positive and negative contributions appears on the two sides of the reattachment line at $\tau =1.0$ . The back-in-time structure of $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ depends on the local flow. The impulses of $\boldsymbol{\varOmega }$ and therefore the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ near the downstream end of the reattachment line are confined by the two near-wall small counter-rotating vortices beneath the primary vortex. In contrast, the impulses near the upstream end of the reattachment line are advected upstream, leave the streamwise vortices above $R_2$ , and enter the attached boundary layer. As a result, an elongated region contributes to the terminal zero-vorticity vector at reattachment. The time history of the interior and boundary terms are shown in figure 20 $(c)$ . The zero target vorticity is mostly due to large cancellation between the azimuthal and the radial contributions. These large terms are dissected in figure 21, which shows a $y$ - $z$ plane, at $x=0.28$ and $\tau =1.5$ . Contours of the forward axial vorticity are plotted in panel ( $a$ ), the inner product $\boldsymbol{\varOmega }^{\perp }\boldsymbol{\cdot }\boldsymbol{\omega }$ is plotted in panel ( $b$ ), and the three entries of the inner product are shown in panels ( $c$ – $e$ ). The axial contribution $\varOmega _x^{\perp }\omega _x$ is finite at locations where the line contours of $\varOmega _x^{\perp }$ and $\omega _x$ overlap in figure 21 $(c)$ . However, the magnitude of $\varOmega _x^{\perp }\omega _x$ is much smaller than the other two terms. The $y$ - $z$ in-plane components of $\boldsymbol{\varOmega }^{\perp }$ are nearly orthogonal to the wall, pointing into the fluid, while the vorticity vectors are towards the wall (figure 21 $d{,}e$ ). The outcome is a negative radial contribution $\varOmega _r^{\perp }\omega _r$ . The in-plane adjoint $\varOmega ^{\perp }_{y}\boldsymbol{e}_y+\varOmega ^{\perp }_{z}\boldsymbol{e}_z$ has a small positive azimuthal component, which samples the positive $\omega _\varphi$ and results in a positive $\varOmega _\varphi \omega _\varphi$ . The oppositely signed $\varOmega _r^{\perp }\omega _r$ and $\varOmega ^{\perp }_\varphi \omega _\varphi$ partially cancel one another, and produce the total interior contribution shown in figure 21 $(b)$ . From the perspective of the forward vorticity evolution, the negative radial component $\omega _r$ is tilted into negative $\omega _x$ by the radial gradient of streamwise velocity. The azimuthal vorticity $\omega _\varphi$ , on the other hand, is tilted into positive $\omega _x$ by the azimuthal gradient of streamwise velocity. The combination of these effects lead to the null vorticity on the reattachment line. We quantitatively justify this argument by simplifying the adjoint-vorticity equation, and demonstrate that the adjoint tilting term is indeed the dominant effect in the early stages of the adjoint evolution.

In analogy to the physical description of vorticity dynamics, here, we introduce an intuitive description of the evolution of the adjoint. To aid the interpretation, we recall that in the inviscid limit the vorticity equation ${\textrm{d}\boldsymbol{\omega }}/{\textrm{d}t}=\boldsymbol{\omega }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}$ is the same as that of an infinitesimal material line element, $\delta \boldsymbol{l}$ . In that same limit, $\boldsymbol{\varOmega }$ satisfies ${\textrm{d}\boldsymbol{\varOmega }}/{\textrm{d}t}=- \boldsymbol{\nabla }\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\varOmega }$ which is the same as the equation for the evolution of an area element $\delta \boldsymbol{S}$ (Xiang et al. 2025) that also undergoes stretching and tilting. These effects are examined with aid of a local coordinate system, $(x',r',\varphi ')$ , where $x'$ is along the reattachment line $R_2$ , $r'$ is orthogonal to the wall and $\varphi '$ is orthogonal to $R_2$ and tangent to the wall. Note that the $r'$ coordinate is nearly aligned with the radial direction $r$ for most of the reattachment line. Following the same simplifications used to obtain (5.3), the initial-time adjoint-vorticity equation reduces to

(5.4) \begin{equation} \frac {\partial \varOmega ^{\perp }_{r'}}{\partial \tau }= \frac {\partial u_{x'}}{\partial r'}\varOmega ^{\perp }_{x'}. \end{equation}

Figure 22. Schematic of vorticity and adjoint-vorticity tilting in linear shear flow. $(a.i$ – $\textit{ii})$ Forward tilting of vorticity $\boldsymbol{\omega }(s)$ into $\boldsymbol{\omega }(t)$ . $(b.i$ – $\textit{ii})$ Forward and backward tilting of the adjoint vorticity $\boldsymbol{\varOmega }^{\perp }$ .

The velocity gradient $ {\partial u_{x'}}/{\partial r'}$ is the dominant shear in the boundary layer, thus a significant temporal growth rate of $\varOmega ^{\perp }_{r'}$ is provided by the adjoint tilting term, which is shown schematically in figure 22. The top row shows the well-established tilting of the forward vorticity (or infinitesimal material line element). Proceeding from figure 22( $a.i$ – $\textit{ii}$ ), $\omega _r'$ is tilted by ${\partial u_{x'}}/{\partial r'}$ and the source term $\omega _{r'}( {\partial u_{x'}}/{\partial r'})$ , thus generating $\omega _{x'}$ in forward time. The second row is first explained in forward time, proceeding from figure 22(b.i)–(b.ii) (replace $\tau =-t$ in (5.4)). An area vector $\boldsymbol{\varOmega }^{\perp }$ with components in both the $(x',r')$ directions is tilted by $-( {\partial u_{x'}}/{\partial r'})$ to generate $\boldsymbol{\varOmega }^{\perp }$ in the $x'$ direction only. In backward time, which is our interest in adjoint analysis, the motion is reversed going from figure 22 $b.\textit{ii}$ – $b.i$ . Stated explicitly, the back-in-time evolution of $\boldsymbol{\varOmega }^{\perp }$ involves $\varOmega ^{\perp }_{x'}$ being tilted by $ {\partial u_{x'}}/{\partial r'}$ to generate $\varOmega _{r'}$ . When viscosity is present, only the back-in-time description of the adjoint-tilting term is well posed, and this conceptual model explains the initial generation of $\varOmega _{r}$ at the reattachment line $R_2$ .

In summary, the zero vorticity along the reattachment line is formed by large cancellation of radial and azimuthal contributions from earlier vorticity within the flow. These radial and azimuthal components are tilted into positive and negative vorticity aligned with the reattachment friction line, and cancel each other at the wall.

5.2. History of vorticity in detached vortical structures

We now turn to the vorticity within the field, and examine two defining vortical features of this flow. Specifically, we consider the primary counter-rotating vortices that are formed on the leeward side, and the more complex necklace vortices that are observed in the wake. The primary vortices are important since they influence the lift force and three-dimensional separation on the body. As for the necklace vortices, they form through a complex vortical interaction in the near wake, and their back-in-time tracing will serve as a final demonstration of the adjoint-vorticity approach.

Figure 23. $(a)$ Locations () of the initial adjoint-vorticity impulses. $(b.i$ – $\textit{iii})$ Interior contributions to the terminal vorticity, visualised using opaque iso-surfaces $\boldsymbol{\varOmega }^{\parallel }\boldsymbol{\cdot }\omega = 1000$ , coloured by $\varOmega ^{\parallel }_x\omega _x$ , at reverse times $\tau =\{1, 2, 3\}$ . $(c.i$ – $\textit{iii})$ Zoomed in view of the marked region in panel $(b.\textit{iii})$ , showing the iso-surface $\boldsymbol{\varOmega }^{\parallel }\boldsymbol{\cdot }\omega = 1000$ , coloured by $(c.i)$ $\varOmega ^{\parallel }_x\omega _x$ , $(c.\textit{ii})$ $\varOmega ^{\parallel }_r\omega _r$ and $(c.\textit{iii})$ $\varOmega ^{\parallel }_\varphi \omega _\varphi$ . In all panels, one of the primary vortices is visualised by translucent iso-surfaces of $Q=2$ . Streamlines in $y$ - $z$ planes at $x=\{-2, -1, 0, 1, 2\}$ are displayed for $z\lt 0$ .

Figure 24. Time histories of the terms in the duality relation, for the origin of the vorticity in the primary vortex. () $\mathcal{I}+\mathcal{B}_N$ ; () $\mathcal{I}$ ; () $\mathcal{B}_N$ ; () $\mathcal{I}_x$ ; () $\mathcal{I}_r$ ; () $\mathcal{I}_{\varphi }$ .

The primary counter-rotating vortices are the most prominent structures in the flow over the prolate spheroid. These vortical structures are observed at the moderate Reynolds number considered in this study (figure 15) and in the super-critical conditions studied by Fu et al. (1994) and Chesnakas & Simpson (1996) where laminar–turbulent transition occurs upstream of the primary separation. The lift and lateral forces on the body are closely related to the axial vorticity $\omega _x$ within the primary vortices (Du & Zaki 2024). Using the adjoint-vorticity framework, we examine the origin of the forward vorticity at the core of one of the primary vortices. A series of adjoint-vorticity impulses are initialised along the centreline of the primary vortex, in the region $0\leqslant x\leqslant 2.5$ as shown in figure 23 $(a)$ . These initial impulses are defined as

(5.5) \begin{equation} \begin{aligned} \boldsymbol{\varOmega }^{\parallel }\left (\boldsymbol{x},T\right ) = \sum _{m=1}^{M}\boldsymbol{e}_{(x_m,T)}\delta (\boldsymbol{x}-\boldsymbol{x}_m). \end{aligned} \end{equation}

The coordinates $\boldsymbol{x}_m$ are the spatial locations of the impulses, which are indexed $m = 1, 2,\ldots ,M$ . The unit vector $\boldsymbol{e}_{(x_m,T)}$ is parallel to the vorticity vector $\boldsymbol{\omega }(x_m,T)$ . With the homogeneous Neumann boundary condition, the corresponding duality relation (2.12) is re-written as

(5.6) \begin{equation} \begin{aligned} \sum _{m=1}^{M}\vert \boldsymbol{\omega }(\boldsymbol{x_m},T)\vert =& \int _{\mathcal{D}}\boldsymbol{\varOmega }^{\parallel }\left (\boldsymbol{x},s\right ) \boldsymbol{\cdot }\boldsymbol{\omega }\left (\boldsymbol{x},s\right ) \textrm{d}V +\int _s^T \oint _{\partial \mathcal{D}} \nu \left ( \boldsymbol{\varOmega }^{\parallel } \boldsymbol{\cdot }\left ( \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\omega } \right )\right ) \,\textrm{d}S \textrm{d}t. \end{aligned} \end{equation}

The meaning of $\boldsymbol{\varOmega }^{\parallel }$ can be readily gleaned from (5.6). The adjoint field $\boldsymbol{\varOmega }^{\parallel }$ samples the earlier vorticity and wall contributions that make up the terminal vorticity at the locations of the impulses, and is useful to find the origin of vorticity in inhomogeneous vortical structures. The isosurfaces of the inner product $\boldsymbol{\varOmega }^{\parallel } \boldsymbol{\cdot }\boldsymbol{\omega }$ are shown in figure 23(b.i–iii) at $\tau =\{1,2,3\}$ . At $\tau =1$ the iso-surface of $\boldsymbol{\varOmega }^{\parallel } \boldsymbol{\cdot }\boldsymbol{\omega }$ captures large positive values of $\varOmega _x^{\parallel }\omega _x$ , because the terminal vorticity is primarily due to earlier $\boldsymbol{\omega }_x$ on this time scale. In fact, the vorticity at the impulse locations only contains an axial contribution as seen from the terms in the duality relation reported in figure 24. With backward time, the interior contributions migrate along the primary vortex towards the attached boundary layer. When the adjoint, and therefore the isosurface $\boldsymbol{\varOmega }^{\parallel } \boldsymbol{\cdot }\boldsymbol{\omega }$ , enters the attached boundary layer, the magnitude of $\varOmega _x^{\parallel }\omega _x$ drops to nearly zero, as shown in figure 23(b.ii–iii). Correspondingly the duality history in figure 24 shows significant drop of $\mathcal{I}_x$ near the end of the adjoint evolution, and the original generation of the primary vortex by $\mathcal{I}_\varphi$ near the leading edge becomes visible. Longer tracking would relate the terminal vorticity to the wall-vorticity flux as shown in § 4.1.

The generation of the streamwise vorticity in the primary vortex was also explored by Lighthill (1963), who anticipated that this vorticity arises on a body surface due to the excessive curvature of wall-friction lines relative to external potential-flow streamlines. Lighthill (1963) attributed this excessive curvature to an imbalance between the reduced centrifugal force of the near-wall flow and the external pressure gradient acting on the fluid elements. In contrast, our analysis uncovers a different mechanism: the streamwise vorticity observed within the primary vortices originates from the tilting of vorticity components from the $y$ - $z$ plane into the $x$ -direction, rather than from pre-existing $x$ vorticity in the boundary layer. This distinction establishes a new perspective on the dynamics of vorticity generation in such flows.

Figure 25. $(a)$ Vortical structures around and in the wake of the prolate spheroid. Translucent surfaces are $Q=2$ , coloured by $\omega _x$ . Two visualisation planes, $V_1$ and $V_2$ , are marked. $(b)$ Contours of $\omega _x$ on $V_1$ . $(c)$ Contours of $\omega _z$ on $V_2$ which is aligned with the $x$ - $y$ plane. $(d)$ Bottom view of the near wake, showing the iso-surface $Q=2$ coloured by $\omega _z$ . The grey arrow points to the target necklace vortex whose history will be analysed. Panels $(c)$ and $(d)$ share the same $x$ range.

The transitional nature and three-dimensionality of the flow over the spheroid poses difficulty for a clear analysis of the near wake. Two key vortical structures can be identified near the trailing edge of the spheroid in the visualisations in figure 25. The shown plane $V_1$ captures the first element, which is the counter-rotating primary vortices discussed above. The second element is the vortex shedding caused by the top and bottom trailing-edge detached boundary layers. This shedding features large negative and positive spanwise vorticity $\omega _z$ , as shown in figure 25(c,d). The strong, localised spanwise vorticity in the $V_2$ plane (figure 25 $c$ ) corresponds to the necklace-shaped vortices wrapped around the primary vortex pair (figure 25 d). This type of secondary structure that contains vorticity orthogonal to the primary vortices is frequently encountered in vortex-dominated flows, including flow over the prolate spheroid (Jiang et al. 2016) and vortex–boundary-layer interaction (Mole, Skillen & Revell 2022). The origin of the necklace can be anticipated to be the boundary layer which contains spanwise vorticity that is stretched by the primary vortices. We verify and quantify this connection using the adjoint-vorticity approach. The target vortex is marked in figure 25 $(d)$ . We adopt the adjoint initial condition

(5.7) \begin{equation} \boldsymbol{\varOmega }^{v}(\boldsymbol{x},T)=\begin{cases} {\boldsymbol{\omega }}/{\vert \boldsymbol{\omega }\vert }, \qquad & \boldsymbol{x}\in \mathcal{R} ,\\ 0, \qquad & \text{otherwise}, \end{cases} \end{equation}

where $\mathcal{R}$ is the region occupied by the target vortex. At each point within $\mathcal{R}$ , a unit vector is defined pointing in the direction of the local forward vorticity. With this adjoint initial condition at initial time

(5.8) \begin{equation} \int \boldsymbol{\varOmega }^{v}\boldsymbol{\cdot }\boldsymbol{\omega }(\boldsymbol{x},T)\textrm{d}V=\int _{\mathcal{R}}\vert \boldsymbol{\omega }(\boldsymbol{x},T)\vert \,\textrm{d}V. \end{equation}

The back-in-time evolution of the adjoint vorticity from this initial condition identifies the origin of the total amount of $|\boldsymbol{\omega }|$ inside the target vortex tube.

Figure 26. $(a)$ Time histories of the terms in the duality relation, for the origin of the vorticity in the necklace vortex. () $\mathcal{I}+\mathcal{B}_N$ ; () $\mathcal{I}$ ; () $\mathcal{B}_N$ ; () $\mathcal{I}_x$ ; () $\mathcal{I}_y$ ; () $\mathcal{I}_{z}$ . $(b$ – $d)$ Top, side and bottom views of (translucent iso-surface) instantaneous vortical structures in the wake, visualised using $Q=2$ . Opaque iso-surfaces show the interior contribution $\boldsymbol{\varOmega }^v\boldsymbol{\cdot }\boldsymbol{\omega }=0.05$ to the terminal vorticity of the necklace vortex at (right to left) reverse times $\tau =\{0, 1.2, 2.4, 3.6, 4.8\}$ , coloured by $\omega _z$ .

The histories of the terms in the duality relation and iso-surfaces of the inner product $\boldsymbol{\varOmega }^{v}\boldsymbol{\cdot }\boldsymbol{\omega }$ at different backward times are shown in figure 26. The three interior contributions $\{\mathcal{I}_x, \mathcal{I}_y, \mathcal{I}_z\}$ oscillate, while the total amount of vorticity is conserved with good accuracy. The inner product is initially distributed along a necklace vortex tube at $x\approx 7$ , and is advected upstream and towards the bottom of the primary vortices in reverse time. By $\tau =3.6$ the iso-surface of the inner product $\boldsymbol{\varOmega }^{v}\boldsymbol{\cdot }\boldsymbol{\omega }$ is split into two. One part returns to the rear end of the spheroid surface, which corresponds to the vorticity contribution from the shedding process depicted in figure 25 $(c)$ , and the other part is entrained into the primary vortices. The adjoint-vorticity analysis thus demonstrates, objectively, that the origin of the necklace vortices is not solely due to the trailing-edge boundary layer $(71\,\%)$ , and that the terminal vorticity in the necklace structure has a similarly relevant contribution from the large-scale primary vortices $(29\,\%)$ . This point is affirmed in figure 27, where we visualise isosurfaces of $\boldsymbol{\varOmega }^{v}\boldsymbol{\cdot }\boldsymbol{\omega }$ in the $y$ - $z$ plane along with contours of the forward vorticity field. Here, we describe the forward evolution as discovered using the adjoint approach. A clear two-stage process can be identified for the formation of the necklace vortices. First, the vorticity from the primary vortices is combined with the detached shear layer at the spheroid trailing edge (figure 27 $a$ – $b$ ). Subsequently the shed shear layer is stretched in the spanwise direction and wraps around the primary vortices (figure 27 $b$ – $c$ ). This stretching effect is itself induced by the primary vortices, thus strengthening the necklace vortex.

Figure 27. Visualisation of the interior contribution to the vorticity in the necklace vortex. Contours show $\omega _x$ in $y$ - $z$ planes. Iso-surfaces correspond to $\boldsymbol{\varOmega }^v\boldsymbol{\cdot }\boldsymbol{\omega }=0.05$ , and are coloured by ${\varOmega _z}^v{\omega _z}$ . ( $a$ – $c$ ) The three planes correspond to three different $x$ -locations and reverse times, starting at $(c)$ $x=6.7$ at $\tau =0$ and moving upstream to $(c)$ $x=3.0$ at $\tau =4.8$ . The three planes align with the dashed lines in figure 26( $c$ ).

6. Conclusions

In this work we examine back-in-time origin of vorticity in separated flows over immersed bodies. The approach relies on the recently introduced adjoint-vorticity equations, which are equivalent to the stochastic Lagrangian interpretation of vorticity (Xiang et al. 2025). By evolving the adjoint-vorticity equations back in time, we compute the volume density of mean deformation of the earlier vorticity that contributes to the target value, through stretching and tilting. The formulation also identifies the surface contribution, either due to the forward vorticity when homogeneous Dirichlet conditions are enforced on the adjoint variable or due to the wall-vorticity flux when homogeneous Neumann conditions are adopted on the adjoint. Only the latter boundary condition can be adopted when tracing back the origin of wall vorticity, which is of particular interest when attempting to discover the precursors of separation.

Three examples are considered for different purposes, including laminar flows over a sphere at $Re=200$ and $300$ based on the diameter, and transitional flow over a prolate spheroid at minor axis $Re=3000$ . Forward direct numerical simulations are performed for these flows, and the adjoint-vorticity equations are solved back in time to study the origin of the vorticity at points of interest.

For the flow over a sphere with $Re=200$ , we find that the zero vorticity at two-dimensional separation originates from the cancellation of positive and negative vorticity fluxes at the wall, which is in apparent, or partial agreement with Lighthill’s theory (Lighthill 1963). The analysis demonstrated that a wall contribution from downstream of separation accounts for $20\,\%$ of the positive vorticity flux into the domain, which is appreciable and essential for achieving a net zero vorticity at separation. This downstream contribution completes Lighthill’s picture. In flow over a sphere at $Re=300$ , the origin of the streamwise vorticity in the shed hairpin vortices is quantitatively traced back to the azimuthal vorticity inside the boundary layer. The tilting process of the original azimuthal vorticity into the streamwise component is captured exactly by the terms that make up the duality relation.

In the more complex example of the flow over a prolate spheroid, the null vorticity along the three-dimensional separation line is realised by the cancellation of local wall-vorticity flux and interior vorticity from far upstream. A distinctive feature of the three-dimensional separation is that the contribution from the interior vorticity from downstream of the separation surface has its origin upstream of separation. In other words, the adjoint field from downstream of separation crosses the separation surface. This feature is caused by strong vorticity tilting that is not active in the two-dimensional separation. In addition, we examined the origin of the vorticity in the detached vortical structures that are shed from the spheroid body. First, the streamwise vorticity in the primary vortices on the leeward side originates from the tilting of azimuthal vorticity in the attached boundary layer near the leading edge. This newly established connection addresses an open question in the literature, where it was hypothesised that the primary vortices may be due to the streamwise vorticity flux at the wall. Second, in the downstream wake, the necklace vortices are formed by the shedding of spanwise vorticity from the trailing edge, and by a stretching effect induced by the primary vortices. The accurate back-in-time analysis of the vorticity within the transitional wake demonstrates the capability of the adjoint-vorticity approach in the analysis of complex vortex-dominated flows.