2.1 Numerical methods

The fluid velocity field $\boldsymbol{u}(\boldsymbol{x},t)$ of an incompressible viscous flow is governed by the Navier–Stokes (NS) equations

(2.1) $$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u},\quad \\ \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\quad \end{array}\right.\end{eqnarray}$$

where $\boldsymbol{x}$ denotes spatial coordinates, $t$ the time, $p$ the pressure, $\unicode[STIX]{x1D70C}$ the density and $\unicode[STIX]{x1D708}$ the kinematic viscosity.

The DNS of decaying HIT is performed to solve (2.1) in a periodic box of side $L=2\unicode[STIX]{x03C0}$ using a standard pseudo-spectral method (see Rogallo Reference Rogallo1981; Yang et al. Reference Yang, Pullin and Bermejo-Moreno2010). The computational domain $\unicode[STIX]{x1D6FA}$ is discretized on uniform grid points $N^{3}$ . Aliasing errors are removed using the two-thirds truncation method with the maximum wavenumber $k_{max}\approx N/3$ . The Fourier coefficient of the velocity is advanced in time using a second-order Adams–Bashforth method, and the time step is chosen to ensure that the Courant–Friedrichs–Lewy (CFL) number is less than $0.5$ for numerical stability and accuracy.

The initial velocity $\boldsymbol{u}_{0}\equiv \boldsymbol{u}(\boldsymbol{x},t=0)$ for decaying HIT is a random Gaussian field with the initial total energy $E_{0}\equiv \int |\boldsymbol{u}_{0}|^{2}/2\,\text{d}\boldsymbol{x}=1$ and a prescribed energy spectrum (e.g. Kraichnan Reference Kraichnan1970; Ishida, Davidson & Kaneda Reference Ishida, Davidson and Kaneda2006)

(2.2) $$\begin{eqnarray}E(k,t=0)\sim k^{4}\text{e}^{-2(k/k_{p})^{2}},\end{eqnarray}$$

where $k=|\boldsymbol{k}|$ denotes the wavenumber magnitude. The wavenumber $k_{p}=4$ where $E(k,t=0)$ peaks was selected to ensure that the fully developed HIT field has a broad range of length scales, and it is comparable to the choices of $k_{p}$ in previous DNS (e.g. Hosokawa & Yamamoto Reference Hosokawa and Yamamoto1986; Elghobashi & Truesdell Reference Elghobashi and Truesdell1992; Mansour & Wray Reference Mansour and Wray1994).

Additionally, in order to illustrate the identification of twisted vortex tubes, we carry out another DNS of TG flows (Taylor & Green Reference Taylor and Green1937; Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983) with the initial velocity

(2.3) $$\begin{eqnarray}\boldsymbol{u}_{0}=(\sin x\cos y\cos z,-\cos x\sin y\cos z,0)\end{eqnarray}$$

at various Reynolds numbers $Re=1/\unicode[STIX]{x1D708}$ . The initial set-up parameters of the DNS of decaying HIT, as well as the TG flows at four $Re$ , are summarized in table 1.

Table 1.

Set-up of DNS cases.

2.2 DNS statistics

Important statistics in the DNS of decaying HIT at six typical times are summarized in table 2, including the total turbulent kinetic energy $E_{tot}=\int E(k)\,\text{d}k$ , mean dissipation rate $\unicode[STIX]{x1D700}=2\unicode[STIX]{x1D708}\int k^{2}E(k)\,\text{d}k$ , root-mean-square velocity $u^{\prime }=(2E_{tot}/3)^{1/2}$ , Taylor–Reynolds number $Re_{\unicode[STIX]{x1D706}}=u^{\prime }\unicode[STIX]{x1D706}_{T}/\unicode[STIX]{x1D708}$ with the Taylor micro-length scale $\unicode[STIX]{x1D706}_{T}=(15\unicode[STIX]{x1D708}{u^{\prime }}^{2}/\unicode[STIX]{x1D700})^{1/2}$ , Kolmogorov length scale $\unicode[STIX]{x1D702}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D700})^{1/4}$ , Kolmogorov time scale $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D700})^{1/2}$ , spatial resolution $k_{max}\unicode[STIX]{x1D702}$ and integral length scale $L_{e}=\unicode[STIX]{x03C0}/(2{u^{\prime }}^{2})\int E(k)/k\,\text{d}k$ . We remark that the spatial resolution in the present DNS is always greater than two, which satisfies the common criterion $k_{max}\unicode[STIX]{x1D702}>1.5$ (Pope Reference Pope2000) for resolving the smallest scales in HIT.

Table 2.

Statistics in the DNS of decaying HIT at six typical times.

From the temporal evolution of $E_{tot}$ , $\unicode[STIX]{x1D700}$ and $E(k)$ in the DNS of decaying HIT in figures 1(a), 1(b) and 2(a), we generally divide the flow evolution into two stages. In the early stage, the dissipation increases and the flow field evolves from an artificially constructed random field to the fully developed turbulent state. In the later stage, the dissipation decreases and the turbulent flow gradually decays towards a stationary flow. In figure 1(b), the dissipation peaks around $t\approx 2$ , implying that the flow field has the maximum total enstrophy with the most intensive vortical structures. In figure 1(a), the total energy decays with the decay law $E_{tot}\sim t^{-10/7}$ (Kolmogorov Reference Kolmogorov1941a ). In figure 2(a), the energy spectrum is gradually broadened in the early stage and has a short inertial region with the five-thirds law $E(k)\sim k^{-5/3}$ and a range of length scales in the HIT field at $t=2$ . Subsequently the entire spectrum decays in the late stage after $t=10$ .

