2.1 Flow governing equations

In a 3-D fluid flow, the Eulerian velocity $\boldsymbol {u}$ is governed by the NS equations, including the continuum equation

(2.1)\begin{equation} \frac{\partial{\rho}}{\partial{t}} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{u}) = 0, \end{equation}

the momentum equation

(2.2)\begin{equation} \frac{\partial{\boldsymbol{u}}}{\partial{t}} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} = \boldsymbol{\nabla} \varPi + \boldsymbol F, \end{equation}

the energy equation and the equation of state. Here, $\varPi = -{p}/{\rho }+\varXi$ denotes a generalized potential with the pressure $p$, the density $\rho$ and the potential $\varXi$ of conservative body forces, and

(2.3)\begin{equation} \boldsymbol{F} = \frac{1}{\rho}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol\tau) + p\boldsymbol{\nabla}\left(\frac{1}{\rho}\right) + \boldsymbol{f} \end{equation}

denotes a generalized non-ideal force term, where $\boldsymbol \tau$ is the viscous stress tensor and $\boldsymbol f$ the non-conservative external body force per unit mass.

An ideal flow, which is governed by the Euler equations, has $\boldsymbol F=0$. In a non-ideal flow, the three terms on the right-hand side of (2.3) denote the viscous term in most real flows, the baroclinic term, e.g. in compressible flows and combustion, and the non-conservative body force term, e.g. the Lorentz force in magnetohydrodynamic (MHD) flows, respectively.

Taking the curl of (2.2) yields the transport equation

(2.4) \begin{equation} \frac{\partial{\boldsymbol{\omega}}}{\partial{t}}-\boldsymbol{\nabla}\times(\boldsymbol{u}\times\boldsymbol{\omega}) = \boldsymbol{\nabla}\times\boldsymbol{F}, \end{equation}

for the vorticity $\boldsymbol {\omega }\equiv \boldsymbol {\nabla }\times \boldsymbol {u}$. The Helmholtz theorem illustrates a frozen-in nature of $\boldsymbol {\omega }$ in ideal flows with $\boldsymbol {F}=0$, i.e. we can track a vortex surface convected by $\boldsymbol u$. However, the Helmholtz theorem breaks down in non-ideal flows with $\boldsymbol \nabla \times \boldsymbol F\neq 0$ (Reference TruesdellTruesdell, 1954). As sketched in figure 1(a), if we track a surface which is a vortex surface at the initial time using $\boldsymbol u$, the surface is no longer a vortex surface at later times. This is a major obstacle to developing a Lagrangian formulation of the vortex dynamics in real flows.

Figure 1.

Sketch of the surface evolution convected by the (a) physical velocity and (b) virtual circulation-preserving velocity from an initial time to a later time in a viscous flow. The initial vortex tube consists of vortex lines. The vortex-surface tracking breaks down in (a) whereas it holds in (b).

2.2 The VSF equation and exact solutions

The VSF $\phi _v$ is a 3-D globally smooth scalar field. In the flow visualization using the VSF or other vortex criteria, the three-component vorticity is converted into a scalar. In this way, an isosurface of such a scalar with an appropriate isocontouring value can represent typical vortical structures, which suggests a dimension reduction of flow information.

In the definition of the VSF, the vorticity is tangent at every point (except for vorticity nulls) on an isosurface of $\phi _v$. Some typical VSF isosurfaces, e.g. vortex surfaces, are shown in figure 2. Thus the VSF satisfies the constraint

(2.5)\begin{equation} \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi_v = 0. \end{equation}

Note that it is non-trivial to construct or evolve the VSF under this constraint both theoretically and numerically.

Figure 2.

Isosurfaces of the exact VSF: (a) vortex column in the integrable ABC flow, (b) vortex ring and (c) knotted vortex tube. Some vortex lines are integrated from points on the isosurfaces.

From the theorem of Reference ErtelErtel (1942), the evolution equation of the VSF can be written in a Lagrangian conserved form in an ideal flow (Reference Hao, Xiong and YangHao et al., 2019). As sketched in figure 1(b), the evolution of vortex surfaces in viscous flow is conceptually equivalent to tracking vortex surfaces using the virtual circulation-preserving velocity.

For given $\boldsymbol {\omega }$, the VSF solution can be obtained by solving (2.5), but the exact VSF solution only exists in special cases. From the Frobenius theorem, a simple flow field with the vanishing helicity density $h\equiv \boldsymbol u\boldsymbol {\cdot } \boldsymbol {\omega }$ can have non-unique exact VSF solutions (Reference He and YangHe & Yang, 2016). By solving characteristic equations of $\boldsymbol {\omega }$, several globally smooth exact VSF solutions were obtained in highly symmetric flow fields, such as the Taylor–Green (TG) and Kida–Pelz fields (Reference Boratav and PelzBoratav & Pelz, 1994; Reference Brachet, Meiron, Orszag, Nickel, Morf and FrischBrachet et al., 1983).

Note that the classical potentials of Reference ClebschClebsch (1859) also satisfy the VSF constraint (2.5), but it is much more strict to obtain exact Clebsch potentials than the VSF, and most of the non-trivial Clebsch potentials have singularities (Reference Nore, Abid and BrachetNore, Abid, & Brachet, 1997; Reference Xiong and YangXiong & Yang, 2020). For a flow field with $h\neq 0$, the classical Clebsch potential does not exist, but as illustrated in figure 2(a), the exact VSF can still exist in some special cases, such as the integral Arnold–Beltrami–Childress (ABC) flow (Reference Dombre, Frisch, Greene, Henon, Mehr and SowardDombre et al., 1986).

Furthermore, the exact VSF was obtained in the knotted vorticity field (see figure 2c) with non-trivial topology and local twisting vortex lines (Reference Shen, Yao, Hussain and YangShen, Yao, Hussain, & Yang, 2023; Reference Xiong and YangXiong & Yang, 2019a, Reference Xiong and Yang2020). The complexity of such fields can be characterized by the finite helicity (Reference MoffattMoffatt, 1969) that is the integral of $h$ over a domain bounded by a vortex surface. This construction is related to generalized Clebsch representations, which will be discussed in § 5.

2.3 Numerical construction of the VSF

Since it is generally hard to obtain an exact VSF solution for a flow field with chaotic vortex lines, alternatively, numerical VSF solutions to the VSF constraint (2.5) are sought using two approaches. One is to solve the characteristic equations of $\boldsymbol {\omega }$ such as finding the global first integral of a set of nonlinear ordinary differential equations (Reference Katsanoulis, Kogelbauer, Kaundinya, Ault and HallerKatsanoulis et al., 2023; Reference Yang and PullinYang & Pullin, 2010), which needs developing tailored numerical schemes to deal with ill-posed problems in solving a linear system. The other is to solve (2.5) with an artificial derivative term with respect to a pseudo-time $\tau$ as (Reference Yang and PullinYang & Pullin, 2011a)

(2.6)\begin{equation} \frac{\partial{\phi_{v}}}{\partial{\tau}} + \boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi_{v} = 0. \end{equation}

This advection equation is transported by a frozen vorticity at a particular physical time. The local error of a numerical VSF solution is quantified by

(2.7)\begin{equation} \lambda_{\omega}\equiv \frac{\boldsymbol{\omega}\boldsymbol{\cdot} \boldsymbol{\nabla} \phi_v}{|\boldsymbol{\omega}||\boldsymbol{\nabla}\phi_v|}, \end{equation}

and the volume-averaged VSF error $\epsilon _\phi \equiv \langle |\lambda _\omega | \rangle \rightarrow 0$ has been demonstrated with $\tau \rightarrow \infty$ (Reference Yang and PullinYang & Pullin, 2011a). Namely, $\epsilon _\phi$ converges below a threshold, e.g. $5\,\%$, at a large $\tau$ in solving (2.6).

