Applications of the vortex-surface field to flow visualization, modelling and simulation

Abstract

We review the progress on the applications of the vortex-surface field (VSF). The VSF isosurface is a vortex surface consisting of vortex lines. Based on the generalized Helmholtz theorem, the VSF isosurfaces of the same threshold at different times have strong coherence. As a general flow diagnostic tool for studying vortex evolution, the numerical VSF solution is first constructed in a given flow field by solving a pseudo-transport equation driven by the instantaneous frozen vorticity, and then the VSF evolution is calculated by the two-time method. From the database of numerical simulations or experiments, the VSF can elucidate mechanisms in the flows with essential vortex dynamics, such as isotropic turbulence, wall flow transition, flow past a flapping plate and turbulence–flame interaction. The characterization of VSFs reveals the correlation between robust statistical features and the critical quantities needed to be predicted in engineering applications, such as the friction coefficient in transition, thrust in bio-propulsion and growth rate in interface instability. Since the VSF evolution captures the essential Lagrangian-based dynamics of vortical flows, it inspires novel numerical methods on cutting-edge hardware, e.g. graphic and quantum processors.

Impact Statement

Vortex dynamics plays a crucial role in the fundamental problems of turbulence and transition and in the engineering applications of high-speed flight, combustion and bio-inspired robotics. On the other hand, there is no consensus on the vortex identification, and there is a lack of a general framework to convert the qualitative observation of coherent structures into predictive models. The vortex-surface field serves as a general flow diagnostic tool. It paves an avenue to investigate the mechanism of energy cascade and the evolution of coherent structures in turbulence and transition from the Lagrangian viewpoint, and inspires a series of simple models and novel simulation methods in engineering applications.

1. Introduction

We review the progress on the theory and applications of the vortex-surface field (VSF). The VSF serves as a general flow diagnostic tool and inspires Lagrangian-based models and simulation methods. In particular, the VSF can be used to tackle the critical issues of vortex dynamics in turbulence and transition.

The fluid motion is governed by the Navier–Stokes (NS) equations. The high nonlinearity in the equations causes the turbulent flow to manifest both order and randomness, e.g. large-scale coherent structures (Adrian, 2007; Lee & Wu, 2008; Robinson, 1991) in shear turbulence, and small-scale velocity fluctuations satisfying the Gaussian distribution remote from the wall (Pope, 2000). The co-existence of ordered and random natures makes the turbulence problem notoriously difficult.

The studies on turbulence can be roughly classified into two schools based on structures and statistics. The former one prefers analysis from certain equations, such as the dynamic system (Holmes, Lumley, Berkooz, & Rowley, 2012), vortex model/dynamics (Pullin & Saffman, 1998) and topological fluid dynamics (Moffatt, 2021); the latter one prefers statistical analysis from stochastic processes, such as non-equilibrium statistical physics (Frisch, 1995; He, Jin, & Yang, 2017; Sreenivasan, 1999). Both turbulence schools have achieved remarkable progresses, but meanwhile they have obvious weaknesses. The structure-based studies can provide qualitative and quantitative descriptions (Lozano-Durán & Jiménez, 2014; Robinson, 1991) but lack a universal approach to developing predictive models. The statistical theory reveals statistical features of small-scale turbulence but lacks the information of large-scale motions in practical flows. Therefore, the present turbulence theory can only deal with relatively ideal cases of turbulence, without a theoretical framework fusing structural and statistical studies.

The transition is another classical challenging problem in fluid dynamics (Eckhardt, Schneider, Hof, & Westerweel, 2007; Mullin, 2011; Reynolds, 1883), and it is also a critical issue in aerospace applications. In high-speed flight, the flow on the aerospace vehicle has a rapid transition from the laminar to turbulent state (Fedorov, 2011; Zhong & Wang, 2012). The emergence of coherent structures often accompanies the surge of the skin friction and aerodynamic heating (Lee & Jiang, 2019), causing bumps and ablation on aero-vehicles (Leyva, 2017). In such complex transitional processes, it is crucial to effectively identify the coherent structures and their evolution, and accurately characterize their impact on the transport of momentum and energy. On the other hand, the direct measurement of aerodynamic forces and heating becomes extremely difficult in the experiments of supersonic/hypersonic transition in wall-bounded flows (Hakkinen, 2004; Schetz, 2010; Zhu, Yuan, Zhang, & Lee, 2013).