The evolution of TG flows also has two stages similar to the case of decaying HIT (see Brachet et al. Reference Brachet, Meiron, Orszag, Nickel, Morf and Frisch1983; Yang & Pullin Reference Yang and Pullin2011). Additionally, as $Re$ increases, the energy spectrum in TG flows at the two-stage separation point $t=7.5$ in figure 2(b) is broadened and approaches the fully developed turbulent state with the five-thirds law.

Figure 1.

Evolution of (a) total energy and (b) mean dissipation rate in HIT.

Figure 2.

Energy spectra in (a) HIT at different times and (b) TG flows at $t=7.5$ .

3.1 VSF constraint and pseudo-transport equation

The VSF $\unicode[STIX]{x1D719}_{v}$ is defined as a globally smooth scalar field whose isosurface is a vortex surface consisting of vortex lines. The constraint of the VSF (Yang & Pullin Reference Yang and Pullin2010) is

(3.1) $$\begin{eqnarray}{\mathcal{C}}_{v}\equiv \unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}=0,\end{eqnarray}$$

where $\unicode[STIX]{x1D74E}\equiv \unicode[STIX]{x1D735}\times \boldsymbol{u}$ is the vorticity. The construction of VSF is equivalent to finding a solution of $\unicode[STIX]{x1D719}_{v}$ satisfying (3.1) for a given, instantaneous $\unicode[STIX]{x1D74E}$ . For a vorticity field with vanishing helicity density $h\equiv \boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D74E}=0$ , we can obtain exact VSF solutions from (3.1) (He & Yang Reference He and Yang2016). Nonetheless, for a general vorticity field with $h\neq 0$ , the non-trivial analytic solution of (3.1) may not exist, so we have to seek an approximate solution of (3.1) using numerical methods.

Based on the definition of VSF, the deviation of an approximate VSF solution $\unicode[STIX]{x1D719}_{v}$ from an exact VSF is defined as the cosine of the angle between $\unicode[STIX]{x1D74E}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}$ as (Yang & Pullin Reference Yang and Pullin2010)

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}\equiv \frac{\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}}{|\unicode[STIX]{x1D74E}||\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|}.\end{eqnarray}$$

For a given instantaneous vorticity field $\unicode[STIX]{x1D74E}$ , the goal of constructing an approximate VSF solution is to minimize the volume-averaged VSF deviation $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ by varying $\unicode[STIX]{x1D719}_{v}$ .

To obtain an approximate VSF solution, Yang & Pullin (Reference Yang and Pullin2011) developed the pseudo-transport equation

(3.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{v}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}+\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}=0,\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},\quad 0<\unicode[STIX]{x1D70F}\leqslant T_{\unicode[STIX]{x1D70F}}\end{eqnarray}$$

driven by the given ‘frozen’ vorticity, where $T_{\unicode[STIX]{x1D70F}}$ denotes the largest pseudo-time in the VSF calculation. The initial condition $\unicode[STIX]{x1D719}_{v0}\equiv \unicode[STIX]{x1D719}_{v}(\boldsymbol{x},\unicode[STIX]{x1D70F}=0)$ of (3.3) can be an arbitrary scalar field with large $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ . Although the exact VSF solution is not unique (see Yang & Pullin Reference Yang and Pullin2010), we found that the statistical geometry of converged VSF solutions is insensitive to the choice of $\unicode[STIX]{x1D719}_{v0}$ if the given vorticity field is chaotic (see Xiong & Yang Reference Xiong and Yang2017), e.g. in turbulence. The sensitivity of the VSF solution to initial conditions is further discussed in appendix A. Moreover, the boundary conditions for (3.3) should be treated carefully except for the periodic boundary condition that can be naturally applied (see Xiong & Yang Reference Xiong and Yang2017).

3.2 Numerical construction of VSFs with the local optimization

Although the VSF deviation of the solution of (3.3) can be converged at large pseudo-time $T_{\unicode[STIX]{x1D70F}}$ , the converged solution tends to have nearly singular structures with a very large scalar gradient owing to the persistent straining of $\unicode[STIX]{x1D719}_{v}$ driven by the frozen vorticity in turbulence or highly chaotic flows (see Pullin & Yang Reference Pullin and Yang2014). This dilemma hinders the effective VSF construction and visualization in complex flow fields.

Thus we propose a new hybrid constraint

(3.4) $$\begin{eqnarray}{\mathcal{C}}_{h}=(1-\unicode[STIX]{x1D701})|\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}+\unicode[STIX]{x1D701}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\end{eqnarray}$$

to balance the small VSF deviation with the desired smoothness of VSF solutions, where $\unicode[STIX]{x1D701}$ is a weighting factor to characterize the preference for the minimization of $|{\mathcal{C}}_{v}|$ or $|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|$ . The construction of a non-trivial normalized VSF solution is then equivalent to an optimization problem as

(3.5a,b ) $$\begin{eqnarray}\min {\mathcal{C}}_{h},\quad \text{s.t.}\langle \unicode[STIX]{x1D719}_{v}\rangle =0\quad \text{and}\quad \text{Var}(\unicode[STIX]{x1D719}_{v})=1,\end{eqnarray}$$