The isosurface of a numerical VSF solution with small $\epsilon _\phi$ outperforms other vortex identification methods in identifying vortex surfaces. Figure 3 compares visualizations of vortex tubes using the conventional $|\boldsymbol {\omega }|$ and the VSF in homogeneous isotropic turbulence (HIT). In figure 3(a), the isosurface of $|\boldsymbol {\omega }|$ significantly misaligns with the vortex lines, so the tube-like structures do not accurately capture the vortex tubes with large $\epsilon _\phi \approx 30\,\%$. In figure 3(b), the isosurface of $\phi _v$ displays real vortex tubes with small $\epsilon _\phi \approx 5\,\%$.

Figure 3.

Comparison of vortex identification methods: (a) isosurface of $|\boldsymbol {\omega }|$ with $\epsilon _\phi \approx 30\,\%$ and (b) isosurface of the VSF colour coded by $|\boldsymbol {\omega }|$ with $\epsilon _\phi \approx 5\,\%$ in HIT. Some vortex lines are integrated from points on the isosurfaces.

On the other hand, the gradient of $\phi _v$ can be very steep due to the advection nature of (2.6), causing the numerical VSF to be unresolved with finite grid resolution. Thus, the numerical regulation is required to smooth the VSF solution using the weighted essentially non-oscillatory (WENO) (Reference Jiang and ShuJiang & Shu, 1996) and targeted essentially non-oscillatory (Reference Fu, Hu and AdamsFu, Hu, & Adams, 2017) schemes, which can keep the high-order accuracy at the smooth region and meanwhile add an artificial diffusivity at the region with large $\boldsymbol {\nabla } \phi _v$.

There are some additional techniques to keep the smoothness of the VSF solution and accelerate the VSF calculation. The local optimization (Reference Xiong and YangXiong & Yang, 2019b) balances the accuracy and smoothness of numerical VSF solutions by minimizing a hybrid constraint of the VSF error and gradient in the pseudo-evolution. For shear or highly symmetric flows in which vortex lines have a dominant direction passing through a particular plane, such as the streamwise–wall-normal plane in channel/boundary-layer flows and the symmetric plane in TG flows, the VSF constraint (2.5) can be added on the planes to accelerate the convergence in solving VSF solutions (Reference Xiong and YangXiong & Yang, 2017). The combination of the numerical methods above can construct the VSF for any complex flow field.

Moreover, the VSF construction has been applied to experimental data of the tomographic particle image velocimetry (tomo-PIV) (Reference ScaranoScarano, 2012). As shown in figure 4(a), Reference Liu and YangLiu and Yang (2023) used the boundary-constraint method to construct the VSF from the instantaneous tomo-PIV velocity data in the wake flow of a ramp vortex generator (Reference Ye, Schrijer and ScaranoYe, Schrijer, & Scarano, 2016) at moderate $Re$. Despite finite experimental noise and the under-resolved velocity, the VSF construction can have a small error $\epsilon _v\approx 5\,\%$.

Figure 4.

Isosurfaces of the numerical VSF obtained from (a) tomo-PIV data for the flow past a vortex generator and (b) simulation data for the flow past a wing–body configuration. Some vortex lines are integrated from points on the isosurfaces. The surfaces are colour coded by the wall distance in (a) and the vorticity magnitude in (b).

Since the pseudo-initial condition $\phi _{v0}$ for solving (2.6) is arbitrary, the uniqueness is an issue of VSF solutions. For simple flows, the most smooth, time-invariant solution in the velocity field, e.g. the wall distance in the Poiseuille laminar flow (Reference Zhao, Yang and ChenZhao, Yang, & Chen, 2016a) was used, so that the VSF evolution effectively illustrates the structural evolution in a simple flow under small perturbations. For chaotic flows, the statistics of scalar structures of $\phi _v$ appear to be insensitive to the choice of $\phi _{v0}$, perhaps due to the strange attractor of the characteristic equation of $\boldsymbol {\omega }$ (Reference Xiong and YangXiong & Yang, 2017). To avoid the arbitrariness in choosing $\phi _{v0}$, a spatially delta-correlated noisy function was used for $\phi _{v0}$ in the VSF construction for HIT (Reference Xiong and YangXiong & Yang, 2019b).

2.4 Evolution of the VSF

The VSF constructed for an instantaneous flow field can be used as an initial VSF to compute the VSF evolution. The two-time method, combining the physical and pseudo-time evolutions, is effective in calculating the VSF evolution in 3-D non-ideal flows (Reference Hao, Xiong and YangHao et al., 2019; Reference Peng and YangPeng & Yang, 2018; Reference Yang and PullinYang & Pullin, 2011a). Hence, the VSF can serve as a general post-processing tool for analysing a temporally resolved data series of flow fields.

As sketched in figure 5, each time step is divided into two substeps. The prediction–correction approach is similar to the projection method (Reference Kim and MoinKim & Moin, 1985) and the level-set method (Reference Osher and FedkiwOsher & Fedkiw, 2001). In the prediction substep, the temporary VSF $\phi _v^*$ is driven by the physical velocity $\boldsymbol u$ as

(2.8)\begin{equation} \frac{\partial{\phi_v^ * (\boldsymbol{x},t)}}{\partial{t}} + \boldsymbol{u}(\boldsymbol{x},t)\boldsymbol{\cdot}\boldsymbol{\nabla}\phi_v^ * (\boldsymbol{x},t) = 0,\quad t>0. \end{equation}

Figure 5.

Sketch of the two-time method for calculating VSF evolution.

In the correction substep, $\phi ^*_v$ is driven by the frozen vorticity at a fixed physical time $t$ and evolved in pseudo-time $\tau$ as

(2.9)\begin{equation} \frac{\partial{\phi_v^ * (\boldsymbol{x},t;\tau)}}{\partial{\tau}} + \boldsymbol{\omega}(\boldsymbol{x},t) \boldsymbol{\cdot}\boldsymbol{\nabla}\phi_v^ * (\boldsymbol{x},t;\tau) = 0,\quad 0<\tau\leq T_{\tau}, \end{equation}

where the initial condition is $\phi _v(\boldsymbol x,t;\tau =0) = \phi _v^*(\boldsymbol x, t)$. It is a crucial step to dynamically keep the VSF constraint (2.5) satisfied in the VSF evolution. At the end of the correction sub step, $\phi _v(\boldsymbol x,t)=\phi _v^*(\boldsymbol x,t;\tau =T_\tau )$ is updated, where $T_\tau$ is the maximum pseudo-time and is typically less than 100 times ${\rm \Delta} t$.

The numerical schemes for solving (2.8) and (2.9) are described in detail in Reference Yang and PullinYang and Pullin (2011a). Temporal integrations are marched by the second-order total-variation-diminishing Runge–Kutta scheme (Reference Gottlieb and ShuGottlieb & Shu, 1998). The time steps ${\rm \Delta} t$ and ${\rm \Delta} \tau$ were chosen to ensure the corresponding Courant–Friedrichs–Lewy conditions based on $\boldsymbol {u}$ and $\boldsymbol {\omega }$, respectively. The convection terms were typically treated by the fifth-order WENO scheme (Reference Jiang and ShuJiang & Shu, 1996) and the numerical diffusion in the WENO scheme serves as a numerical dissipative regularization to remove small-scale, nearly singular scalar structures. The first VSF evolution was reported in TG flows (Reference Yang and PullinYang & Pullin, 2011a), which reveals a continuous timewise deformation and topological change of vortex surfaces in figure 6. More VSF evolutions will be detailed in § 3.

Figure 6.