A striking feature of turbulence is the existence of multi-scale vortical structures with strong nonlinear interactions (Davidson, 2004; Pullin & Saffman, 1998). These structures, occupying a small fraction of volume, can have an impact on the momentum transport and generation of the turbulent kinetic energy, so they are described as the sinews of turbulent motion (Moffatt, Kida, & Ohkitani, 1994). To reduce the complexity on the mathematical treatment of the NS equations, we can take the curl of the velocity and then obtain the vorticity and its transport equation. The vorticity equation removes the non-locality due to the pressure term in the momentum equation, so the vorticity tends to concentrate in local regions and forms organized structures (Jimenez, Wray, Saffman, & Rogallo, 1993; She, Jackson, & Orszag, 1990). Based on the attached eddy hypothesis (Townsend, 1976), the statistically self-similar coherent structures inspire the structure-based models in wall turbulence (Chung & Pullin, 2009; Liu & Zheng, 2021; Marusic & Monty, 2019).

The global vorticity structures are integrated from the vorticity field, e.g. its integral line and surface are vortex line and surface, respectively. The vortex surface (also called vorticity surface Wu, Ma, & Zhou, 2015a) is tangent to the vorticity everywhere, and it generally has tube/sheet-like shapes. The importance of the vorticity is rooted in two fundamental theorems in fluid mechanics. First, the Helmholtz vorticity theorem (Helmholtz, 1858) and its extension (Hao, Xiong, & Yang, 2019) provide an elegant description of the motion of vortex surfaces. The theorem pointed out that the vortex surface is a material surface, as frozen in the velocity, in an inviscid flow. By contrast, the vortex surface can have topological changes (Kida & Takaoka, 1994; Yao & Hussain, 2022) in a viscous flow, so convoluted multi-scale vortical structures can be generated in flows at a large Reynolds number ($Re$). Second, the vorticity can be transformed into the velocity via the Biot–Savart theorem, and thus the global geometry and topology of vortex surfaces reflect evolutionary features of the velocity.

However, there is no consensus on the best method to identify vortical structures. There is even no clear mathematical definition for ‘vortex’. Most studies described the vortex as a bundle of vortex lines (Jimenez et al., 1993; She et al., 1990), or adopted the vortex criteria based on the Eulerian local velocity (Chong, Perry, & Cantwell, 1990; Hunt, Wray, & Moin, 1988; Jeong & Hussain, 1995; Liu, Gao, Tian, & Dong, 2018; Zhou, Adrian, Balachandar, & Kendall, 1999). The vortex criteria are capable of identifying the ‘vortex core’ – a region with strong local rotational motion. On the other hand, the evolution of such identified structures lacks a simple evolution equation, so one may only qualitatively describe the identified result and cannot develop a predictive model based on the identified structures. Additionally, there is no objective standard to choose the vortex criterion and the isocontouring threshold, so such vortex visualization is prone to subjectivity (Epps, 2017; Wu, 2018).

For comparison, the vortex surface in a Lagrangian framework contains the flow history of vortex dynamics, but its study was restricted to conceptional discussions due to the difficulty in the construction of vortex surfaces and the breakdown of the Helmholtz theorem in viscous flows. Tracking the vortex surface was believed to be impossible under the violation of the Helmholtz theorem, which is one of the major obstacles to understanding the timewise evolution of vortical structures in turbulence and transition.

The development of the VSF (Yang & Pullin, 2010) reshapes the structure identification and characterization based on the vorticity. The VSF is a three-dimensional (3-D) scalar field, whose isosurface is a vortex surface consisting of vortex lines. It is straightforward to extend the idea of VSF to identify other integral surfaces, e.g. the magnetic surface (Hao & Yang, 2021) and stream surface (Katsanoulis, Kogelbauer, Kaundinya, Ault, & Haller, 2023), in various vector fields in physics. The VSF method builds a quantitative framework to study the Lagrangian vortex dynamics, which extends the Helmholtz theorem to viscous flows and has been applied to various flow applications. Being distinguished from the Eulerian vortex criteria, which only consider the instantaneous local velocity, the VSF contains the history of vortex stretching and twisting in the view of Lagrangian-like tracking. This overcomes the weakness on the lack of time coherence and causality between the structures identified at different times in the existing Eulerian methods. In addition, the recently developed identification methods of flow structures also stressed the Lagrangian tracking (Haller, 2015; Zhu & Xi, 2019), non-locality (Wang & Peters, 2006) and objectivity (Haller, 2005).