where $\text{Var}(f)\equiv \langle (f-\langle f\rangle )^{2}\rangle$ denotes the variance of an arbitrary function $f$ in a volume.

With this key improvement, each computational step in the implementation of VSF construction is divided into two sub-steps: (i) solving the pseudo-transport equation (3.3); (ii) applying the local optimization (3.5). First, a temporary VSF $\unicode[STIX]{x1D719}_{v}^{\ast }$ is evolved as (3.3) in pseudo-time, where $\unicode[STIX]{x1D719}_{v}^{\ast }$ is advanced in $\unicode[STIX]{x1D70F}$ using the third-order TVD Runge–Kutta method, and the convection term is approximated by the fifth-order weighted essentially non-oscillatory (WENO) scheme (Jiang & Shu Reference Jiang and Shu1996). The pseudo-time step should be small enough to ensure that the CFL number based on vorticity is less than 0.5 for numerical stability.

Subsequently, we minimize ${\mathcal{C}}_{h}$ to obtain a local optimized solution $\tilde{\unicode[STIX]{x1D719}}_{v}$ , where the weighting factor in (3.4) is chosen as

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D70F})=(1-\unicode[STIX]{x1D70F}/T_{\unicode[STIX]{x1D70F}})^{2}.\end{eqnarray}$$

In other words, we prefer smooth VSF solutions for small $\unicode[STIX]{x1D70F}$ and accurate solutions for large $\unicode[STIX]{x1D70F}$ . It is noted that the specific form of (3.6) is ad hoc based on both the accuracy and smoothness of VSF solutions in numerical experiments. In general, if $\unicode[STIX]{x1D701}=1$ or $\unicode[STIX]{x1D701}$ is closer to unity than (3.6), $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ and $\sqrt{\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle }$ tend to slightly increase and decrease at large $\unicode[STIX]{x1D70F}$ , respectively, but geometries of corresponding VSF isosurfaces are almost identical. Furthermore, appendix B provides a two-dimensional example to illustrate the effect of the local optimization and the choice of $\unicode[STIX]{x1D701}$ on VSF calculation.

The numerical algorithm of minimizing the hybrid constraint ${\mathcal{C}}_{h}$ using the variational method is described in detail below. The minimization of ${\mathcal{C}}_{h}$ in (3.4) by varying $\unicode[STIX]{x1D719}_{v}^{\ast }$ is equivalent to a local optimization problem satisfying

(3.7a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\Vert {\mathcal{C}}_{h}(\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}),\boldsymbol{x})\Vert _{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a})}=0\quad \text{and}\quad \frac{\unicode[STIX]{x2202}^{2}\Vert {\mathcal{C}}_{h}(\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}),\boldsymbol{x})\Vert _{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a})^{2}}>0.\end{eqnarray}$$

Here, $\boldsymbol{a}$ denotes the coordinates at which $\unicode[STIX]{x1D719}_{v}^{\ast }$ is varied in the optimization, and $\Vert f\Vert _{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}}\equiv \int _{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}}|f|\,\text{d}\boldsymbol{x}$ denotes the $L_{1}$ norm of a function $f$ over a finite subdomain $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}$ containing $\boldsymbol{a}$ .

The temporary scalar field $\unicode[STIX]{x1D719}_{v}^{\ast }$ is discretized on uniform grid points $N_{\unicode[STIX]{x1D719}}^{3}$ in computational domain $\unicode[STIX]{x1D6FA}$ . Coordinates for the uniform grids are $\boldsymbol{x}_{\boldsymbol{i}}\equiv (x_{i},y_{j},z_{k})=(i\unicode[STIX]{x0394}x,j\unicode[STIX]{x0394}x,k\unicode[STIX]{x0394}x)$ with grid spacing $\unicode[STIX]{x0394}x=2\unicode[STIX]{x03C0}/N_{\unicode[STIX]{x1D719}}$ , which can be expressed as a point set

(3.8) $$\begin{eqnarray}G=\{\boldsymbol{x}_{\boldsymbol{i}}\mid i,j,k\in \{0,1,\ldots ,N_{\unicode[STIX]{x1D719}}-1\}\}.\end{eqnarray}$$

In the local optimization, $\unicode[STIX]{x1D719}_{v}^{\ast }$ is varied on coarse staggered grids. Coordinates for the staggered grids are $\boldsymbol{a}_{\boldsymbol{m}}\equiv (a_{m},b_{n},c_{l})=(m\unicode[STIX]{x0394}x,n\unicode[STIX]{x0394}x,l\unicode[STIX]{x0394}x)$ . The $N_{\unicode[STIX]{x1D719}}^{3}/(n_{s}+1)$ variational points are selected as a subset of $G$ as

(3.9) $$\begin{eqnarray}G^{\prime }=\{\boldsymbol{a}_{\boldsymbol{ m}}\mid m,n,l\in \{0,1,\ldots ,N_{\unicode[STIX]{x1D719}}-1\}\text{ and }m+n+l\equiv 0~(\text{mod}~n_{s}+1)\},\end{eqnarray}$$