Evolution of the VSF isosurface in the TG flow (Reference Yang and PullinYang & Pullin, 2011a) with $Re=400$ at (a) $t=0$, (b) $t=3$, (c) $t=5$ and (d) $t=7$, from initial large-scale vortex blobs to small-scale highly twisted vortex tubes. Attached vortex lines are integrated from the surfaces and colour coded by $|\boldsymbol {\omega }|$.

Additionally, the criteria of tracking vortex surfaces (Reference Han and YangHan & Yang, 2022) contain not only the deviation $\epsilon _\phi$ of the numerical VSF solution, characterizing the accuracy for the VSF constraint, but also a Lagrangian conservation metric based on numerical dissipation, characterizing the time coherence for the one-to-one flow map of VSF isosurfaces. The criteria provide a time range when the vortex-surface tracking is satisfactory.

The computational cost of the VSF evolution grows with increasing complexity of the flow field. For turbulent flows, it can be comparable to or larger than that for the flow solver, which is demanding. To further accelerate the VSF calculation, CPU/GPU heterogeneous computing techniques have been employed. The hybrid method assigns the iteration operations with high computational cost to GPU, and achieves significant improvements on computational efficiency, with acceleration ratios up to 50, in tests on channel flows.

The VSF calculation has also been implemented in complex computational domains on unstructured grids based on the finite volume method with the OpenFOAM open-source software (Reference Weller, Tabor, Jasak and FurebyWeller, Tabor, Jasak, & Fureby, 1998). Since OpenFOAM lacks high-order, discontinuity-resolving schemes due to the challenges on stencil collection and memory management (Reference Gärtner, Kronenburg and MartinGärtner, Kronenburg, & Martin, 2020), explicit artificial diffusivities in the order of viscosity were introduced in the right-hand sides of (2.8) and (2.9) (Reference GuGu, 2022). Thus, the VSF calculation is applicable to general CFD softwares as a flow postprocessing module, which is illustrated in figure 4(b).

3.1 Homogeneous isotropic turbulence

Many past studies believed that the interaction of vortical structures in turbulence leads to the energy transport among different scales (Reference Bermejo-Moreno, Pullin and HoriutiBermejo-Moreno, Pullin, & Horiuti, 2009; Reference Cardesa, Vela-Martin and JimenezCardesa, Vela-Martin, & Jimenez, 2017; Reference DavidsonDavidson, 2004). However, there lacks an effective method to identify complete vortex tubes and then to illustrate an intuitive picture of the energy cascade. The phenomenological models of Reference RichardsonRichardson (1922) and Reference KolmogorovKolmogorov (1941) speculated that the energy cascade from large to small scales is through a process of vortex breakup (Reference FrischFrisch, 1995) such as cell division. By contrast, the VSF visualizations in figures 3 and 7 depict complex networks of tangled vortex tubes.

Figure 7.

Comparison of vortex identifications in HIT. (a) The VSF isosurface in HIT. Close-up views of (b) the isosurface of the VSF and (c) the isosurface of $|\boldsymbol {\omega }|$ in the subdomain enclosed by the dashed rectangle in (a). The isosurface is colour coded by $|\boldsymbol {\omega }|$ and some vortex lines are integrated from points on the surface.

The different vortex identification methods can lead to very different routes of scale cascade. The close-up view of the VSF isosurface with vortex lines integrated on the isosurface in figure 7(b) displays stretched spiral vortex tubes with local vorticity intensification, which is hypothesized as an elementary vortical structure in turbulence (Reference LundgrenLundgren, 1982; Reference Pullin and SaffmanPullin & Saffman, 1998). For comparison, the conventional method, the isosurface of $|\boldsymbol {\omega }|$, in figure 7(c) only displays short, disconnected tube-like structures. This visual breakup fails to display complete vortex surfaces due to the intermittent distribution of $|\boldsymbol {\omega }|$ on the twisting vortex tubes. Therefore, the VSF isosurface can elucidate vortex dynamics such as vortex twisting at small scales (Reference She, Jackson and OrszagShe et al., 1990) and rolling-up of vortex sheets (Reference Lindsay and KrasnyLindsay & Krasny, 2001).

Consistent with the vorticity dynamics, the morphology of VSF isosurfaces at different times in decaying TG and HIT flows (Reference Han and YangHan & Yang, 2022; Reference Xiong and YangXiong & Yang, 2019b; Reference Yang and PullinYang & Pullin, 2011a) illustrates how a large-scale blob-like vortex surface gradually evolves into smaller tube- and sheet-like ones. In figure 6, a possible scenario of scale cascade includes the collapse of vortex surfaces, vortex reconnection, the formation and roll-up of vortex tubes, vorticity intensification between anti-parallel vortex tubes and vortex stretching and twisting.

3.2 Transition of wall-bounded flows

The generation of turbulent structures in the transition of wall-bounded flows was described as ‘abrupt and mysterious’ (Reference MullinMullin, 2011). Existing phenomenological models (Reference BarkleyBarkley, 2016; Reference Barkley, Song, Mukund, Lemoult, Avila and HofBarkley et al., 2015) treated this process as a bifurcation of a nonlinear system, whereas the VSF uncovers the detailed vorticity dynamics on the generating mechanism of coherent structures.

The VSF has been applied to the transition of three canonical wall flows with different symmetries of boundary conditions sketched in figure 8. First, the VSF evolution was calculated in the Klebanoff-type (K-type) natural transition in channel flows using the two-time method. It demonstrates the continuous deformation and topological transformation from simple to convoluted shapes during the transition (Reference Zhao, Yang and ChenZhao et al., 2016a). In the Lagrangian view, figure 9 shows that an initial near-wall planar vortex surface in the laminar state evolves into a triangular bulge, and then the edges of the bulge roll up into the signature hairpin-like structure. Subsequently, the interaction of hairpins in the late transition involves the reconnection of vortex lines originating from opposite walls and the self-reconnection of a single $\varOmega$-shaped vortex line at its neck. These reconnections lead to scale cascade as a critical step in the transition (Reference Zhao, Yang and ChenZhao, Yang, & Chen, 2016b). The reconnection location and time were accurately quantified using the minimum distance and the vorticity flux between two VSF isosurfaces.

Figure 8.

Schematic (front view) of the vortex surfaces in the transition in three canonical wall-bounded flows: (a) boundary layer, (b) channel and (c) pipe. Directions of the vorticity tangent to the vortex surfaces are labelled by arrows.

Figure 9.

Evolution of the VSF isosurface in the late K-type transition in channel flow (Reference Zhao, Yang and ChenZhao et al., 2016b). Attached vortex lines are integrated from the surfaces, and colour coded by $|\boldsymbol {\omega }|$. The evolution can be divided into four stages: (a) formation of the triangular bulge at $t=100$, (b) vorticity intensification at $t=106$, (c) generation of hairpin-like structure at $t=110$ and (d) vortex reconnection at $t=114$. The important events are highlighted by red circles.

Second, the reconnection between vortex bulges grown from opposite walls does not occur in the transition of a flat-plate boundary layer without the symmetry in the wall-normal direction. The hairpins regenerate and then evolve into a complex composite structure, i.e. the turbulent spot (Reference Wu, Moin, Wallace, Skarda, Lozano-Durán and HickeyWu et al., 2017), consisting of multiple hairpins. The VSF was calculated from the direct numerical simulation (DNS) data (Reference Sayadi, Hamman and MoinSayadi, Hamman, & Moin, 2013) of a K-type boundary-layer transition using the boundary-constraint method (Reference Xiong and YangXiong & Yang, 2017). Reference Zhao, Xiong, Yang and ChenZhao, Xiong, Yang, and Chen (2018) found a critical mechanism in lateral growth of the turbulent spot. A train of hairpin-like bulges induce sinuous distortion on neighbouring lower wavy vortex surfaces, and then the spanwise propagation of the distortion generates numerous skewed hairpin structures.