The VSF serves not only as a post-processing flow diagnostic tool, but also as a framework to inspire novel flow simulation methods. The representation of Clebsch (1859) links the VSFs and the velocity–vorticity field. Pioneered by Kuznetsov and Mikhailov (1980), Chern, Knöppel, Pinkall, and Schröder (2017); Chern, Knöppel, Pinkall, Schröder, and Weißmann (2016) proposed the spherical Clebsch map that overcomes the difficulty of the classical Clebsch map to represent the velocity field with non-zero helicity, and meanwhile maintain the capability to reconstruct the VSF. Then, several numerical simulation methods based on the spherical Clebsch map have been developed in computer graphics. Chern et al. (2016) applied the Clebsch map to smoke simulation by solving the two-component Schrödinger equation (Madelung, 1927; Sorokin, 2001). Then, Yang et al. (2021) extended the Clebsch method with the gauge method (Saye, 2016) to simulate free-surface flows. Note that these methods for fast and stable computations have notable artificial viscosities, so they have not been applied to the computational fluid dynamics (CFD).

The outline of this paper is as follows. In § 2, we describe the theoretical framework and numerical methods for the VSF. In § 3, we show that the VSF visualizations facilitate the elucidation of critical mechanisms in turbulence, transition and other flow applications. In §§ 4 and 5, we review VSF-based models for predicting aerodynamic forces and mixing rates, and VSF-inspired methods for simulating vortical flows, respectively. Summary and outlook of the VSF study are given in § 6.

2. The VSF method