where $n_{s}$ is the order of the finite difference scheme utilized for evaluating ${\mathcal{C}}_{h}$ , and satisfies $N_{\unicode[STIX]{x1D719}}\equiv 0~(\text{mod}~n_{s}+1)$ . The staggered grids and a specific finite difference scheme define staggered overlapping subdomains $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}$ , so that the variation of $\unicode[STIX]{x1D719}_{v}^{\ast }$ at $\boldsymbol{a}_{\boldsymbol{m}}$ only influences ${\mathcal{C}}_{h}$ in each subdomain. In the schematic diagram in figure 3, each $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}$ is enclosed by red dashed lines, and it only has one variational point at $\boldsymbol{a}_{\boldsymbol{m}}$ marked by the open circle to avoid the interference among the variations of $\unicode[STIX]{x1D719}_{v}^{\ast }$ at multiple points in the subdomain.

The discretized $L_{1}$ norm of ${\mathcal{C}}_{h}(\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{m}}),\boldsymbol{x}_{\boldsymbol{i}})$ over each $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}$ is minimized by varying $\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{m}})$ . A third-order difference scheme with $n_{s}=3$ is applied to compute the gradient of $\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{x}_{\boldsymbol{i}})$ in ${\mathcal{C}}_{h}(\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{m}}),\boldsymbol{x}_{\boldsymbol{i}})$ , e.g. the partial derivative respect to $x$ at specific $\boldsymbol{x}_{\boldsymbol{i}}$ on the left of $\boldsymbol{a}_{\boldsymbol{m}}=\boldsymbol{x}_{\boldsymbol{i}}+\unicode[STIX]{x0394}\boldsymbol{x}$ (see figure 3) with $\unicode[STIX]{x0394}\boldsymbol{x}=(\unicode[STIX]{x0394}x,0,0)$ is discretized as

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{x}_{\boldsymbol{i}})}{\unicode[STIX]{x2202}x}=\frac{1}{\unicode[STIX]{x0394}x}\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{ m}})-\frac{1}{6\unicode[STIX]{x0394}x}[2\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{x}_{\boldsymbol{ i}}-\unicode[STIX]{x0394}\boldsymbol{x})+3\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{x}_{\boldsymbol{ i}})+\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{x}_{\boldsymbol{ i}}+2\unicode[STIX]{x0394}\boldsymbol{x})].\end{eqnarray}$$

The partial derivatives respect to $y$ and $z$ have a similar form to (3.10). Substituting all the discretized partial derivatives in three directions into (3.4) yields the discretized hybrid constraint in a polynomial form as

(3.11) $$\begin{eqnarray}{\mathcal{C}}_{h}(\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{ m}}),\boldsymbol{x}_{\boldsymbol{i}})=\unicode[STIX]{x1D709}_{2}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}})\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{ m}})^{2}+\unicode[STIX]{x1D709}_{1}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}})\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{ m}})+\unicode[STIX]{x1D709}_{0}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}}),\end{eqnarray}$$

where $\unicode[STIX]{x1D709}_{0}$ , $\unicode[STIX]{x1D709}_{1}$ and $\unicode[STIX]{x1D709}_{2}$ are obtained by collecting finite difference coefficients in the discretized partial derivatives such as (3.10). After the discretization, (3.7) is re-expressed as

(3.12) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{m}})}\left[\mathop{\sum }_{\boldsymbol{x}_{\boldsymbol{i}}\in \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}}{\mathcal{C}}_{h}(\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{ m}}),\boldsymbol{x}_{\boldsymbol{i}})\right]=\mathop{\sum }_{\boldsymbol{x}_{\boldsymbol{i}}\in \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}}[2\unicode[STIX]{x1D709}_{2}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}})\unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{a}_{\boldsymbol{ m}})+\unicode[STIX]{x1D709}_{1}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}})]=0,\end{eqnarray}$$

and the corresponding second partial derivative only containing squared terms must be positive. Thus scalar values at all the variational points $\boldsymbol{a}_{\boldsymbol{i}}$ are optimized as

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{v}^{\prime }(\boldsymbol{a}_{\boldsymbol{ m}})=-\frac{1}{2}\mathop{\sum }_{\boldsymbol{x}_{\boldsymbol{i}}\in \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}}\unicode[STIX]{x1D709}_{1}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}})\bigg/\mathop{\sum }_{\boldsymbol{x}_{\boldsymbol{i}}\in \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D6FA}_{\boldsymbol{m}}}\unicode[STIX]{x1D709}_{2}(\boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{a}_{\boldsymbol{m}})\end{eqnarray}$$

to minimize the local hybrid constraint.

Figure 3.

Schematic diagram for the subdivision of the computational domain.

Finally, the VSF solution in the local optimization sub-step is updated as

(3.14) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D719}}_{v}(\boldsymbol{x}_{\boldsymbol{i}})=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D719}_{v}^{\prime }(\boldsymbol{x}_{\boldsymbol{i}}),\quad & \text{if }\boldsymbol{x}_{\boldsymbol{i}}\in G^{\prime },\\ \unicode[STIX]{x1D719}_{v}^{\ast }(\boldsymbol{x}_{\boldsymbol{i}}),\quad & \text{if }\boldsymbol{x}_{\boldsymbol{i}}\in G\setminus G^{\prime }.\end{array}\right.\end{eqnarray}$$

Then the VSF solution is normalized to zero mean and unity variance as

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{v}=\frac{\tilde{\unicode[STIX]{x1D719}}_{v}-\langle \tilde{\unicode[STIX]{x1D719}}_{v}\rangle }{\sqrt{\text{Var}(\tilde{\unicode[STIX]{x1D719}}_{v})}}.\end{eqnarray}$$