Third, under the axisymmetry of boundary conditions, the VSF was applied to the pipe transition with a radial wave-like velocity disturbance in the inlet region (Reference Ruan, Xiong, You and YangRuan, Xiong, You, & Yang, 2022). The VSF visualization showed that vortex surfaces are first corrugated with streamwise elongated bulges. The resultant highly coiled and stretched vortex loops reconnect with each other, and it then triggers successive vortex reconnections via a ‘greedy snake’ mechanism. The streamwise vortex loops consecutively capture the secondary vortex rings pinched off with self-reconnection, forming long helical vortex loops spanning over 10 pipe radii in the streamwise direction.

Moreover, in the pipe transition with the radial initial perturbation, the transition region has forward and reversed hairpins (Reference Wu, Moin, Adrian and BaltzerWu, Moin, Adrian, & Baltzer, 2015b). The VSF visualization in figure 10 shows that the forward and reversed hairpins are generated from the wall and core regions, respectively (Reference Xiong, You, Ruan and YangXiong, You, Ruan, & Yang, 2019). In contrast to the viscous cancellation of vortex reconnection (Reference Kida and TakaokaKida & Takaoka, 1994) in transitional channel flows, the vorticity directions of approaching forward and reversed hairpin-like vortex lines tend to align (see figure 8). This difference causes stronger vortex interactions, so the pipe transition can be more abrupt than that in transitions of boundary layers and channel flows.

Figure 10.

The side view of the VSF isosurface constructed from the DNS of pipe transition (Reference Wu, Moin, Adrian and BaltzerWu et al., 2015b; Reference Xiong, You, Ruan and YangXiong et al., 2019). The inner VSF isosurface near the core region is colour coded by the normalized radial coordinate $r/R$, and the outer VSF isosurface near the wall is translucent and coloured in grey.

In sum, different symmetries in boundary conditions of the canonical wall flows have an impact on the coherent structures generated in the transition. The VSF studies uncover the emergence of ordered structures and their dynamical evolution in a universal framework, which does not need to objectively switch identification methods and adjust contouring thresholds in different transitional stages.

3.3 Vortex bursting

The VSF is effective in identifying the abrupt deformation of vortex surfaces, e.g. vortex bursting (Reference Ji and van ReesJi & van Rees, 2022; Reference Melander and HussainMelander & Hussain, 1994), which was found in aircraft trailing vortices (Reference TombachTombach, 1973) and could be overlooked by Eulerian vortex criteria. Based on the VSF, Reference Shen, Yao, Hussain and YangShen et al. (2023) developed a framework to construct the differential twist that establishes the theoretical relation between the total twisting number and the local twist rate of each vortex surface. This framework characterizes coiling vortex lines and internal structures within a vortex tube.

As a typical example, the dynamics of a vortex ring with differential twist was explored via DNS. Two twist waves with opposite chiralities propagate towards each other along the centreline of the closed vortex tube and then collide when the local twist rate surges. In figure 11, local vortex surfaces are squeezed into a disk-like structure containing coiled vortex lines, leading to vortex bursting. Note that the visualization of Eulerian vortex criteria, e.g. the visual breakup of isosurfaces of $|\boldsymbol {\omega }|$ similar to figure 7(c), cannot identify the complete vortex tube as visualized by the VSF during bursting.

Figure 11.

The VSF isosurface during the vortex bursting of a vortex ring (Reference Shen, Yao, Hussain and YangShen et al., 2023). The isosurface is colour coded by $h$ and some vortex lines are integrated from points on the surface.

3.4 Flows past a plate

A solid object moving in a fluid flow can generate strong vorticity within coherent structures near the boundary in biological locomotion (Reference DabiriDabiri, 2005; Reference Shyy, Aono, Chimakurthi, Trizila, Kang, Cesnik and LiuShyy et al., 2010; Reference Zhang, Huang, Pan, Yang and HuangZhang, Huang, Pan, Yang, & Huang, 2022) and aerial/underwater vehicles (Reference Jiménez, Hultmark and SmitsJiménez, Hultmark, & Smits, 2010; Reference Park, Chang, Huang and SungPark, Chang, Huang, & Sung, 2014). The morphology of the coherent structures largely depends on the flow direction and the trajectory of the moving body. The VSF was used to investigate the vortex dynamics in flows past a plate (Reference Tong, Yang and WangTong, Yang, & Wang, 2021) simulated by the immersed boundary method (Reference Wang and ZhangWang & Zhang, 2011). In the flow past a stationary finite plate at low $Re$, the VSF isosurface displays that near-plate vortex surfaces first roll up from the plate edges and then evolve into hairpin-like structures near the leading edge and semi-ring structures near the plate tips and in the wake. Two typical vortical structures, the leading edge vortex (LEV) and the tip vortex (TIV), are quantitatively distinguished by the vanishing streamwise vorticity on VSF isosurfaces. Based on circulations through cross-sections of vortex surfaces, the lift generated from the LEV was found to be suppressed by the finite growth of TIVs.

The vortex dynamics of flows past a flapping plate was elucidated using the VSF evolution (Reference Tong, Yang and WangTong et al., 2021). The VSF visualization in figure 12 reveals that the spoon-like vortex surface dominated by tip vortex lines is formed and periodically shed into the wake due to the alternating upstroke and downstroke of the flapping plate. The simplification of the finite-domain impulse theory (Reference Kang, Liu, Su and WuKang, Liu, Su, & Wu, 2018; Reference Wu, Ma and ZhouWu et al., 2015a) demonstrates that the force on the plate can be modelled by the sum of the vortical impulse and Lamb-vector integral of the vortex surface enclosing the plate.

Figure 12.

The VSF isosurface in the flow past a flapping plate (Reference Tong, Yang and WangTong et al., 2021). The isosurface is colour coded by $|\boldsymbol {\omega }|$ and some vortex lines are integrated from points on the surface.

For the moving body with complex geometries, the VSF was employed to study the flow around the NASA Common Research Model (Reference Vassberg, Dehaan, Rivers and WahlsVassberg, Dehaan, Rivers, & Wahls, 2008) in a wing–body configuration (Reference GuGu, 2022). From the simulation data with the Spalart–Allmaras model (Reference Spalart and AllmarasSpalart & Allmaras, 1994), the VSF with $\epsilon _\phi < 10\,\%$ was calculated using the two-time method. The VSF isosurface around the wing–body configuration in figure 4(b) shows that most vortices develop along the streamwise direction, forming bulges similar to those in canonical wall flows (see figure 9a). They gradually evolve into small-scale structures and eventually shed off from the wing.

3.5 Compressible flows

Besides the constant-density incompressible flows in fundamental studies, the VSF was also used to analyse various multi-physics coupled flows in practical flows. Due to the characteristic structure arising from additional physical or chemical processes in such flows, the coupling mechanism of these structures and vortices is challenging in the modelling of engineering turbulence. Based on the VSF study in incompressible TG flows (Reference Yang and PullinYang & Pullin, 2011a), evolutions of vortex surfaces have been compared in compressible, reacting or MHD TG flows to elucidate the interaction mechanisms between the vortex surface and the shock wave, flame or electromagnetic field.

In compressible hydrodynamic instability (Reference BrouilletteBrouillette, 2002) and turbulence (Reference Samtaney, Pullin and KosovicSamtaney, Pullin, & Kosovic, 2001; Reference Wang, Shi, Wang, Xiao, He and ChenWang et al., 2012), the effects of the Mach number ($Ma$) on the vortex dynamics are important. For the strong compressibility at large $Ma$, the VSF study found that the vortex dynamics in high-$Ma$ TG flows has three phases (Reference Peng and YangPeng & Yang, 2018). In the early stage, the velocity jump near shocklets generates sinks to shrink vortex volume and distort vortex surfaces. The subsequent vortex reconnection occurs earlier and has a higher reconnection degree for larger $Ma$. In the late stage, the positive dissipation rate and negative pressure work accelerate the loss of kinetic energy and suppress vortex twisting with increasing $Ma$.