The VSF provides a Lagrangian-based framework for the identification, characterization and modelling of flow structures. From the theoretical perspective, the global geometry of the vorticity, as a family of integral surfaces of the vorticity, is represented by the VSF. By introducing the virtual circulation-preserving velocity (Hao et al., 2019), the theorems of Helmholtz (1858) and Ertel (1942) can be conceptually extended to viscous flows, so that the VSF evolution equation can be expressed in a Lagrangian-like conservation form. Thus the VSF isosurfaces of the same threshold at different times have strong coherence for tracking vortex surfaces.

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$ (Truesdell, 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 Ertel (1942), the evolution equation of the VSF can be written in a Lagrangian conserved form in an ideal flow (Hao 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 (He & 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 (Boratav & Pelz, 1994; Brachet et al., 1983).

Note that the classical potentials of Clebsch (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 (Nore, Abid, & Brachet, 1997; Xiong & 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 (Dombre 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 (Shen, Yao, Hussain, & Yang, 2023; Xiong & Yang, 2019a, 2020). The complexity of such fields can be characterized by the finite helicity (Moffatt, 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 (Katsanoulis et al., 2023; Yang & 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 (Yang & 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$ (Yang & 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) (Jiang & Shu, 1996) and targeted essentially non-oscillatory (Fu, 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 (Xiong & 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 (Xiong & 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) (Scarano, 2012). As shown in figure 4(a), Liu 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 (Ye, 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 (Zhao, 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 }$ (Xiong & 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 (Xiong & 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 (Hao et al., 2019; Peng & Yang, 2018; Yang & 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 (Kim & Moin, 1985) and the level-set method (Osher & 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 Yang and Pullin (2011a). Temporal integrations are marched by the second-order total-variation-diminishing Runge–Kutta scheme (Gottlieb & 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 (Jiang & 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 (Yang & 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 (Yang & 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 (Han & 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 (Weller, Tabor, Jasak, & Fureby, 1998). Since OpenFOAM lacks high-order, discontinuity-resolving schemes due to the challenges on stencil collection and memory management (Gä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) (Gu, 2022). Thus, the VSF calculation is applicable to general CFD softwares as a flow postprocessing module, which is illustrated in figure 4(b).

3. The VSF in flow visualization

Based on 3-D numerical/experimental data, the VSF visualization facilitates elucidating mechanisms in the flows with essential vortex dynamics, such as turbulence and transition. The universal framework of VSF evolution illustrates an ordered transformation of coherence structures in transitional flows, in which there is no need to change ad hoc the structure-identification method and isocontour threshold. Moreover, the quantitative VSF study sheds light on the interacting mechanisms between the vortex surface and the shock wave, flame or electromagnetic field in multi-physics coupled flows.

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 (Bermejo-Moreno, Pullin, & Horiuti, 2009; Cardesa, Vela-Martin, & Jimenez, 2017; Davidson, 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 Richardson (1922) and Kolmogorov (1941) speculated that the energy cascade from large to small scales is through a process of vortex breakup (Frisch, 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 (Lundgren, 1982; Pullin & 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 (She et al., 1990) and rolling-up of vortex sheets (Lindsay & Krasny, 2001).

Consistent with the vorticity dynamics, the morphology of VSF isosurfaces at different times in decaying TG and HIT flows (Han & Yang, 2022; Xiong & Yang, 2019b; Yang & 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’ (Mullin, 2011). Existing phenomenological models (Barkley, 2016; Barkley 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 (Zhao 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 (Zhao, 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 (Zhao 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 (Wu et al., 2017), consisting of multiple hairpins. The VSF was calculated from the direct numerical simulation (DNS) data (Sayadi, Hamman, & Moin, 2013) of a K-type boundary-layer transition using the boundary-constraint method (Xiong & Yang, 2017). Zhao, 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 (Ruan, 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 (Wu, 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 (Xiong, You, Ruan, & Yang, 2019). In contrast to the viscous cancellation of vortex reconnection (Kida & 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 (Wu et al., 2015b; Xiong 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 (Ji & van Rees, 2022; Melander & Hussain, 1994), which was found in aircraft trailing vortices (Tombach, 1973) and could be overlooked by Eulerian vortex criteria. Based on the VSF, Shen 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 (Shen 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 (Dabiri, 2005; Shyy et al., 2010; Zhang, Huang, Pan, Yang, & Huang, 2022) and aerial/underwater vehicles (Jiménez, Hultmark, & Smits, 2010; Park, 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 (Tong, Yang, & Wang, 2021) simulated by the immersed boundary method (Wang & 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 (Tong 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 (Kang, Liu, Su, & Wu, 2018; Wu 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 (Tong 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 (Vassberg, Dehaan, Rivers, & Wahls, 2008) in a wing–body configuration (Gu, 2022). From the simulation data with the Spalart–Allmaras model (Spalart & 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 (Yang & 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 (Brouillette, 2002) and turbulence (Samtaney, Pullin, & Kosovic, 2001; Wang 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 (Peng & 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 (Almgren, Bell, Rendleman, & Zingale, 2006) and ignition in inertial confinement fusion (ICF) (Aglitskiy et al., 2010). Peng, 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, Peng, 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 (Peng et al., 2021b) and (b) premixed flame (location marked by the red arrow) propagating in HIT (You 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 (Renard, Thévenin, Rolon, & Candel, 2000), but its study was restricted to 2-D configurations (Poinsot, Veynante, & Candel, 1991). The VSF study illustrated the influence of a propagating premixed flame on 3-D TG vortices (Zhou 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, You, 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 (Davidson, 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) (Hao & 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 (Linton, Dahlburg, & Antiochos, 2001) with opposite chirality. In particular, the MSF can pinpoint the region of the magnetic reconnection (Priest & 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. Hao 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 (Epps, 2017; Wu, 2018). In boundary-layer transition (Sayadi et al., 2013; Wu 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 Zhao 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 (Wu 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 (Peng & 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 (Zhao 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 (Zhao et al., 2018, 2016b).

4. The VSF in flow modelling

The idea of the VSF-based or VSF-inspired modelling is to extract the essential information of the velocity/vorticity field using the VSF or a characteristic scalar, and then explore robust statistics by characterizing scalar features. From the quantified flow mechanism, we model the statistical correlation between the statistical features and the critical physical quantity to be predicted in engineering applications.

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 (Zhao et al., 2018, 2016b). 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 (Lee & 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 (Liu & Shen, 2007; Su & Dahm, 1996).

The multi-scale geometric analysis was applied to tracking a passive scalar (Candes, Demanet, Donoho, & Ying, 2006; Yang, Pullin, & Bermejo-Moreno, 2010). From the band-pass filtering, the geometric features of the scalar were quantified at different scales. Yang 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$ (Zheng, Yang, & Chen, 2016). As shown in figure 14, a model for predicting $C_f$ (Zheng, Ruan, Yang, He, & Chen, 2019) was proposed from multi-scale structural geometry in 2-D images in supersonic boundary-layer transition (He, 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 (Dhawan & Narasimha, 1958) in turbulence modelling (Suzen & Huang, 2000), so it blends the empirical models of $C_f$ in laminar and fully developed turbulent states (White, 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 (Zheng 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 (White, 2006), respectively.

4.2 Aerodynamic forces

Shedding vortices in the far wake, e.g. in the wake of a flying bird (Hubel, Riskin, Swartz, & Breuer, 2010) or swimming fish (Lauder, 2015), can imply the forces exerted on the moving body (Li & Lu, 2012). Based on the VSF study and the vortex impulse theory (Wu, Liu, & Liu, 2018), Tong 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. Tong, 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) (He, 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 (Zhou, 2017a, 2017b), 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, Peng 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 (Peng 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 (Peng 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. The VSF in flow simulation

The VSF evolution captures essential Lagrangian-based dynamics and geometry of vortical flows. This advantage can be used as an important ingredient to develop novel numerical frameworks and methods in CFD. The velocity and multiple VSFs are linked by the Clebsch representation (Clebsch, 1859). Several fast algorithms based on the Clebsch representation instead of the conventional velocity or vorticity have been developed for cutting-edge hardware platforms, e.g. graphic and quantum processors.

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.

Kuznetsov and Mikhailov (1980) proposed the spherical Clebsch map (Chern et al., 2017, 2016) 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 (Bloch, 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$ (Hopf, 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 (Chern et al., 2016; Xiong, Wang, Wang, & Zhu, 2022; Yang et al., 2021). Tao, 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 (Kedia, 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 (Xiong 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}

Yang 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 (Stam, 1999), a staggered grid was adopted to discretize fluid variables and solve the equations using the advection–projection splitting method (Kim & 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 (Yang & Pullin, 2011b; Yang et al., 2010). To maintain the stability of the advection, Yang 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 (Stam, 1999) into the advection of $\boldsymbol \psi$, and the real viscosity is expected to be incorporated in future work.

Yang 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 (Yang 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 (Xiong et al., 2022; Yang et al., 2021). Assuming the fluid is subject to gravity and surface tension and according to the Young–Laplace equation (Batchelor, 1967), $\varPi$ in (5.4) satisfies a jump condition on $\varGamma$. The interface is tracked by a level-set function (Fedkiw & Osher, 2002; Osher & Sethian, 1988), and (5.6) with the jump condition is solved by the ghost fluid method (Fedkiw, 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 (Xiong 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 (Nielsen & Chuang, 2010). To date, quantum computing has been demonstrated to be effective in handling some linear problems (Clader, Jacobs, & Sprouse, 2013; Harrow, Hassidim, & Lloyd, 2009), but it remains challenging to efficiently solve the highly nonlinear NS equations (Liu et al., 2021; Lloyd et al., 2020).

The fluid dynamics of a potential flow was related to a hydrodynamic representation of quantum mechanics via the Madelung transform (Madelung, 1927). By generalizing the Madelung transform (Chern et al., 2016; Meng & 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 (Hao et al., 2019) and has helicity conservation (Meng, 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 (Koch, Wessing, & Alsing, 2019).

In particular, the incompressible HSE with unity density governs the incompressible Schrödinger flow (ISF) (Chern 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 (Meng & Yang, 2023; Tao 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 (Chern et al., 2017, 2016; Tao 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 (Madeira, Caracanhas, dos Santos, & Bagnato, 2020; Müller, Polanco, & Krstulovic, 2021) and the tangle of spiral vortex tubes and sheets in classical turbulence (Cardesa et al., 2017; She et al., 1990; Xiong & 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.

6. Summary and outlook

In the Lagrangian viewpoint, the VSF study paves an avenue to investigate the mechanism of energy cascade and evolution of coherent structures by characterizing vortex surfaces. Being different from the Eulerian vortex criteria based on the instantaneous local flow field, the VSF evolution contains the evolutionary history of vortex stretching and twisting, overcoming the weakness of lacking causality and time coherence of identified flow structures. As a general flow diagnostic tool, the numerical VSF solution can be constructed in any given flow field by solving a pseudo-transport equation driven by the instantaneous frozen vorticity. Then the evolution of VSFs is calculated by the two-time method.

In high-$Re$ flows, the VSF study reveals how a large-scale vortex surface evolves into small-scale structures through dynamic processes of stretching, rolling-up, twisting and reconnection. The large-scale structures have simple, ordered geometry, whereas the small-scale structures have a convoluted geometry manifesting the random nature. The evolution of VSF elucidates that the transition from the ordered to random structures is not an abrupt process but a timewise evolution that can be quantified. Thus, the VSF study embraces the two approaches for studying turbulence from statistical and structural perspectives.

The statistical features in VSF evolution are correlated to the key physical quantities in engineering applications. This link facilitates the development of VSF-based or VSF-inspired models, e.g. the models of the skin friction in supersonic boundary-layer transition based on the inclination angle of small-scale flow structures, the thrust of a flapping plate based on shedding vortex surfaces in the far wake and the mixing rate based on vortex rings with different vorticity-induced velocities in the RMI. In this way, the VSF provides a systematic approach to develop simple algebraic models for complex flows.

The velocity can be represented by VSFs via Clebsch maps. This Lagrangian-based representation involving the flow history inspires several fast algorithms on new-generation graphic and quantum computing hardware. Since vortical structures are well preserved using VSFs, the Clebsch method is capable of simulating a wide range of fluid flows that were impractical for previous methods in computer graphics, including complex vortex filaments, bubble rings and vortex–solid interactions. Furthermore, the HSE provides a possible way of simulating fluid dynamics on quantum computers (Succi, Itani, Sreenivasan, & Steijl, 2023). These progresses enlighten the development of revolutionary CFD methods based on VSFs.

Some future issues of the VSF study are listed below.

1. (i) Application of the VSF in flow diagnostics and simulation. The VSF is expected to combine with the high-resolution large-scale DNS database (Wu et al., 2015b, 2017), 3-D time-resolved near-wall PIV techniques (He et al., 2011; Wang, Gao, Wang, Wang, & Pan, 2019; Zhu et al., 2013), data-driven methods (Duraisamy, Iaccarino, & Xiao, 2019) and advanced GPU and quantum computing hardware. These interdisciplinary studies can extend the VSF to solving more critical problems in engineering sciences, such as the transition on high-speed vehicles, turbulent combustion in aero-engine and swimming/flying of bio-inspired robotics.

2. (ii) Geometric description of flows based on vortex surfaces. A universal Clebsch representation is to be developed to transform a 3-D velocity field into a two-component wave function. The exact expressions in spherical or other non-Euclidian spaces can be used to represent complex vortex surfaces. The constructed wave function serves as the initial condition of the generalized Schrödinger equation. Its evolution characterizes fluid motion and tracks vortex surfaces in turbulence and transition. In such simulations, the regularization method needs to be developed to treat viscous effects. In addition, reduced-order models and subgrid-scale models could be developed based on the relation between vortex geometry/topology and turbulence statistics.

3. (iii) Numerical issues in the VSF method. Although the VSF has been successfully applied to various problems, its computational cost and ease of use are still important issues. The additional cost of the VSF calculation in solving (2.6), a simple advection equation, needs multiple pseudo-time steps and appropriate boundary conditions for $\phi _v$. Thus, it is desired to improve the efficiency in the VSF calculation using GPU and to adapt the VSF solver to various complex boundary conditions and unstructured grids in open-source software.