3.3 Assessment of VSF construction

We construct VSFs from instantaneous vorticity fields in DNS of decaying HIT at six typical times listed in table 2. The number of grid points $N_{\unicode[STIX]{x1D719}}^{3}$ for the VSF, as listed in table 3, can be larger than required $N^{3}$ for DNS, and $\unicode[STIX]{x1D74E}$ from DNS is interpolated onto the finer VSF grids by the cubic spline interpolation for further VSF construction. In general, required $N_{\unicode[STIX]{x1D719}}^{3}$ increases with decreasing the tolerance of $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ or increasing the complexity of vorticity.

Table 3.

The number of scalar grid points and CPU hours in the construction of VSFs in decaying HIT at six typical times.

The initial condition $\unicode[STIX]{x1D719}_{v0}$ is a spatial delta-correlated field in which the scalar value at each point is generated independently as a random number satisfying the uniform distribution from 0 to 1. Although a tailored smooth $\unicode[STIX]{x1D719}_{v0}$ should achieve faster convergence in VSF construction, the present choice avoids the subjective selection of $\unicode[STIX]{x1D719}_{v0}$ .

Figure 4 displays the plane cuts of VSF solution contours at three pseudo-times during VSF construction for the HIT field at $t=20$ . In the pseudo-evolution, the scalar field evolves from a very noisy initial random field through a fine-scale scalar field at $\unicode[STIX]{x1D70F}/T_{v}=1$ to a relatively smooth field at $\unicode[STIX]{x1D70F}/T_{v}=4$ when the VSF solution has been nearly converged. Here, $T_{v}\equiv L_{\unicode[STIX]{x1D6FA}}/\langle |\unicode[STIX]{x1D74E}|\rangle$ is the characteristic pseudo-time scale (Xiong & Yang Reference Xiong and Yang2017), where the characteristic length scale of the computational domain $L_{\unicode[STIX]{x1D6FA}}=2\sqrt{3}\unicode[STIX]{x03C0}$ is the diagonal length of the DNS cube.

Figure 4.

Plane cuts of the contour of VSF solutions at three pseudo-times during the VSF construction for HIT at $t=20$ . (a)  $\unicode[STIX]{x1D70F}/T_{v}=0$ ; (b)  $\unicode[STIX]{x1D70F}/T_{v}=1$ ; (c)  $\unicode[STIX]{x1D70F}/T_{v}=4$ .

We compare the convergence and smoothness of VSF solutions with and without local optimization. In the VSF construction for instantaneous $\unicode[STIX]{x1D74E}$ of decaying HIT at $t=2$ , 20 and 60 with $N_{\unicode[STIX]{x1D719}}=512$ , pseudo-evolutions of $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ and $\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle$ in normalized pseudo-time $\unicode[STIX]{x1D70F}/T_{v}$ are shown in figures 5(a) and 5(b), respectively. In general, $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ and $\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle$ of VSF solutions have converged before $\unicode[STIX]{x1D70F}/T_{v}=10$ , so we set $T_{\unicode[STIX]{x1D70F}}=10T_{v}$ . Additionally, $\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle$ without local optimization becomes larger at higher $Re_{\unicode[STIX]{x1D706}}$ . By contrast, the pseudo-transport with local optimization significantly reduces the scalar gradient to ensure both satisfactory VSF deviation $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle \leqslant 10\,\%$ and smoothness of VSF solutions, although the optimization strategy of $\unicode[STIX]{x1D701}(\unicode[STIX]{x1D70F})$ in (3.6) slightly sacrifices the convergence rate of VSF solutions for small pseudo-times.

Figure 5.

The pseudo-evolution of (a)  $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ and (b)  $\sqrt{\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle }$ in decaying HIT at $t=2$ , 20 and 60 with $N_{\unicode[STIX]{x1D719}}=512$ and with or without local optimization.

The averaged VSF deviation can be further reduced by increasing the grid resolution for VSF (see Yang & Pullin Reference Yang and Pullin2011). Taking the VSF construction for the most fluctuating HIT field at $t=2$ for example, $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ is reduced from $18\,\%$ with $N_{\unicode[STIX]{x1D719}}^{3}=256^{3}$ to $6\,\%$ with $N_{\unicode[STIX]{x1D719}}^{3}=1024^{3}$ in figure 6(a), and $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ for converged numerical VSF solutions is proportional to $N_{\unicode[STIX]{x1D719}}^{-0.8}$ . At the meantime, $\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle$ is increased from $10$ with $N_{\unicode[STIX]{x1D719}}^{3}=256^{3}$ to $28$ with $N_{\unicode[STIX]{x1D719}}^{3}=1024^{3}$ in figure 6(b), indicating that the high-resolution VSF solution reveals finer-scale structures of vortex surfaces.

Figure 6.

The pseudo-evolution of (a)  $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ and (b)  $\sqrt{\langle |\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{v}|^{2}\rangle }$ in decaying HIT at $t=2$ on different VSF grids with local optimization.

We remark that the VSF construction is computationally expansive, which needs to be improved in future work. The computational time depends on the number of VSF grid points, the tolerance of $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ and the complexity of vorticity (see Xiong & Yang Reference Xiong and Yang2017). Typical CPU hours for constructing VSFs in decaying HIT are summarized in table 3.