The VSF was then applied to the Richtmyer–Meshkov instability (RMI) accelerated by a weak incident shock. The RMI is important in the supernova explosion (Reference Almgren, Bell, Rendleman and ZingaleAlmgren, Bell, Rendleman, & Zingale, 2006) and ignition in inertial confinement fusion (ICF) (Reference Aglitskiy, Velikovich, Karasik, Metzler, Zalesak, Schmitt and ObenschainAglitskiy et al., 2010). Reference Peng, Yang, Wu and XiaoPeng, Yang, Wu, and Xiao (2021a) elucidated the effect of the secondary baroclinic vorticity (SBV) in the RMI. Two major mechanisms of the single-mode RMI, the primary baroclinic vorticity (PBV) and the pressure perturbation, are distinguished by VSF-based models. They found that the interfacial growth is first driven by the PBV. Subsequently, the SBV is generated by the misalignment between the density gradient across the interface and the pressure gradient produced by the PBV-induced velocity. Driven by the induced velocities of the PBV and SBV, the acceleration of the trough and deceleration of the crest of the vortex surface lead to spike- and bubble-like structures.

Furthermore, Reference Peng, Yang and XiaoPeng, Yang, and Xiao (2021b) investigated the SBV effect on the energy cascade in the mixing induced by the multi-mode RMI (see figure 13a). With the VSF-based models, they found that the effect of the SBV peaks at a critical time when the vortex reconnection widely occurs in the mixing zone. Before the critical time, spikes and bubbles evolve almost independently, and the variation of the kinetic energy spectrum induced by the SBV has the $-1$ scaling law at intermediate wavenumbers. This SBV effect causes the slope of the total energy spectrum at intermediate wavenumbers to evolve towards $-$3/2 at the critical time. Subsequently, the SBV effect diminishes and the energy spectrum decays to the $-$5/3 law.

Figure 13.

The VSF isosurfaces in multi-physics coupled flows: (a) multi-mode RMI after a shock accelerating an interface (Reference Peng, Yang and XiaoPeng et al., 2021b) and (b) premixed flame (location marked by the red arrow) propagating in HIT (Reference You, Lu and YangYou et al., 2020). The isosurface is colour coded by $|\boldsymbol {\omega }|$ and some vortex lines are integrated from points on the surface.

3.6 Combustion

For combustion, the flame–vortex interaction is important in various combustion applications (Reference Renard, Thévenin, Rolon and CandelRenard, Thévenin, Rolon, & Candel, 2000), but its study was restricted to 2-D configurations (Reference Poinsot, Veynante and CandelPoinsot, Veynante, & Candel, 1991). The VSF study illustrated the influence of a propagating premixed flame on 3-D TG vortices (Reference Zhou, You, Xiong, Yang, Thévenin and ChenZhou et al., 2019). The vortex surfaces merge into a bulky structure in the burnt region, whereas they roll up into stretched and twisted small-scale vortex tubes in the unburnt region, as in non-reacting flows. This finding suggests that the heat release near a turbulent premixed flame significantly alters the fluid density and viscosity, leading to anisotropic velocity statistics.

As shown in figure 13(b) for turbulent premixed combustion, Reference You, Lu and YangYou, Lu, and Yang (2020) found that a tangle of twisted vortex tubes on the unburnt side is stretched along the streamwise direction near the flame front due to the thermal expansion, and the small-scale vortex tubes gradually merge into large-scale bulky structures on the burnt side. By analysing the enstrophy transport equation and the Lumley triangle of Reynolds stresses, they concluded that the anisotropy of velocity fields increases from the unburnt to the burnt side. The variation of the geometry of vortex surfaces near the flame front is highly correlated to the anisotropic statistics of the fluctuating velocity.

3.7 Magnetohydtodynamic flows

The transport equations for the vorticity in neutral fluid flows and the magnetic induction in MHD flows have the same form (Reference DavidsonDavidson, 2001). Due to the analogy between vorticity and magnetic dynamics, the VSF can be naturally extended to analyse the magnetic surface as the magnetic-surface field (MSF) (Reference Hao and YangHao & Yang, 2021). The MSF is rooted in the Alfvén theorem – an analogue of the Helmholtz vorticity theorem to illustrate the frozen-in nature of magnetic fields. The MSF study characterized the evolution of a pair of orthogonal helical flux tubes (Reference Linton, Dahlburg and AntiochosLinton, Dahlburg, & Antiochos, 2001) with opposite chirality. In particular, the MSF can pinpoint the region of the magnetic reconnection (Reference Priest and ForbesPriest & Forbes, 2000), where a knot cascade of magnetic field lines through stepwise reconnections was identified.

Additionally, the Lorentz force in MHD flows behaves as a strong non-ideal force in (2.3). Its magnitude can be comparable to that of the convection term in practical MHD flows. Reference Hao, Xiong and YangHao et al. (2019) found the spurious deformation due to the Lorentz force can be eliminated using an approximate virtual velocity. The vortex surfaces driven by the virtual velocity accelerates the vortex reconnection and significantly increases the dissipation of kinetic energy.

3.8 Remarks on the VSF visualization

Besides the time coherence issue discussed in § 2.2, the geometric features identified in turbulence and transition strongly depend on the identification method and the isocontour threshold. These selections have obvious subjectivity and cause many debates on coherent structures in past decades (Reference EppsEpps, 2017; Reference WuWu, 2018). In boundary-layer transition (Reference Sayadi, Hamman and MoinSayadi et al., 2013; Reference Wu, Moin, Wallace, Skarda, Lozano-Durán and HickeyWu et al., 2017), the isosurface of the streamwise velocity or passive scalar can only show streaky structures, and the isosurface of the Eulerian vortex criterion can only show vortex cores. As compared among the isosurfaces of the streamwise velocity, the $Q$ criterion and VSF in a boundary-layer transition in figure 13 in Reference Zhao, Xiong, Yang and ChenZhao et al. (2018), isocontour thresholds of the streamwise velocity and $Q$ have to be carefully adjusted to capture complete vortex surfaces. Therefore, the whole picture of the vortex dynamics in transition, e.g. interactions of streaks and hairpin vortices (Reference Wu, Moin, Wallace, Skarda, Lozano-Durán and HickeyWu et al., 2017), may not be fully represented using a single classical method.

In the VSF framework, the very different coherent structures can be considered as different geometries of the same vortex surfaces at different times and locations. For example, the laminar shear flow corresponds to wall-parallel planar surfaces, the $\varLambda$-like structure corresponds to the triangular bulge from the wall-normal disturbed vortex surface, the hairpin corresponds to the hairpin bulge from rolling up and stretching of the triangular bulge and the streak corresponds to the wavy surface downstream of the hairpin ‘legs’.

Additionally, the fixed VSF isocontour threshold, in principle, identifies a particular vortex surface in the evolution in an ideal flow. This guides that the VSF should be normalized to ensure the fixed VSF isocontour threshold for the isosurface which encloses the same mass in the evolution of the VSF in viscous flows (Reference Peng and YangPeng & Yang, 2018), which mitigates the numerical dissipation in the regularization in solving (2.8) and (2.9).

Note that the VSF is space filling. The selection of the VSF threshold follows several criteria. For the flows without a dominant direction of $\boldsymbol {\omega }$, e.g. TG flow and HIT, the median value of $\phi _v$ can be used to show the vortex surface with both large volume and distinguished geometry. For transitional wall flows, the morphology of vortex surfaces varies with the initial wall distance, and the VSF threshold is proportional to the wall distance in the laminar stage (Reference Zhao, Yang and ChenZhao et al., 2016a). The VSF isosurface near to or remote from the wall with a very small or large threshold has minor deformation, while the VSF isosurface with the threshold $O(0.1)\sim O(1)$ of the boundary-layer thickness have significant deformation to characterize the evolution of coherent structures (Reference Zhao, Xiong, Yang and ChenZhao et al., 2018, Reference Zhao, Yang and Chen2016b).

4.1 Friction drag

The prediction of the skin-friction coefficient $C_f$ in boundary layers is crucial for the design of high-speed vehicles and propulsion systems. The VSF study found that $C_f$ grows with the deformation and reconnection of vortex surfaces in transitional wall flows (Reference Zhao, Xiong, Yang and ChenZhao et al., 2018, Reference Zhao, Yang and Chen2016b). This correlation supports the development of structure-based models of $C_f$. For practical interest, the VSF deviation of a passive scalar/tracers used in experimental visualization, e.g. dye, smoke and hydrogen bubbles (Reference Lee and WuLee & Wu, 2008), can be calculated by (2.7). If $\epsilon _\phi$ is small, the VSF study can be extended to experimental visualization, and has a potential to improve the vorticity reconstruction from the scalar field (Reference Liu and ShenLiu & Shen, 2007; Reference Su and DahmSu & Dahm, 1996).

The multi-scale geometric analysis was applied to tracking a passive scalar (Reference Candes, Demanet, Donoho and YingCandes, Demanet, Donoho, & Ying, 2006; Reference Yang, Pullin and Bermejo-MorenoYang, Pullin, & Bermejo-Moreno, 2010). From the band-pass filtering, the geometric features of the scalar were quantified at different scales. Reference Yang and PullinYang and Pullin (2011b) found the similarity of inclination and sweeping angles of coherent structures in wall turbulence, and statistical correlation between the growths of the inclination angle and $C_f$ (Reference Zheng, Yang and ChenZheng, Yang, & Chen, 2016). As shown in figure 14, a model for predicting $C_f$ (Reference Zheng, Ruan, Yang, He and ChenZheng, Ruan, Yang, He, & Chen, 2019) was proposed from multi-scale structural geometry in 2-D images in supersonic boundary-layer transition (Reference He, Yi, Zhao, Tian and ChenHe, Yi, Zhao, Tian, & Chen, 2011). In this model, the averaged inclination angle is a typical geometric feature of vortical structures. This feature is robust in a range of $Re$ and $Ma$ and is correlated to the intermittency factor (Reference Dhawan and NarasimhaDhawan & Narasimha, 1958) in turbulence modelling (Reference Suzen and HuangSuzen & Huang, 2000), so it blends the empirical models of $C_f$ in laminar and fully developed turbulent states (Reference WhiteWhite, 2006) to estimate $C_f$ in the transitional stage.

Figure 14.

Predictive model of the skin-friction coefficient based on the multi-scale decomposition of experimental images in supersonic boundary-layer transition (Reference Zheng, Ruan, Yang, He and ChenZheng et al., 2019): (a) inclination angle of small-scale structures after multi-scale decomposition and (b) friction estimation based on the inclination angle. Here, $l_s$ denotes the characteristic scale in terms of the boundary-layer thickness $\delta$; $c_{fL}$ and $c_{fT}$ denote empirical models of $C_f$ in laminar and fully developed turbulent states (Reference WhiteWhite, 2006), respectively.

4.2 Aerodynamic forces

Shedding vortices in the far wake, e.g. in the wake of a flying bird (Reference Hubel, Riskin, Swartz and BreuerHubel, Riskin, Swartz, & Breuer, 2010) or swimming fish (Reference LauderLauder, 2015), can imply the forces exerted on the moving body (Reference Li and LuLi & Lu, 2012). Based on the VSF study and the vortex impulse theory (Reference Wu, Liu and LiuWu, Liu, & Liu, 2018), Reference Tong, Yang and WangTong et al. (2021) developed a time-averaged thrust model from near-wake discrete vortex surfaces. In this model, a particular vortex surface is chosen as the integral domain in the finite-domain impulse theory, so that the force on the moving body is only determined by the vortical impulse and Lamb-vector integral of the vortex surface enclosing the body. Furthermore, the mean thrust based on two arbitrary vortex surfaces in the far wake is estimated from the linear impulse decay of periodically shedding vortex surfaces.

The VSF study suggested that the state of a moving body in experimental investigation and practical applications can be inferred from wake structures. Reference Tong, Wang and YangTong, Wang, and Yang (2022) reported a comparative study of theoretical and data-driven models for estimating forces from velocity data in the wake of a plate with a range of angles of attack. The theoretical model estimates forces on a flat plate from cross-sectional velocity data in the far wake. This algebraic model incorporates the local momentum deficit and pressure variation. The data-driven model estimates forces from the convolutional neural network (CNN) (Reference He, Zhang, Ren and SunHe, Zhang, Ren, & Sun, 2016) by regarding the velocity field on a series of cross-sections as images. The model performance indicated that the optimized CNN with the integration of physical information can identify important flow regions and learn empirical physical laws.

4.3 Turbulent mixing

The mixing of the interface induced by the RMI (Reference ZhouZhou, 2017a, Reference Zhou2017b), quantified by the peak-to-valley amplitude of the mixing zone, is particularly important in the ICF. Inspired by the SBV mechanism from the VSF study on the RMI, Reference Peng, Yang, Wu and XiaoPeng et al. (2021a) developed a predictive model of spike and bubble growth rates using the motion of viscous vortex rings. The circulation of the vortex ring is modelled with the SBV effect. The modelling process is sketched in figure 15. The feature of almost independently evolving spikes and bubbles is characterized by the vortex rings with opposite vorticity-induced velocities. The motion of each vortex ring is modelled by the dual polar coordinates moving with the vortex ring. This model was validated by several DNS and experimental data sets of the single-mode RMI with various initial conditions.

Figure 15.

Schematic for the modelling of spikes and bubbles in the multi-mode RMI (Reference Peng, Yang and XiaoPeng et al., 2021b). (a) Close-up view of spikes and bubbles; (b) spikes and bubbles moving by the induced velocities of vortex rings; (c) vortex-ring modelling.

This model was then extended to multi-mode RMI to estimate the mixing width (Reference Peng, Yang and XiaoPeng et al., 2021b), and was validated using numerical simulations of the RMI with different modes of initial perturbations. It captures the nonlinear growth of the mixing width before the self-similar growth stage.

5.1 Spherical Clebsch representation

The Clebsch representation encodes the velocity and vorticity into multiple scalar potentials $(\varphi _1,\varphi _2,\varphi _3)$ as

(5.1a,b)\begin{align} \boldsymbol u = \varphi_1 \boldsymbol \nabla \varphi_2 - \boldsymbol \nabla \varphi_3,\quad \boldsymbol \omega = \boldsymbol \nabla \varphi_1 \times \boldsymbol \nabla \varphi_2, \end{align}

where $\varphi _1$ and $\varphi _2$ satisfy the VSF constraint (2.5), so this Clebsch map contains the flow geometry. However, the classical Clebsch map in (5.1a,b) is unable to describe a flow field with non-zero helicity, restricting its application in general 3-D flows.

Reference Kuznetsov and MikhailovKuznetsov and Mikhailov (1980) proposed the spherical Clebsch map (Reference Chern, Knöppel, Pinkall and SchröderChern et al., 2017, Reference Chern, Knöppel, Pinkall, Schröder and Weißmann2016) to express the flow field with non-zero helicity, i.e. the velocity