5.1 Isocontour values and VSF deviations

As demonstrated in §§ 3.3 and 4, we are able to construct relatively accurate and smooth VSFs in HIT by improving the VSF calculation with local optimization (3.5), and the VSF isosurface can effectively identify coherent and twisted vortex tubes. These indicate that the VSF is promising to identify the very convoluted vortex surfaces in fully developed HIT.

Isosurfaces of the VSF solution and vorticity magnitude in decaying HIT at six typical times listed in table 3 are shown in figures 9 and 10. The isosurfaces are colour coded by $|\unicode[STIX]{x1D74E}|$ , and the VSF isosurfaces with $|\unicode[STIX]{x1D74E}|<0.2|\unicode[STIX]{x1D74E}|_{max}$ are cut off for clarity, where $|\unicode[STIX]{x1D74E}|_{max}$ denotes the maximum vorticity magnitude in each instantaneous $|\unicode[STIX]{x1D74E}|$ .

Figure 9.

Isosurfaces of $\unicode[STIX]{x1D719}_{v}$ and $|\unicode[STIX]{x1D74E}|$ in decaying HIT at $t=0$ , 1 and 2. The isocontour values and averaged VSF deviations of $\unicode[STIX]{x1D719}_{v}$ and $|\unicode[STIX]{x1D74E}|$ are listed in table 4. All the isosurfaces are colour coded by $|\unicode[STIX]{x1D74E}|$ , and the VSF isosurfaces with $|\unicode[STIX]{x1D74E}|<0.2|\unicode[STIX]{x1D74E}|_{max}$ are cut off for clarity. (a)  $\unicode[STIX]{x1D719}_{v}$ and (b)  $|\unicode[STIX]{x1D74E}|$ at $t=0$ ; (c)  $\unicode[STIX]{x1D719}_{v}$ and (d)  $|\unicode[STIX]{x1D74E}|$ at $t=1$ ; (e)  $\unicode[STIX]{x1D719}_{v}$ and (f)  $|\unicode[STIX]{x1D74E}|$ at $t=2$ .

Figure 10.

Same as figure 9 except $t=10$ , 20 and 60. Some vortex lines are integrated on the isosurfaces at $t=20$ . (a)  $\unicode[STIX]{x1D719}_{v}$ and (b)  $|\unicode[STIX]{x1D74E}|$ at $t=10$ ; (c)  $\unicode[STIX]{x1D719}_{v}$ and (d)  $|\unicode[STIX]{x1D74E}|$ at $t=20$ ; and (e)  $\unicode[STIX]{x1D719}_{v}$ and (f)  $|\unicode[STIX]{x1D74E}|$ at $t=60$ .

The probability density function (p.d.f.) of the space-filling, normalized $\unicode[STIX]{x1D719}_{v}$ evolves from the initial double-delta distribution towards the standard Gaussian distribution in the pseudo-transport with numerical diffusion and local optimization, and it has a strong peak at zero with a slight skewness in figure 11(a). Thus we use the isocontour level $\unicode[STIX]{x1D719}_{v}=0$ to display the typical VSF isosurface which almost occupies the largest space of all the VSF isosurfaces. We remark that the subjective selection of isocontour levels is an important issue suffered by many vortex-identification methods. As listed in table 4, we fix the isocontour level $\unicode[STIX]{x1D719}_{v}=0$ in the visualization of VSFs, whereas we have to tune the threshold of normalized $|\unicode[STIX]{x1D74E}|$ to display the structures with similar characteristic sizes in figures 9 and 10.

Figure 11.

Statistics of VSF solutions in decaying HIT at $t=2$ , 20 and 60. (a) Probability density function of $\unicode[STIX]{x1D719}_{v}$ ; (b) volume average of $|\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|$ conditioned on $|\unicode[STIX]{x1D74E}|$ .

Table 4.

Comparisons between the VSF and vorticity magnitude at the isocontour value selection and averaged VSF deviation in the vortex identifications in decaying HIT at six typical times.

On the faithful identification of vortex tubes and sheets, the VSF, due to its own definition, is superior to other vortex-identification methods. Table 4 compares the averaged VSF deviations between the VSF solution and vorticity magnitude in decaying HIT at six typical times, and $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ for $\unicode[STIX]{x1D719}_{v}$ is approximately five times less than that for $|\unicode[STIX]{x1D74E}|$ . As an example, some vortex lines are integrated on the isosurfaces at $t=20$ in figure 10(c,d). The large VSF deviation $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle \approx 30.8\,\%$ for the isosurface of $|\unicode[STIX]{x1D74E}|$ implies a significant misalignment between the isosurface and vortex lines in figure 10(d). The isosurface of $|\unicode[STIX]{x1D74E}|$ only displays strong vorticity regions, but loses directional information of vorticity, so it is hard to use in Lagrangian-like structure tracking based on the Helmholtz vorticity theorem owing to the large deviation from the vortex surfaces (see Yang & Pullin Reference Yang and Pullin2011). For comparison, the isosurface of $\unicode[STIX]{x1D719}_{v}$ displays the vortex tubes or sheets with a relatively small VSF deviation $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle \approx 5.5\,\%$ , so vortex lines are almost attached to the isosurface in figure 9(c), and the VSF isosurface colour coded by $|\unicode[STIX]{x1D74E}|$ captures both the magnitude and directional information of vorticity. Additionally, the geometry of the isosurfaces of other vortex-identification criteria based on the velocity gradient tensor is very similar to the isosurface of $|\unicode[STIX]{x1D74E}|$ .