(5.2)\begin{equation} \boldsymbol u = \hbar \mbox{Re}(\boldsymbol \nabla\bar{\boldsymbol \psi}\mathrm{i} \boldsymbol \psi) = \hbar(a\boldsymbol \nabla b-b \boldsymbol \nabla a + c\boldsymbol \nabla d-d\boldsymbol \nabla c) \end{equation}

is expressed by a quaternion $\boldsymbol \psi =a+b\mathrm {i} +c\mathrm {j}+d\mathrm {k}$ with a normalization constraint $\bar {\boldsymbol \psi }\boldsymbol \psi =1$, where $a$, $b$, $c$ and $d$ are real-valued functions, $\mathrm {i}$, $\mathrm {j}$ and $\mathrm {k}$ are imaginary units, $\hbar$ is a parameter for quantization of $\boldsymbol u$ and $\mbox {Re} (\boldsymbol \psi )$ and $\bar {\boldsymbol \psi }$ are the real part and the complex conjugate of $\boldsymbol \psi$, respectively. The vorticity Clebsch map, also known as the Bloch or spin vector (Reference BlochBloch, 1946), is expressed as $\boldsymbol s = s_1\mathrm {i} +s_2\mathrm {j}+s_3\mathrm {k}$ and linked with $\boldsymbol \psi$ by the Hopf fibration $\boldsymbol s =\bar {\boldsymbol \psi }\mathrm {i} \boldsymbol \psi$ (Reference HopfHopf, 1931). Thus, the vorticity is re-written as

(5.3) \begin{equation} \boldsymbol \omega = \frac{\hbar}{2}(s_{1}\boldsymbol \nabla s_{2}\times \boldsymbol \nabla s_{3} + s_{2} \boldsymbol \nabla s_{3}\times\boldsymbol \nabla s_{1} + s_{3}\boldsymbol \nabla s_{1}\times\boldsymbol \nabla s_{2}), \end{equation}

in terms of $\boldsymbol s$, where the Clebsch potentials $s_i,\ i=1,2,3$ are VSFs.

Figure 16 summarizes the transformations among four representations of $\boldsymbol \psi$, $\boldsymbol u$, $\boldsymbol \omega$ and $\boldsymbol s$ for the same vortex-ring field in different spaces. The two-component or quaternion wave function (upper left) in $\mathbb {S}^3$ is represented by a complex structure that cannot be directly visualized in $\mathbb {R}^3$. The spin vector (upper right) in $\mathbb {S}^2$ can be obtained via the Hopf fibration from $\boldsymbol \psi$. The velocity (lower left) in $\mathbb {R}^3$ can be calculated from a given wave function by (5.2), which describes a vortex ring (green tube) with circular stream lines (brown lines). The vorticity (lower right) can be obtained from $\boldsymbol u$ by its definition or the spin vector by (5.3).

Figure 16.

Transformations among the wave function $\boldsymbol \psi$, velocity $\boldsymbol u$, vorticity $\boldsymbol \omega$ and spin vector $\boldsymbol s$ for a vortex ring in different spaces.

Similar to the VSF construction, the construction of Clebsch potentials for a given flow field is challenging (Reference Chern, Knöppel, Pinkall, Schröder and WeißmannChern et al., 2016; Reference Xiong, Wang, Wang and ZhuXiong, Wang, Wang, & Zhu, 2022; Reference Yang, Xiong, Zhang, Feng, Liu and ZhuYang et al., 2021). Reference Tao, Ren, Tong and XiongTao, Ren, Tong, and Xiong (2021) developed an analytical method to construct $\boldsymbol \psi$ that can be transformed into a knotted velocity field. Taking the flow field of a trefoil vortex for example, $\boldsymbol \psi$ is constructed by two polynomials $P = \alpha ^k$ and $Q = \alpha ^3+\beta ^2$ representing the geometry of the vortex tube (Reference Kedia, Foster, Dennis and IrvineKedia, Foster, Dennis, & Irvine, 2016) with $\alpha =2(x+\mathrm {i} y)f(R)/(1+R^{2})$ and $\beta =\mathrm {i} +2(z-\mathrm {i} )f(R)/(1+R^{2})$, where $R=\sqrt {x^{2}+y^{2}+z^{2}}$ is the distance from the knot centre and the function $f(R) = \exp [-(R/3)^8]$ has $f(0)=1$ and $\textrm {lim}_{R\rightarrow \infty }f(R)=0$. Then, $\boldsymbol \psi$ is obtained from the normalization and pressure projection of $P+Q\mathrm {j}$. Figure 17 plots the isosurfaces of $s_1$ of two trefoil vortex tubes with different contour and helicity values.

Figure 17.

Isosurfaces of $s_1$ of initial knotted vortex tubes with (a) $P = \alpha ^1$, $Q = \alpha ^3+\beta ^2$, $s_1 = 0.92$ and (b) $P = \alpha ^5$, $Q = \alpha ^3+\beta ^2$, $s_1 = 0.25$. Some vortex lines colour coded by $h$ are plotted on the surfaces.

5.2 Computer graphics of fluid motion

In the Clebsch method, the evolution of $\boldsymbol u$ is transformed into the evolution of $\boldsymbol \psi$ using (5.2). For an incompressible, constant-density inviscid flow with $\boldsymbol F = 0$ in (2.2), the momentum conservation (2.2) becomes (Reference Xiong, Wang, Wang and ZhuXiong et al., 2022)

(5.4)\begin{equation} \frac{\partial \boldsymbol \psi}{\partial t} + \boldsymbol u \boldsymbol{\cdot} \boldsymbol \nabla \boldsymbol \psi = \frac{\mathrm{i}}{\hbar}\left(\varPi + \frac {|\boldsymbol{u}|^2}{2}\right)\boldsymbol \psi, \end{equation}

and the incompressible condition is equivalent to

(5.5)\begin{equation} \mbox{Re}( \nabla^2\bar{\boldsymbol \psi}\mathrm{i} \boldsymbol \psi) = 0. \end{equation}

Reference Yang, Xiong, Zhang, Feng, Liu and ZhuYang et al. (2021) developed a fast algorithm to solve (5.4) and (5.5) with notable numerical viscosity, which is suitable for fluid animation in computer graphics. Similar to the stable fluid method (Reference StamStam, 1999), a staggered grid was adopted to discretize fluid variables and solve the equations using the advection–projection splitting method (Reference Kim and MoinKim & Moin, 1985). The grid stores $\boldsymbol \psi$ and the gauge variable $\varPi$ at the centre of each cell, and $\boldsymbol u$ at the centre of each cell face. As illustrated in figure 16, the wave function extracts vortex lines and surfaces analytically, so the Clebsch method can preserve the vortical structure at long times in the flow simulation with lower grid resolution and computational cost than conventional CFD methods.

The advection of $\boldsymbol \psi$ inherits the chaotic behaviour as evolving a Lagrangian-based VSF (Reference Yang and PullinYang & Pullin, 2011b; Reference Yang, Pullin and Bermejo-MorenoYang et al., 2010). To maintain the stability of the advection, Reference Yang, Xiong, Zhang, Feng, Liu and ZhuYang et al. (2021) used a blending method to add the numerical viscosity for $\boldsymbol u$ into the evolution of $\boldsymbol \psi$. After the advection, the divergence-free conditions on $\boldsymbol u$ and $\boldsymbol \psi$ are both enforced by solving the Poisson equation

(5.6)\begin{equation} {\rm \Delta} t \nabla^2 \varPi = {-}\boldsymbol \nabla \boldsymbol{\cdot} \boldsymbol u, \end{equation}

for $\varPi$ in (5.4) and projecting $\boldsymbol u +{\rm \Delta} t\boldsymbol {\nabla } \varPi \to \boldsymbol u$ and $\boldsymbol \psi \exp (\mathrm {i}\varPi {\rm \Delta} t/ \hbar ) \to \boldsymbol \psi$, respectively. The blending and projection introduce the numerical viscosity generated in the semi-Lagrangian scheme (Reference StamStam, 1999) into the advection of $\boldsymbol \psi$, and the real viscosity is expected to be incorporated in future work.

Reference Yang, Xiong, Zhang, Feng, Liu and ZhuYang et al. (2021) demonstrated that the Clebsch method can resolve the fluid motion with complex vortical structures and boundary conditions on coarse grids. In smoke simulations, figure 18(a) shows that arrays of coherent vortices are generated after a high-speed smoke plume passes an obstacle, and figure 18(b) shows the reconnection of two vortex rings in an oblique collision. These examples with fine vortex-filament dynamics were difficult to simulate using a standard grid-based solver with the same grid resolution.

Figure 18.

Visualizations of flow simulations based on the Clebsch method (Reference Yang, Xiong, Zhang, Feng, Liu and ZhuYang et al., 2021). (a) Smoke passing multiple obstacles with the grid size $192 \times 192 \times 384$; (b) evolution of oblique ring collision with the grid size $128 \times 128 \times 128$.

The Clebsch method has been extended to simulate free-surface vortical flows (Reference Xiong, Wang, Wang and ZhuXiong et al., 2022; Reference Yang, Xiong, Zhang, Feng, Liu and ZhuYang et al., 2021). Assuming the fluid is subject to gravity and surface tension and according to the Young–Laplace equation (Reference BatchelorBatchelor, 1967), $\varPi$ in (5.4) satisfies a jump condition on $\varGamma$. The interface is tracked by a level-set function (Reference Fedkiw and OsherFedkiw & Osher, 2002; Reference Osher and SethianOsher & Sethian, 1988), and (5.6) with the jump condition is solved by the ghost fluid method (Reference Fedkiw, Aslam, Merriman and OsherFedkiw, Aslam, Merriman, & Osher, 1999).

In addition, the Clebsch method is efficient in simulating a broad range of free-surface flows with strong vortices, such as bubble rings, horseshoe vortices, sink swirls and free-surface wake vortices. Figure 19(a) shows the vortex shedding from a moving paddle on a coarse grid, and figure 19(b) shows the reconnection of two bubble rings. By contrast, the vortices may rapidly dissipate with the regular finite difference methods, followed by the fragmentation of the interfacial vortices and the bubble ring.

Figure 19.

Visualizations of flow simulations based on the Clebsch method (Reference Xiong, Wang, Wang and ZhuXiong et al., 2022). (a) Vortex shedding from a moving paddle with the grid size $128\times 32 \times 64$; (b) reconnection of bubble rings with the grid size $256 \times 256 \times 256$.

5.3 Quantum computing of fluid dynamics

Quantum computing has emerged to be the next disruptive technology, which can have much less computational cost than the traditional computing (Reference Nielsen and ChuangNielsen & Chuang, 2010). To date, quantum computing has been demonstrated to be effective in handling some linear problems (Reference Clader, Jacobs and SprouseClader, Jacobs, & Sprouse, 2013; Reference Harrow, Hassidim and LloydHarrow, Hassidim, & Lloyd, 2009), but it remains challenging to efficiently solve the highly nonlinear NS equations (Reference Liu, Kolden, Krovi, Loureiro, Trivisa and ChildsLiu et al., 2021; Reference Lloyd, Palma, Gokler, Kiani, Liu, Marvian and PalmerLloyd et al., 2020).

The fluid dynamics of a potential flow was related to a hydrodynamic representation of quantum mechanics via the Madelung transform (Reference MadelungMadelung, 1927). By generalizing the Madelung transform (Reference Chern, Knöppel, Pinkall, Schröder and WeißmannChern et al., 2016; Reference Meng and YangMeng & Yang, 2023), the compressible/incompressible flow with finite vorticity and dissipation can be described by the hydrodynamic Schrödinger equation (HSE)

(5.7)\begin{equation} \mathrm{i}\hbar\frac{{\partial}\boldsymbol{\psi}}{{\partial} t} = \left(-\frac{\hbar^2}{2}\nabla^2 + \varPi + \frac{\hbar^2}{8\rho^2}|\boldsymbol{s}|^2 \right)\boldsymbol{\psi}. \end{equation}

This equation is equivalent to (2.1) and (2.2) with an external force $\boldsymbol F = - ({\hbar }^2/4\rho )\boldsymbol {\nabla }\boldsymbol {s}\boldsymbol {\cdot }[\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\nabla }\boldsymbol {s}/\rho )]$. The corresponding flow satisfies the virtual Helmholtz theorem (Reference Hao, Xiong and YangHao et al., 2019) and has helicity conservation (Reference Meng, Shen and YangMeng, Shen, & Yang, 2023).

The quantum computing of fluid dynamics based on the HSE can be promising in simulating 3-D turbulent flows. Since the HSE is expressed as a unitary operator on $\boldsymbol \psi$, it is more suitable than the NS equations for quantum computing. The flow governed by the HSE can resemble a turbulent flow consisting of tangled vortex tubes with the $-5/3$ scaling of energy spectrum. A prediction–correction quantum algorithm for solving the HSE was implemented for a simple one-dimensional flow on the quantum simulator Qiskit (Reference Koch, Wessing and AlsingKoch, Wessing, & Alsing, 2019).

In particular, the incompressible HSE with unity density governs the incompressible Schrödinger flow (ISF) (Reference Chern, Knöppel, Pinkall, Schröder and WeißmannChern et al., 2016). Although the HSE without a viscous term and with an external force is different from the NS equations, the ISF resembles the viscous flow in terms of the similar flow statistics and structures. The evolution of the vortex knot, TG vortex, decaying HIT showed the similarities between the ISF and the viscous flow (Reference Meng and YangMeng & Yang, 2023; Reference Tao, Ren, Tong and XiongTao et al., 2021).

The effect of the parameter $\hbar$ in (5.7) is similar to the kinematic viscosity. In general, the length scale of vortices is proportional to $\hbar$ via the vorticity Clebsch mapping (Reference Chern, Knöppel, Pinkall and SchröderChern et al., 2017, Reference Chern, Knöppel, Pinkall, Schröder and Weißmann2016; Reference Tao, Ren, Tong and XiongTao et al., 2021), and the flow stability depends on the value of $\hbar$. The vortex surface in the TG ISF is visualized in figure 20 using the isosurface of a component of $\boldsymbol s$. The entangled vortex tubes can be mapped to a patch on the unit sphere $\mathbb {S}^2$ (or the Bloch sphere) via the vorticity Clebsch mapping. The geometry of vortex surfaces in the ISF is in between the vortex filaments in quantum turbulence (Reference Madeira, Caracanhas, dos Santos and BagnatoMadeira, Caracanhas, dos Santos, & Bagnato, 2020; Reference Müller, Polanco and KrstulovicMüller, Polanco, & Krstulovic, 2021) and the tangle of spiral vortex tubes and sheets in classical turbulence (Reference Cardesa, Vela-Martin and JimenezCardesa et al., 2017; Reference She, Jackson and OrszagShe et al., 1990; Reference Xiong and YangXiong & Yang, 2019b). Therefore, the turbulent ISF manifests the features of both quantum and classical turbulent flows.

Figure 20.

Mapping between the spin vector and the vortex tube in the TG ISF. (a) A patch (blue) on $\mathbb {S}^2$ of the spin vector, where each dot corresponds to a vortex line in $\mathbb {R}^3$; (b) isosurface of a component of the spin vector in $\mathbb {R}^3$, where the surface is colour coded by $h$; (c) close-up view of the vortex tubes with attached vortex lines.