The averaged VSF deviation conditioned on $|\unicode[STIX]{x1D74E}|$ in figure 11(b) implies that the construction of VSFs tends to be more accurate in the region with strong vorticity than in that with weak vorticity, because the visualization of vortex lines (She et al. Reference She, Jackson and Orszag1990; Jimenez et al. Reference Jimenez, Wray, Saffman and Rogallo1993) depicts that, in general, the vortex lines with large $|\unicode[STIX]{x1D74E}|$ are fairly straight while the vorticity vectors with small $|\unicode[STIX]{x1D74E}|$ have chaotic orientations.

5.2 Geometry of vortex surfaces in HIT

There is still considerable debate on the geometry of vortical structures in HIT (see Davidson Reference Davidson2004). Multi-scale diagnostic tools illustrate that structures of $|\unicode[STIX]{x1D74E}|$ evolve in a continuous distribution from blob-like and moderately stretched tube-like shapes at large scales, through a predominance of tube-like structures at intermediate scales, to highly stretched sheet-like structures at small scales (Ishihara et al. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007; Bermejo-Moreno, Pullin & Horiuti Reference Bermejo-Moreno, Pullin and Horiuti2009; Leung, Swaminathan & Davidson Reference Leung, Swaminathan and Davidson2012).

In figures 9 and 10, the VSF isosurfaces with strong and weak $|\unicode[STIX]{x1D74E}|$ , in general, tend to be tube-like and sheet-like, respectively. This observation is qualitatively consistent with that from isosurfaces of $|\unicode[STIX]{x1D74E}|$ . On the other hand, the VSF reveals more details of twisting and spirals of vortical structures than $|\unicode[STIX]{x1D74E}|$ . Plane cuts of the contour of $\unicode[STIX]{x1D719}_{v}$ in figure 12 depict the spiral structures in the contour of $\unicode[STIX]{x1D719}_{v}$ which are similar to passive scalars in HIT in Brethouwer, Hunt & Nieuwstadt (Reference Brethouwer, Hunt and Nieuwstadt2003) and Yang et al. (Reference Yang, Pullin and Bermejo-Moreno2010).

Figure 12.

Plane cuts of VSF and vorticity magnitude contours in decaying HIT. (a)  $\unicode[STIX]{x1D719}_{v}$ and (b)  $|\unicode[STIX]{x1D74E}|$ at $t=2$ ; (c)  $\unicode[STIX]{x1D719}_{v}$ and (d)  $|\unicode[STIX]{x1D74E}|$ at $t=20$ ; (e)  $\unicode[STIX]{x1D719}_{v}$ and (f)  $|\unicode[STIX]{x1D74E}|$ at $t=60$ .

The theoretical analysis in Yang et al. (Reference Yang, Pullin and Bermejo-Moreno2010) suggests that a passive scalar field $\unicode[STIX]{x1D719}$ with vanishing diffusivity in stationary HIT tends to be attracted towards a VSF, but $\langle |\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D714}}|\rangle$ is converged around 25 % owing to small but persistent viscous effects. The multi-scale geometric analysis of $\unicode[STIX]{x1D719}$ characterizes a geometrical progression from blobs through tubes to sheet-like structures with decreasing length scale in a space of reduced geometrical parameters. Compared with the statistical geometry of instantaneous $|\unicode[STIX]{x1D74E}|$ in HIT obtained by Bermejo-Moreno et al. (Reference Bermejo-Moreno, Pullin and Horiuti2009), $\unicode[STIX]{x1D719}$ tends to exhibit more prevalent sheet-like shapes at intermediate and small scales, which agrees with the visual difference between the VSF and $|\unicode[STIX]{x1D74E}|$ in figures 9, 10 and 12.

The generic spiral patterns of $\unicode[STIX]{x1D719}_{v}$ driven by frozen vorticity and the passive scalar driven by evolving velocity imply that the bundle of vortex lines can be highly twisted at large $Re_{\unicode[STIX]{x1D706}}$ , which is hypothesized as an elementary vortical structure in turbulence and is a key ingredient in the spiral vortex model for obtaining $E(k)\sim k^{-5/3}$ (see Lundgren Reference Lundgren1982; Pullin & Saffman Reference Pullin and Saffman1998). For comparison, the contours of large $|\unicode[STIX]{x1D74E}|$ in figure 12 are concentrated in blob-like regions, which correspond to the three-dimensional tube-like structures in figures 9 and 10 and appear not to illustrate the fine, internal structures in vortex tubes.

5.3 Implications of the morphology of vortex surfaces for energy cascade

The energy cascade is hypothesized to be related to the morphology of vortical structures at different scales. The isosurfaces of $\unicode[STIX]{x1D719}_{v}$ and $|\unicode[STIX]{x1D74E}|$ from $t=0$ to $2$ are shown in figure 9, when the flow evolves from a Gaussian random field into fully developed HIT. At $t=0$ , the isosurfaces of $\unicode[STIX]{x1D719}_{v}$ and $|\unicode[STIX]{x1D74E}|$ exhibit very different geometries, i.e. curved sheet-like structures and scattered blob-like structures in figures 9(a) and 9(b), respectively. The sheet-like vortex surfaces are rolled up into tube-like structures at $t=1$ , which is consistent with the observations of decaying HIT in Vincent & Meneguzzi (Reference Vincent and Meneguzzi1994) using vorticity vectors and with the transitional stage in TG flows in Yang & Pullin (Reference Yang and Pullin2011) using VSF.

The vortex tubes are then stretched and their characteristic scales decrease until the viscous effect becomes dominant to suppress the scale cascade. The smallest scale of vortex surfaces in figure 9 occurs around $t=2$ , which also corresponds to the time of the maximum dissipation rate in figure 1(b). The differences of the morphology of vortex surfaces from $t=0$ to $2$ are also highlighted in the close-up views in figure 13. Instead of the isolated helical vortex tubes in TG flows in figure 7 and scattered short tubes of $|\unicode[STIX]{x1D74E}|$ in figure 9(f), the vortex tubes and sheets in HIT generally connect and tangle with neighbouring ones to constitute a complex network in figure 9(e). Subsequently, the scale of vortex surfaces begins to increase during the late evolution stage after $t=10$ in figure 10 owing to the viscous effect (see Hao et al. Reference Hao, Xiong and Yang2019), when the turbulent kinetic energy has diminished below $1\,\%$ of the initial value. Moreover, appendix A explains that the statistical geometry of VSF isosurfaces is insensitive to the initial condition of (3.3), so the tangle of vortex tubes appears to be a global attractor in the chaotic vorticity system of fully developed HIT.

Figure 13.

Close-up views of the VSF isosurfaces (from left to right: $t=0$ , 1 and 2) around the central region in figures 9(a), (c) and (e).

The deformation of VSF isosurfaces is accompanied by intensive stretching, folding and twisting of the attached vortex lines with increasing curvature and torsion. The local curvature and torsion of vortex lines are calculated from instantaneous $\unicode[STIX]{x1D74E}$ as

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D705}=\frac{|\unicode[STIX]{x1D74E}\times (\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D74E})|}{|\unicode[STIX]{x1D74E}|^{3}}\end{eqnarray}$$

and

(5.2) $$\begin{eqnarray}{\mathcal{T}}=\frac{\unicode[STIX]{x1D74E}\boldsymbol{\cdot }[\unicode[STIX]{x1D735}(\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D74E})]\boldsymbol{\cdot }[\unicode[STIX]{x1D74E}\times (\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D74E})]}{|\unicode[STIX]{x1D74E}\times (\unicode[STIX]{x1D74E}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D74E})|^{2}},\end{eqnarray}$$

respectively. Figure 14 plots the temporal evolution of $\unicode[STIX]{x1D705}$ and $|{\mathcal{T}}|$ averaged over the isosurface of $\unicode[STIX]{x1D719}_{v}=0$ . The averaged curvature and torsion peak at $t\approx 2$ , which coincides with the occurrence of the maximum distortion of VSF isosurfaces in figures 9 and 10 and the maximum dissipation rate in figure 1(b). Hence, we speculate that the stretching and twisting of vortex surfaces/lines are correlated to the energy cascade in HIT, which is also supported by the VSF characterization of filtered HIT fields in appendix C.

Figure 14.

Temporal evolution of the curvature and absolute value of torsion averaged over the VSF isosurface of $\unicode[STIX]{x1D719}_{v}=0$ .

Additionally, the understanding of vortex morphology depends on the structure-identification method. For example, in figure 10(a,b) at $t=10$ , the isosurface of $\unicode[STIX]{x1D719}_{v}$ displays a complex network of vortex surfaces as a tangle of long vortex tubes, whereas the isosurface of $|\unicode[STIX]{x1D74E}|$ displays clusters of short tube-like structures as ‘vortex worms’. Sometimes the energy cascade is characterized by the visual ‘breakdown’ of isosurfaces of $|\unicode[STIX]{x1D74E}|$ , as shown in the right columns of figures 9 and 10 and illustrated by blue patches in figure 15. However, this breakdown of isosurfaces of $|\unicode[STIX]{x1D74E}|$ is elusive in terms of vortex dynamics, because physically reasonable deformation of vortex lines and surfaces in energy cascade should be due to continuous processes such as stretching, twisting and reconnection, as illustrated by red patches in figure 15.

Figure 15.

A schematic diagram of different views of the scale cascade of vortical structures in turbulence. Blue and red patches denote the isosurfaces of vorticity magnitude and VSF, respectively. Red lines with arrows denote vortex lines, dashed black lines with arrows denote the velocity induced by the vortex lines attached to vortex surfaces and grey arrows denote the background straining motion.

Therefore, as sketched in figure 15, the isosurface of $\unicode[STIX]{x1D719}_{v}$ appears to depict a more comprehensive picture than the isosurface of $|\unicode[STIX]{x1D74E}|$ of vortex dynamics in energy cascade, e.g. vortex twisting at small scales and rolling-up of vortex sheets. A similar visual difference between isosurfaces of $\unicode[STIX]{x1D719}_{v}$ and $|\unicode[STIX]{x1D74E}|$ is observed in TG flows in figure 7. For comparison, HIT has less isolated twisted vortex tubes in figures 9 and 10, because vortex surfaces in HIT are superposed and tangled together with strong interactions.