2.1. Initial condition

The vortex filament of a trefoil knot is described via the parametric curve:

(2.1)\begin{equation} \boldsymbol{X}(\theta) = \begin{pmatrix} R_{min}\left(\sin(\theta)+2\sin(2\theta)\right) \\ R_{min}\left(\cos(\theta)-2\cos(2\theta)\right) \\ -R_{min}\sin(3\theta) \end{pmatrix} \end{equation}

with the local knot radius given by

(2.2)\begin{equation} R(\theta) = \sqrt{R_{min}^2 \sin ^2(3 \theta)+R_{min}^2(5-4\cos(3\theta)) }. \end{equation}

Here $R_{min}$ is the minimum radius, and also the knot half-thickness in the propagation direction (figure 1). The mean radius, defined as

(2.3)\begin{equation} \bar{R}=\frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}R(\theta) \,\textrm{d} \theta, \end{equation}

will be used as a characteristic length scale in the present study, where we also choose the knot half-thickness to be $R_{min}$, resulting in $R_{max}=3R_{min}$, $\bar {R}=2.243R_{min}$ or $\bar {R}=0.748R_{max}$. All of the geometrical parameters used in the present study are illustrated and summarized in figure 1 and table 1.

Table 1. Geometrical parameters defining the knotted vortex filament and initial effective viscous core radius, $r_e$.

The velocity field induced by the vortex filament is initialized by numerically integrating the Biot–Savart law:

(2.4)\begin{equation} \boldsymbol{u}(\boldsymbol{x})=-\frac{\varGamma}{4{\rm \pi}} \int{K_v\frac{\left(\boldsymbol{x}-\boldsymbol{X}(\theta)\right)\times \boldsymbol{t}(\theta)}{\left|\boldsymbol{x}-\boldsymbol{X}(\theta)\right|^3}\textrm{d}\theta}, \end{equation}

where $\boldsymbol {t}(\theta )$ is the unit vector tangent to the filament, $\varGamma$ is the circulation and $K_v$ is a regularization function taken from Vatistas’ vortex-core model (Vatistas, Kozel & Mih Reference Vatistas, Kozel and Mih1991). The expression for $K_v$ is proposed by Van Hoydonck, Bakker & Van Tooren (Reference Van Hoydonck, Bakker and Van Tooren2010) for the vortex-core tangential velocity as a function of the radial distance from the core centre:

(2.5)\begin{equation} K_v=\frac{\left|\boldsymbol{x}-\boldsymbol{X}(\theta)\right|^3} {{\left(\left|\boldsymbol{x}- \boldsymbol{X}(\theta)\right|^{4}+{r_c}^{4}\right)}^{{3}/{4}}}. \end{equation}

Here $r_c$ is the nominal vortex tube radius, which we relate to the effective core radius $r_e$ via $r_e = 2r_c$. With this choice, a vortex tube of diameter $2r_e$ accounts for 95 % of the vorticity carried by the tube. With this kernel function choice, a Gaussian profile can be obtained in the initialized vortex-core structure. Core thickness is a very crucial factor affecting the overall dynamics as well as the small scales for this problem. The purpose of this paper is to report the numerical simulation results of the evolution of a thin-core vortex knot, whose thickness is inspired by Irvine's experimental investigations of a knotted vortex (Kleckner & Irvine Reference Kleckner and Irvine2013).

Simulations are performed in a triply periodic cubic domain with dimensions $\varOmega =[0,L]^3$, with the trefoil knot vortex initialized at the centre ($0.5L,0.5L,0.5L$) with $L=30R_{min}$ for AMR simulation and $L=15R_{min}$ for LES. Appendix A provides the analysis on the sensitivity of the box size, showing that $L=15R_{min}$ is a sufficient box size for the current study. All simulations are run with circulation $\varGamma =0.02$ m$^2$ s$^{-1}$ and $r_e/ \bar {R}=0.059$ (see table 1).

The circulation is set to $\varGamma =0.02$ m$^2$ s$^{-1}$ for all cases and changes in Reynolds number $Re_{\varGamma }=\varGamma /\nu$ are achieved by solely changing the viscosity. The ratio of specific heats $\gamma =1.4$, reference density $\rho _{ref} = 1.0$ and reference temperature $T_{{ref}}=1.0$ are used to compute the reference pressure $p_{{ref}}=\rho c^2/ \gamma$. In the following text, a non-dimensional time $t^*=t\varGamma /\bar {R}^2$ is used to display transient data.

2.2. Governing equations

The flow motion considered is assumed to be governed by the set of compressible Navier–Stokes equations,

(2.6)\begin{equation} \mathcal{NS}(\boldsymbol{{w}})=\frac{\partial{\boldsymbol{{w}}}}{\partial{t}}+\boldsymbol{{\nabla}}\boldsymbol{{\cdot}} \left[\boldsymbol{\mathsf{F}}_{{c}}(\boldsymbol{{w}})-\boldsymbol{\mathsf{F}}_{{v}}(\boldsymbol{{w}},\boldsymbol{{\nabla}}\boldsymbol{{w}})\right]=\boldsymbol{{0}}, \end{equation}

where $\boldsymbol {{w}}=(\rho ,\rho \boldsymbol {{U}},\rho E)^{\mathrm {T}}$ is the vector of conserved variables $\rho$, $\boldsymbol {{U}}$ and $E$, density, velocity and total energy, respectively, and $(\boldsymbol {{\nabla }}\boldsymbol {{w}})_{ij} = \partial {w_i}/\partial {x_j}$ its gradient. The viscous and convective flux tensors $\boldsymbol {\mathsf {F}}_{{c}},\boldsymbol {\mathsf {F}}_{{v}}\in \mathbb {R}^{5\times 3}$ read as

(2.7a,b)\begin{equation} \boldsymbol{\mathsf{F}}_{{c}}=\begin{pmatrix} \rho\boldsymbol{{U}}^{\mathrm{T}}\\ \rho\boldsymbol{{U}}\otimes\boldsymbol{{U}}+ p\boldsymbol{\mathsf{I}}\\ (\rho E+p )\boldsymbol{{U}}^{\mathrm{T}} \end{pmatrix}\quad \mathrm{and}\quad \boldsymbol{\mathsf{F}}_{{v}}=\begin{pmatrix} \boldsymbol{{0}}\\ \boldsymbol{\tau}\\ \boldsymbol{\tau}\cdot\boldsymbol{{U}}-\lambda\boldsymbol{{\nabla}} T^{\mathrm{T}} \end{pmatrix}, \end{equation}

where $T$ is the temperature, $p$ is the pressure, $\lambda$ is the thermal conductivity of the fluid and $\boldsymbol {\mathsf {I}} \in \mathbb {R}^{3\times 3}$ is the identity matrix. For a Newtonian fluid, we have

(2.8)\begin{equation} \boldsymbol{\tau}=2\mu\boldsymbol{\mathsf{S}}, \end{equation}

where $\mu$ is the dynamic viscosity and

(2.9)\begin{equation} \boldsymbol{\mathsf{S}}=\tfrac{1}{2}\left[\boldsymbol{{\nabla}}\boldsymbol{{U}}+\boldsymbol{{\nabla}}\boldsymbol{{U}}^{\mathrm{T}}-\tfrac{2}{3} \left(\boldsymbol{{\nabla}}\boldsymbol{{\cdot}}\boldsymbol{{U}}\right)\boldsymbol{\mathsf{I}}\right].\end{equation}

The ideal gas law is considered for the closure of the system of equations, namely,

(2.10)\begin{equation} p=(\gamma-1)\left(\rho E-\tfrac{1}{2}\rho\boldsymbol{{U}}\cdot\boldsymbol{{U}}\right), \end{equation}

where $\gamma$ is the heat capacity ratio.

2.3. Adaptive mesh refinement computational framework: the VAMPIRE code

The DNS of the knotted vortex reconnection presented in this work are performed with the AMR high-order compact finite difference code VAMPIRE (Zhao & Scalo Reference Zhao and Scalo2020). The AMR framework of the VAMPIRE code is based on Paramesh (MacNeice et al. Reference MacNeice, Olson, Mobarry, De Fainchtein and Packer2000). For each simulation, before time advancement, a sensor estimating the discretization error of the sampled initial condition is ran on all the blocks to determine whether the mesh should be locally refined; such a process is repeated until convergence of the AMR tree structure. A sample of the resulting adaptively refined mesh at time $t^*=0$, as well as subsequent times, is shown in figure 2 for $Re_\varGamma =6000$.

Figure 2. Grid following the vortex transport in the DNS simulation with AMR. The upper row and lower row show different views of the simulation results for $Re_\varGamma =6000$.

In this work our goal is to apply the AMR code to fully resolve the vortex reconnection dynamics. Therefore, we choose to use a vorticity-based sensor, defined in the non-dimensional form as

(2.11)\begin{equation} {f_\omega } = \frac{{ |\boldsymbol{{\nabla}} (\omega )| \varDelta}}{{\max \left( {{{\left\| \omega \right\|}_\infty }, {\omega _{\textrm{{ref}}}}} \right)}}, \end{equation}

where $\varDelta$ is the local grid spacing, and $\omega$ is the magnitude of vorticity, i.e. $\omega := |\boldsymbol {\omega }|$. This expression is the normalized gradient of vorticity by a local indicator of the vorticity magnitude. From a numerics perspective to avoid division by zero, an extra term $\omega _{{ref}}$ is added to the denominator of the expression.

This reference vorticity value is expected to be a measure for the global level of vorticity magnitude. The sensor is expected to refine the grid in the region where the vortical features are significant, while coarsening the grid at the regions far away from the vortex ring, where the vorticity vanishes. More details of the numerics can be found in Zhao & Scalo (Reference Zhao and Scalo2020). The AMR structure is updated every time step. All of the AMR simulation results presented in this paper are run at ${CFL}=0.25$.

2.4. Large-eddy simulation framework on Cartesian grid

In the present study, the compressible LES formalism introduced by Lesieur & Metais (Reference Lesieur and Metais1996) and coworkers (Lesieur & Comte Reference Lesieur and Comte2001; Lesieur, Métais & Comte Reference Lesieur, Métais and Comte2005) is adopted yielding the following set of filtered compressible Navier–Stokes equations:

(2.12)\begin{equation} \mathcal{NS}(\bar{\boldsymbol{{w}}}) = \boldsymbol{{\nabla}} \boldsymbol{{\cdot}} \boldsymbol{\mathsf{F}}_{{SGS}} (\bar{\boldsymbol{{w}}},\boldsymbol{{\nabla}}\bar{\boldsymbol{{w}}}). \end{equation}

Here $\bar {\boldsymbol {{w}}}=(\bar {\rho },\bar {\rho }\widetilde {\boldsymbol {{U}}},\bar {\rho }\widetilde {E})^{\mathrm {T}}$ is the vector of filtered conservative variables. The LES equations are obtained by applying a low-pass filter to the Navier–Stokes equations (Leonard Reference Leonard1974). The spatial filtering operator applied to a generic quantity $\phi$ reads as

(2.13)\begin{equation} \bar{\phi}(\boldsymbol{{x}},t)=g (\boldsymbol{{x}})\star\phi (\boldsymbol{{x}}), \end{equation}

where $\star$ is the convolution product and $g(\boldsymbol {{x}})$ is a filter kernel related to a cutoff length scale $\bar {\varDelta }$ in physical space (Sagaut Reference Sagaut2006). The compressible case requires density-weighted (or Favre) filtering operators, defined as

(2.14)\begin{equation} \widetilde{\phi}=\frac{\bar{\rho\phi}}{\bar{\rho}}. \end{equation}

The SGS tensor $\boldsymbol {\mathsf {F}}_{{SGS}}$ is the result of the filtering operation, it encapsulates the dynamics of the unresolved SGS, and is modelled here using the eddy-viscosity assumption:

(2.15)\begin{equation} \boldsymbol{\mathsf{F}}_{{SGS}}(\bar{\boldsymbol{{w}}},\boldsymbol{{\nabla}}\bar{\boldsymbol{{w}}}) =\begin{pmatrix} \boldsymbol{{0}}\\ 2\mu_{t}\bar{\boldsymbol{\mathsf{S}}}\\ -\dfrac{\mu_{t}C_{p}}{Pr_{t}}\boldsymbol{{\nabla}}\widetilde{T}^{\mathrm{T}} \end{pmatrix}. \end{equation}

Here $\bar {\boldsymbol {\mathsf {S}}}$ is the shear stress tensor computed from (2.9) based on the Favre-filtered velocity $\widetilde {\boldsymbol {{U}}}$, $Pr_{t}$ is the turbulent Prandtl number, which is set to 0.5 (Erlebacher et al. Reference Erlebacher, Hussaini, Speziale and Zang1992), $C_{p}$ is the heat capacity at constant pressure of the fluid and $\mu _{t}$ is the eddy viscosity, which depends on the chosen subgrid model.

In the present work, the CvP-Smagorinsky closure is adopted, which yields accurate results for transitional and high-Reynolds-number flows. The LES methodology employs the CvP approach introduced by Chapelier et al. (Reference Chapelier, Wasistho and Scalo2018), which has been validated in the context of transitional helical vortices simulations (Chapelier et al. Reference Chapelier, Wasistho and Scalo2019). The expression for eddy viscosity for this closure reads as

(2.16)\begin{equation} \mu_t=\rho f(\sigma) (C_{\mathrm{S}}\bar{\varDelta})^2\sqrt{2S_{ij}S_{ij}}, \end{equation}

where $C_{\mathrm {S}}=0.172$ is the Smagorinsky constant, $\bar {\varDelta }$ is the local (linear) grid size and $f(\sigma )$ is the CvP turbulence sensor built from the ratio of the test-filtered to grid-filtered enstrophy,

(2.17)\begin{equation} \sigma=\dfrac{\;\;\widehat{\bar{\xi}}\;\;}{\;\;\bar{\xi}\;\;}, \end{equation}

where $\bar {\xi }= \bar {\boldsymbol {\omega }} \cdot \bar {\boldsymbol {\omega }}/2$. Near-unitary values of such ratio, $\sigma \simeq 1$, imply a low degree of local spectral broadening of the flow (for example, in regions of transitional turbulence or coherent laminar vortices) and, hence, result in the local attenuation of SGS dissipation.

The compressible, Favre-filtered Navier–Stokes equations are solved on a simple uniform Cartesian grid using a sixth-order compact-finite-difference scheme solver originally written by Nagarajan, Lele & Ferziger (Reference Nagarajan, Lele and Ferziger2003), under continued development at Purdue University. The solver is based on a staggered grid arrangement, providing superior accuracy and robustness compared with a fully collocated approach (Lele Reference Lele1992). The time integration is performed using a third-order Runge–Kutta scheme.

For the presently considered computations, the speed of sound is arbitrarily selected as $c_{ref}=5V_{{max}}$, where $V_{{max}}$ is the maximum velocity magnitude at the initial condition. The maximum Mach number is achieved during reconnection and it does not exceed 0.25 in all cases simulated. This choice yields near-incompressible flow dynamics. Simulations are initialized with unitary density and temperature fields, equal to the reference values $\rho _{ref} = 1.0$ and $T_{{ref}}=1.0$. The gas constant is then defined by setting the reference pressure to $p_{{ref}}=\rho _{ref} c_{ref}^2/ \gamma$, where $\gamma =1.4$ is the ratio of specific heats. Levels of density, temperature and pressure are dimensionless and physically irrelevant, as the compressibility effects on the simulated hydrodynamics are negligible.

3.1. Grid convergence analysis of AMR calculations

A grid sensitivity analysis has been carried out for all Reynolds numbers considered (figure 3, table 2) based on the globally integrated enstrophy

(3.1)\begin{equation} \varXi = \iiint_\varOmega{ \xi\,\textrm{d} V}. \end{equation}

This quantity is highly sensitive to small-scale dynamics and a grid-converged enstrophy evolution is therefore indicating that the resolution is sufficient to capture the small scales generated during the reconnection process. We hereafter describe that state of the AMR tree by reporting the block tree depth, $d$, and number of grid points per block $N_b$. Grid configurations used for the two circulation-based Reynolds numbers considered in this paper include 8 and 9 levels with $12^3$, $18^3$, $24^3$, $30^3$, $36^3$ grid points per block. As observed from the simulation results, with the refinement sensor specified in § 2.3, all the regions containing the vortex core are initially refined to the finest level; therefore, the resolution provided by the grid configuration can be represented by the number of points inside the initial vortex-core diameter ($2r_e/\varDelta$). The enstrophy curves from the two finest simulations in figure 3 overlap each other. An AMR grid with 9 refinement levels, $12^3$-points and $30^3$-points per block provide grid-converged global enstrophy for the ${Re}_{\varGamma } = 2000$ and 6000 cases, respectively. Figure 3 shows that, for both Reynolds numbers, starting from the same initial condition, the global enstrophy decreases first due to viscous dissipation. The dissipation rate is relatively stronger for ${Re}_{\varGamma } = 2000$. During reconnection, both enstrophy curves increase and reach their peak, followed by a gradual decay. The peak value of global enstrophy for ${Re}_{\varGamma } = 6000$ is significantly higher than the ${Re}_{\varGamma } = 2000$ case, and even higher than its initial value, in this case due to the production of very small vortical scales during the reconnection. Whether the latter constitutes the occurrence of turbulence in a state of quasi-equilibrium or not, a proper dynamic LES model tasked with simulating such a highly unsteady and inhomogeneous flow should not vitiate the underlying coherent vortex dynamics by adding unnecessary amounts of eddy viscosity.

Figure 3. Grid convergence analysis for AMR simulations: $Re_{\varGamma }=2000$ (green), $Re_{\varGamma }=6000$ (blue), showing the evolution of the resolved volume-integrated enstrophy $\varXi$ (or, equivalently, normalized total viscous dissipation, $\epsilon$). The results in the plot include the six AMR simulations at $Re_{\varGamma }=6000$ and the three AMR simulations for $Re_{\varGamma }=2000$ (see table 2) with increasing levels of grid refinement as shown by the arrow.

Table 2. Summary of flow and grid parameters for the (single-mesh-level, i.e. $d=1$) Cartesian CvP-LES and the AMR runs, executed on $d=8$ and $d=9$ mesh levels. Here $N^3$ indicates the total number of points per AMR block, or in the whole computational domain for the Cartesian LES, $2r_e/\varDelta$ is an estimate of the number of grid points per initial effective vortex diameter, and $N_{{tot}}$ is the total degrees of freedom for each simulation in the unit of millions. For AMR simulations, $N_{{tot}}$ is taken at $t^*=3.0$ right after the reconnection.

3.2. Grid sensitivity study for LES with Cartesian grid

A grid sensitivity analysis of the LES computations has been carried out for all Reynolds numbers considered (figure 4, table 2). Computational grids with $N^3= 120^3,\ 192^3,\ 288^3$ and $432^3$ points are considered, corresponding to, respectively, 3.4, 3.5, 5.1 and 7.6 points inside the initial vortex-core diameter ($2r_e$). The adequacy of the adopted grid resolution is also assessed by monitoring the level of SGS, or modelled, dissipation (figure 4). Such quantities are not intended to provide a rigorous statistical representation of the state of turbulence (when and if present at all), since the flow is highly inhomogeneous and unsteady. Such quantities will be merely used to infer the sensitivity of the resolved flow to the grid and the robustness of the adopted CvP-LES closure. As discussed in the following, the flow evolution is in fact laminar (albeit highly unsteady) before and after reconnection, with near-equilibrium turbulence occurring only during reconnection and at higher Reynolds numbers (figure 5), posing an important numerical modelling challenge for the CvP-Smagorinsky closure. As indicated in table 2, the CvP LES conducted in this study are computationally inexpensive compared with DNS with AMR; however, as shown in the rest of the paper, the CvP-LES correctly predicts key features of the flow, such as the jump in global helicity discussed in § 6.1, demonstrating its suitability as a predictive tool for high-Reynolds-number vortical flow evolution. Figure 5 also shows the comparison between one-dimensional spectra extracted in the $z$ direction from the CvP-LES and AMR simulations. The overlap in the low wavenumber range indicates that the CvP-LES correctly captures the effects of the unresolved small scales onto the resolved large scale despite the highly inhomogeneous and non-equilibrium nature of turbulence in this flow.

Figure 4. Grid sensitivity analysis for LES: $Re_{\varGamma }=2000$ (green), $Re_{\varGamma }=6000$ (blue) showing the evolution of the resolved volume-integrated enstrophy, $\varXi$ (a), normalized total dissipation, $\epsilon _{total}$ (b), and SGS to total dissipation ratio (c). Grid resolutions considered are the finest three datasets from each simulation: coarse $( \cdots )$; intermediate $(\textrm {- - -})$; fine $(\textrm {------})$ (see data in table 2). Reference data shown in circles ($\circ$) are taken from the finest DNS data.

Figure 5. Grid sensitivity analysis for Cartesian CvP-LES (black) and AMR (green/blue) simulations: $Re_{\varGamma }=2000$ (a), $Re_{\varGamma }=6000$ (b) showing the one-dimensional energy spectra in the $z$ direction at the time of first reconnection ($t_1$) defined as the time when the maximum dissipation is achieved (see figure 2). Grid resolutions considered are the finest three datasets from each simulation: coarse $( \cdots )$; intermediate $(\textrm {- - -})$; fine $(\textrm {------})$ (see data in table 2). Red dotted line shows a reference $-$5/3 slope. Here $L_z$ is the domain size in $z$ direction.

The volume-integrated resolved enstrophy (figure 4a) increases monotonically with Reynolds numbers and grid resolution, as expected. The peak in total dissipation (figure 4b), which comprises the modelled $\epsilon _{SGS}$ and the resolved, occurring upon reconnection, varies modestly for all grids considered at any given Reynolds number, indicating a healthy response of the CvP-LES closure. Velocity-fluctuation intensity production only occurs upon vortex reconnection, responsible for the generation of small-scale vortical structures (see § 5), where a non-vanishing SGS energy content is detectable even for the lowest $Re_{\varGamma }$ at the highest grid resolution available (figure 4c). An acceptable level of grid convergence on the total dissipation before and after reconnection is achieved for both Reynolds numbers, which can be seen from figure 4(c) in the comparison with the grid-converged AMR results. The high levels of subgrid dissipation at the beginning of the computation (${\sim }70\,\%$ of the total dissipation for the coarser grid) can be explained as the CvP model detects marginally resolved initial vortex cores and proceeds to make them smoother by the action of subgrid dissipation. This process does not impact the subsequent total dissipation evolution ($t^*>1$) which is remarkably close to the DNS results for both Reynolds numbers considered.

5.1. Anti-parallel vortex pair $\langle I \rangle$

Vortex regions undergoing the highest rate of axial stretching align right before reconnection in an anti-parallel fashion (figure 11, $t^*=2.42$) and viscously interact. The vorticity in the outer layers of the vortex tube undergoes cancellation through viscous diffusion, forming hairpin vortices acting as bridges (figure 11, $t^*=2.62$); the small-scale structures left over from the reconnection process, named threads, form weak topological connections with the two newly formed vortex rings (figure 11, $t^*=2.81$), and eventually dissipate.

These dynamics are similar to those observed in Crow instability (Crow Reference Crow1970), which is the most studied scenario for anti-parallel vortex reconnection, investigated in the context of contrail vortices behind an aircraft. The difference with respect to the knotted vortex lies in the source of perturbations. The Crow instability is a mutual induction sinusoidal instability through which the two anti-parallel vortices interact. However, in the present study, the reconnection process is triggered by the vortex self-advection and distortion.

5.2. Hairpin vortex formation via bridging $\langle II \rangle$

Figure 13 illustrates the detail of the reconnection process for $Re_{\varGamma }=2000$ and $Re_{\varGamma }=6000$. While the interaction during the anti-parallel reconnection results in a cancellation of the outer layer vorticity, enhancement of the inner layer one (in the points of contact), and thinning of the vortex tubes (structure 1), the remaining vorticity unaffected by the cancellation is responsible for initiating the bridging process (structure 2 at $t^*=3.16$), yielding intermediate hairpin vortex structures on the two extremities of the thinning anti-parallel structures – the leading and following extremity, based on the self-propagation direction of the knot (i.e. along ${+}z$). Parts of newly formed hairpin vortices merge into structure 3, spanning both the anti-parallel structure and early hairpin vortex (structure 2). With structures 1 and 2 gradually dissipating, structure 3 becomes dominant with a vortex-core thickness comparable to the original structure and ultimately responsible for completing the reconnection process. This type of reconnection was named bridging by Kida & Takaoka (Reference Kida and Takaoka1994), to be distinguished from the vortex cancellation process described in Fohl & Turner (Reference Fohl and Turner1975) and Oshima & Asaka (Reference Oshima and Asaka1977). The various steps involved in this process are clearly visible at $Re_{\varGamma }=2000$, while at $Re_{\varGamma }=6000$, the reconnection results in the generation of smaller scale structures. The same process occurs on both ends of the anti-parallel structure, resulting in the initialization of axial flow and the travelling of helical structures, or Kelvin waves, on the small and large vortex rings after separation. The curvature of the bridges causes the repulsion from the reconnection zone due to the induced velocity field, which further stretches the threads.

Figure 13. Visualization of the generation of hairpin vortices from $Q$-isosurfaces at $Q\bar {R}^4/\varGamma ^2=96$ for ${Re}_{\varGamma }=2000$ and ${Re}_{\varGamma }=6000$ from DNS/AMR data, coloured by normalized helicity density; threads created during the bridging process, intermediate hairpin vortex structure the final dominant hairpin structure are illustrated.

5.3. Thread structures $\langle III \rangle$

Threads are generated at the beginning of the reconnection process as shown in figures 14 and 15 and figure 16 at around $t^*=2.90$ for $Re_{\varGamma } = 2000$ and $t^*=2.63$ for $Re_{\varGamma } = 6000$. The thread structure is the byproduct of the circulation transfer process from one vortex segment to its anti-parallel neighbour, thus, its generation coincides with the formation of the bridge comprising the hairpin vortex. Similarly to the observation reported by Hussain & Duraisamy (Reference Hussain and Duraisamy2011), the two anti-bridges curl back and retract from one another shortly after they are formed. As the threads in between are stretched, the induced velocity by the bridges tend to change the curvature of the threads, further distancing them from one another, which can be noticed especially at the thread-bridge junctions at $Re_{\varGamma } = 2000$. For the higher Reynolds number case $Re_{\varGamma } = 6000$, the stretching is more intense so that the anti-parallel structure is squeezed into a thinner film and subsequently breaks down into finer scales due to Kelvin–Helmholtz instability. In contrast to the case $Re_{\varGamma } = 2000$ where threads are side by side, two major thread structures in $Re_{\varGamma } = 6000$ can be identified, one on the top and the other at the bottom, with the thin-film vorticity remaining in between, as shown in figure 15. The majority of the high vorticity is concentrated at the thread on the top only, and at a high Reynolds number features significant subsequent local reconnections (i.e. cascaded reconnection), which results in small-scale vorticity surrounding the thread on the top. The structure of the small-scale vortices generated by this reconnection is similar to the cascaded bridgelet structure discussed by Yao & Hussain (Reference Yao and Hussain2020b) (C-structure in figure 7(e) of the cited paper), in which the mechanism of the cascade reconnection has been outlined for the $Re_{\varGamma }=9000$ case. The difference with the current work is that the strain field in the current paper is asymmetric, resulting in a more asymmetric evolution of the bridgelets (as shown on the lower right of figure 13). The reversal of curvature near the thread-bridge junctions at $Re_{\varGamma } = 6000$ can also be observed, but it manifests itself in the form of bifurcations of the threads together with the generation of numerous smaller scale structures (as seen in the lower right parts of figure 13).

Figure 14. Isosurfaces of vorticity magnitude at various times during reconnection: $\omega = 0.3\varGamma /{\rm \pi} r_{e,0}^2$ coloured in transparent grey and $\omega = 1.0\varGamma /{\rm \pi} r_{e,0}^2$ coloured by red at $Re_{\varGamma } = 2000$.

Figure 15. Isosurfaces of vorticity magnitude at specific time steps during reconnection: $\omega = 0.5\varGamma /{\rm \pi} r_{e,0}^2$ coloured in transparent grey and $\omega = 3.0\varGamma /{\rm \pi} r_{e,0}^2$ coloured by red at $Re_{\varGamma } = 6000$.

Figure 16. Visualization of the threads (marked by dashed rectangles) evolution using $Q$-isosurfaces with $Q\bar {R}^4/\varGamma ^2=323$ for ${Re}_{\varGamma }=6000$ DNS coloured by normalized pressure fluctuations (structure $\langle III \rangle$) in figure 11.

Along with the separation into leading small and trailing large vortex rings, the threads become stretched in the self-propagation direction $z$. At the same time, the threads begin to diverge along the large vortex ring due to the convection of the flow around and along vortex tubes ($t^*=3.80$).

The main vortex ring structure keeps rotating around the $z$-axis and the longitudinally stretched threads separate into two parts: one part is convected by the triangular-shaped vortex rings, then wraps around the bridges. This part will stay connected to the reconnection points and rotates together with the two separate rings around the $z$ direction. The other part loses the original connection to the main vortex ring structures and stays at its old position in the $x$–$y$ plane, as seen in the dashed boxes at $t^*=4.52$ and $5.22$. The large vortex ring will rotate across such stagnating threads, which serves as a source of perturbation that excites more helical structures at the outer layers of vortex rings, as discussed in the next section.

Furthermore, additional asymmetric features of the reconnection in this problem include the relative rotation between two bridges, as can be seen in figure 10(c–e). At $Re_{\varGamma } = 6000$, the lower bridge, which is part of the larger trailing ring after separation, undergoes a torsional displacement with respect to the threads (i.e. with respect to the Y’-axis defined in the moving local coordinate in figure 10), while the upper bridge on the small ring remains nearly orthogonal to the threads through the entire reconnection process, and the vertical film formed by the threads remains straight. In contrast, at $Re_{\varGamma } = 2000$ the two side-by-side threads are further distorted by the relative motions between bridges (see figure 10b). This is due to the asymmetric configuration of the vortex tube and the asymmetric and non-uniform distribution of the axial flow packets at the onset of reconnection.

5.4. Helical structures and Kelvin waves $\langle IV \rangle$

The hairpin vortex structures generated after reconnection are responsible for exciting helix-shaped structures (figure 17), corresponding to two Kelvin-wave packets propagating away from the reconnection region, in opposite directions and with opposite handedness. These structures, coloured in the figure by the helicity density (5.1), are advected by the axial flow, which occurs when the flow velocity vector inside a vortex core is aligned with the vortex line, thus, the positive values denote a propagation along the direction of the vorticity vector (vice versa, if negative).

Figure 17. $Q$-isosurfaces with $Q\bar {R}^4/\varGamma ^2=323$ from DNS results at $Re_\varGamma =6000$. Panels (a,b) show the axial flow direction before reconnection; (c,d,b) depicts the Kelvin-wave packet generated and colliding on the leading small ring, and (f–i) show the Kelvin-wave packets travelling and colliding on the large vortex ring after reconnection; (h) shows when the Kelvin-wave packets collide on the large vortex ring, an annular structure is formed; (j) shows the centreline axial flow energy spectrum before, during and after the Kelvin-wave collision event happened on the large vortex ring as shown in (g–i), respectively.

The helical structures travelling along each single side of the triangular-shaped vortex ring periodically collide with each other, then continue to travel past one another. Each collision event exhibits features of vortex bursting as defined by Spalart (Reference Spalart1998), who investigated this phenomenon in the context of contrail vortices in airplane wakes. Vortex bursting is also observed in the reconnection of anti-parallel vortex tubes (Crow instability), after which the axial flow is driven away from the reconnection region as observed in Misaka et al. (Reference Misaka, Holzäpfel, Hennemann, Gerz, Manhart and Schwertfirm2012) and Van Rees et al. (Reference Van Rees, Hussain and Koumoutsakos2012).

When the helical structures collide for the first time, the axial velocities cancel each other and form an annular structure (visualized in figure 17(h) by a local cut plane across the vortex ring). Similar structures and associated dynamics were observed and discussed by Melander & Hussain (Reference Melander and Hussain1994). When such an annular structure propagates along the vortex tube, it develops additional twists, which represent higher helical modes. Spalart (Reference Spalart1998) argues that vortex bursting does not compromise the coherent structure if the circulation-preservation principle is respected, and that experimental smoke visualizations cannot properly track the vortical structure during bursting.

In the present study such collisions occur as the helical, or Kelvin, waves are travelling along the vortex loop with nonlinear interactions. From a spectral energy content standpoint, the collision process transfers the energy from low wavenumbers to higher ones. One way to check this effect is to perform a Fourier analysis of the axial flow energy spectrum on one vortex loop, defined as

(5.2)\begin{equation} E_{\tilde{w}}(k,t)=\frac{1}{\ell (t)}\int_\ell \rho \tilde{w}^2|_{r=0} \exp({-j2{\rm \pi} k/\ell\tilde{z}})\,\textrm{d}\tilde{z}, \end{equation}

where $\tilde {z}$ is the position along the axial direction of the vortex line, and the axial flow velocity $\tilde {w}(\tilde {z},r,t)$ is defined as

(5.3)\begin{equation} \tilde{w}(\tilde{z},r,t) = \boldsymbol{u} \cdot {\tilde{\boldsymbol{z}}}, \end{equation}

where $\boldsymbol {{\tilde {z}}}$ is the unit tangential vector of the vortex line, which is normal to the local cut plane, and the boundary of the vortex tube cross-sectional area. Figure 17(j) shows the axial flow energy spectrum before, during and after the first three simultaneous Kelvin-wave collisions on the large trailing ring (figure 17g–i), depicting that the energy is suppressed during the collision, and redistributed to a higher wavenumber after the collision. Also the Kelvin waves are dispersing as they propagate and collide, resulting in a homogeneous distribution of wave energy along the vortex loops. On the leading small ring, similar but more intense events occur, because of the smaller size of the ring. The evolution of helical structures on the small ring can also be observed and the positive and negative helicity parts alternate rapidly due to the shorter length allowing propagation. The small scales generated on the two separated vortex rings, and especially on the leading smaller ring (figure 17d–i), are the result of Kelvin-wave collisions (which result in an energy cascade into smaller scales), and the interaction between the main vortex tube and the threads wrapping around it (as explained in § 5.3). A similar phenomenon can be found in Van Rees et al. (Reference Van Rees, Hussain and Koumoutsakos2012), where energy cascade from nonlinear Kelvin-wave interaction, as well as deformation of the vortex lines on the resulting vortex rings after reconnection is observed. The energy cascade from Kelvin waves is conjectured to be different from Kolmogorov cascade and has been widely discussed by Kozik & Svistunov (Reference Kozik and Svistunov2004), Vinen, Tsubota & Mitani (Reference Vinen, Tsubota and Mitani2003) and Yepez et al. (Reference Yepez, Vahala, Vahala and Soe2009) for quantum turbulence and superfluid turbulence. In Van Rees et al. (Reference Van Rees, Hussain and Koumoutsakos2012) the thread structure stays in the colliding plane because of symmetry and does not wrap around the vortex rings, whereas in the current work, the wrapping of the threads around the vortex tube and the resulting interaction (pairing and merging) further deforms the vortex tubes and adds to the complexity of the small scales generated.

Overall, there are similarities between the reconnection dynamics of topologically complex thin-core vortex knots and of a symmetric parallel reconnection as discussed by previous works, such as Hussain & Duraisamy (Reference Hussain and Duraisamy2011), Van Rees et al. (Reference Van Rees, Hussain and Koumoutsakos2012) and Yao & Hussain (Reference Yao and Hussain2020b), etc. An asymmetric knot configuration, however, results in much more complex fine-scale vortex structures and non-trivial global helicity dynamics. Moreover, the asymmetry increases with increased Reynolds number, which can be reflected by the curvature change shown in figure 9.

6.1. Total helicity

For the two different Reynolds numbers considered in this study, before reconnection, the total helicity is approximately constant with a value $H/ \varGamma ^2 \approx -3.3$, in good agreement with the results of Kida & Takaoka (Reference Kida and Takaoka1988) and Scheeler et al. (Reference Scheeler, Kleckner, Proment, Kindlmann and Irvine2014) who report values of $H/ \varGamma ^2 \approx 3.28$ and $H/ \varGamma ^2 \approx 3.26$, respectively. The difference in sign is merely due to a trivial difference in the handedness of the trefoil knot, specifically, a change in the sign of the $z$ component in (2.1).

As shown in figure 18, for $Re_\varGamma =2000$ (laminar reconnection), a maximum in helicity is reached and a steady decay observed after that. For higher Reynolds numbers, where turbulent production occurs, a significant and more abrupt increase in helicity is observed, followed by a steady increase in helicity. Figure 18 also shows the total helicity results from the finest LES. The results indicate that the LES is capable of capturing the trend of changing of total helicity correctly, including the abrupt helicity production during reconnection and the gradual decay or increase after reconnection, even though no specific helicity-preserving SGS model is employed in this work.

Figure 18. Direct numerical simulation and companion LES data for the evolution of volume-integrated, or total, helicity for $Re_\varGamma =2000$ (green) and $Re_\varGamma =6000$ (blue) for the vortex knot configuration considered in the present study. Direct numerical simulation results are shown by solid lines and the finest LES by circles ($\circ$).

To further study this peculiar behaviour, the following helicity budgets are considered:

(6.4)\begin{equation} \frac{\textrm{d}}{\textrm{d} t}H = \iiint -2\nu \boldsymbol{\omega} \cdot \left({\boldsymbol{\nabla}} \times \boldsymbol{\omega} \right) \textrm{d} \varOmega. \end{equation}

Here $\boldsymbol {{\omega }} \cdot (\boldsymbol {{\nabla }} \times \boldsymbol {\omega } )$ is the super-helicity density. Just like kinetic energy, the global helicity is a conserved quantity in an inviscid flow. However, viscous effects may entail an increase or decrease of the helicity depending on the right-hand side of (6.4). Figures 19(a) and 19(b) show the total helicity budget at both Reynolds numbers considered in the DNS. The excellent agreement between the independently evaluated left- and right-hand sides of (6.4) indicates that the evolution of global helicity is well resolved in both cases. The closure of the global helicity budgets, in fact, requires all scales of the flow to be properly resolved, especially due to the three nested spatial derivative operations present in the super-helicity term on the right-hand-side of (6.4). In the LES with CvP-Smagorinsky closure, the eddy viscosity $\mu _t$ contributes to the viscous helicity dynamics, as shown in the budget equation

(6.5)\begin{equation} \frac{\textrm{d}}{\textrm{d} t}H = \underbrace{\iiint -2\nu \boldsymbol{\omega} \cdot \left( \boldsymbol{{\nabla}} \times \boldsymbol{\omega} \right) \textrm{d} \varOmega}_\textrm{resolved helicity dissipation} + \underbrace{\iiint -\frac{2\mu_t}{\rho} \boldsymbol{\omega} \cdot \left( \boldsymbol{{\nabla}} \times \boldsymbol{\omega} \right) \textrm{d} \varOmega}_\textrm{modelled SGS helicity dissipation}, \end{equation}

where the helicity $H$ is here intended as based on the resolved velocity and vorticity, and $\boldsymbol {\omega }$ is the resolved vorticity vector. Figures 19(c) and 19(d) show the total helicity budget for both Reynolds numbers considered in the LES. The SGS contribution to the helicity dissipation adds on to the resolved contribution, allowing the LES to predict a peak in total rate of helicity change $\textrm {d} H(t)/\textrm {d} t$ matching reasonably well with the DNS results (figure 18). As opposed to the DNS/AMR results, the helicity budgets from the LES are not exactly numerically closed, due to the very nature of the adopted LES approach, which exhibits a high wavenumber energy build up. Moreover, the staggered grid arrangement entails various interpolation operations which further propagate numerical error, resulting in a loss of accuracy in the direct evaluation of the super-helicity on the right-hand side of (6.5). These numerical errors, however, vanish when refining the grid, and become negligible only once and if DNS levels of grid resolution are achieved. An explicitly filtered LES approach is expected to yield an exact numerical closure of the budgets, and will be considered in future work. The reader is also directed to previous research (Li et al. Reference Li, Meneveau, Chen and Eyink2006) on the suitability of SGS closure in vortex dominated turbulent flows.

The finest LES resolution for the $Re_\varGamma =2000$ and $Re_\varGamma =6000$ cases is 288$^3$ and 432$^3$, respectively. When the grid resolution of the LES is increased, the relative importance of the resolved viscous contribution to the total helicity dissipation rate increases with respect to the SGS contribution.

Furthermore, it is interesting to investigate the spatial distribution of helicity and helicity production (super-helicity), and their evolution with time. This is depicted in figures 20 and 21 showing the $z$ direction distribution of helicity and production of helicity density averaged in the $xy$ plane, where $z$ is the propagation direction of the self-induced motion of vortex knot. The helicity and production of helicity density average in the $xy$ plane are defined as

(6.6)\begin{equation} \left\langle h ^*\right\rangle_{xy} = \frac{1}{L^2}\int_0^L \int_0^L h(x,y,z,t) \bar{R}^3/\varGamma^2\,{\textrm{d}x}\,{\textrm{d}y} \end{equation}

and

(6.7)\begin{equation} \left\langle \textrm{d} h^*/\textrm{d} t^* \right\rangle_{xy} = \frac{1}{L^2}\int_0^L \int_0^L -2\nu \boldsymbol{\omega} \cdot \left( {\boldsymbol{\nabla}} \times \boldsymbol{\omega} \right) \bar{R}^5/\varGamma^3 \,{\textrm{d}x}\,{\textrm{d}y},\end{equation}

respectively. Figure 20 shows that when the trefoil knot reaches reconnection, the trailing part, which becomes the larger trailing ring, carries the majority of the helicity; on the other hand, the leading part, that becomes the leading smaller ring later, carries almost zero helicity. From figure 17(b), it can be seen that the leading part presents a nearly triangular shape with a distribution of helicity with a nearly zero average. Figure 21 shows that the location of the abrupt helicity increase during the reconnection (reflected in figure 18) at the reconnection regions between the leading and trailing rings. After reconnection, for both the $Re_{\varGamma }=2000$ and $Re_{\varGamma }=6000$ case, the helicity decay is found to be occurring mainly on the large trailing ring (see figure 21). This decay is likely due to the competing effects of Kelvin-wave collision and viscous dissipation.

Figure 19. Global helicity budget at different Reynolds numbers from the DNS and finest LES available: (a) $Re_{\varGamma }=2000$ DNS, in green; (b) $Re_{\varGamma }=6000$ DNS, in blue; (c) $Re_{\varGamma }=2000$ LES, in green; (d) $Re_{\varGamma }=6000$ LES, in blue. The symbols correspond to: $(\circ )$ resolved helicity production term (first term on the right-hand side of (6.4) and (6.5)); $(\triangle )$ SGS helicity production term (second term in (6.5)); $(\textrm {------})$ time derivative of resolved global helicity $H(t)$ (left-hand side of (6.4) or (6.5)) obtained by differentiating the time series of $H(t)$ collected at every Navier–Stokes time step.

Figure 20. Plane-averaged helicity as a function of $z$ (vortex propagation direction) and $t$. The white dashed lines show the positions of vortex rings estimated by the planed-averaged vorticity magnitude, with velocities $v_a, v_{b1}$ and $v_{b2}$ defined in figure 6 and calculated the same way as in figure 7. The DNS and LES results are compared for each Reynolds number.

Figure 21. Plane-averaged helicity production $\textrm {d} h^*/\textrm {d} t^*$ as a function of $z$ (vortex propagation direction) and $t$ from the DNS results. The white dashed lines show the positions of vortex rings estimated by the planed-averaged vorticity magnitude, with velocities $v_a, v_{b1}$ and $v_{b2}$ defined in figure 6 and calculated the same way as in figure 7.

At higher Reynolds numbers, reconnection and Kelvin-wave collisions tend to produce higher mode helical structures that are responsible for more intense local helicity dissipation. Kida & Takaoka (Reference Kida and Takaoka1988) analysed two reconnection processes, bridging and vortex cancellation from the view point of helicity evolution, claiming that bridging slows down the rate of total helicity change. On the contrary, in the case of a knotted vortex, the bridge structures (see § 5.2) accelerate the rate of change of helicity.

Any helicity variation is indeed a result of reconnection of vortex lines in the presence of viscous effects. During a strong reconnection event, as in the case of a knotted vortex, the vortex tubes are abruptly flattened by the strain field, forcing the vortex lines to reconnect rapidly, which results in rapid changes in the total helicity. An example of weak reconnection is viscous dissipation of twist helicity, as discussed by McGavin & Pontin (Reference McGavin and Pontin2019), which is due to slippage of the vortex lines within the tubes. Arguing why total helicity should change is more straightforward from a topological standpoint. The reconnection event experienced by the trefoil knot transforms the topology from a stable one to an unstable one; as can be seen from figure 18, helicity is conserved up until reconnection. As discussed by Scheeler et al. (Reference Scheeler, Kleckner, Proment, Kindlmann and Irvine2014), before reconnection the helicity is stored in the form of writhe helicity, which is a stable form of helicity storage. Whereas after reconnection, a large part of helicity is converted into twist helicity, which is more sensitive to viscous dissipation. The net negative change in total helicity, as observed in this flow, is due to the local positive helicity density near the reconnection spot being viscously dissipated more intensely than the negative helicity in other parts of the domain. This is shown by the disappearance of negative helicity patches for $Re_\varGamma =2000$ and $Re_\varGamma =6000$ DNS in figure 20. This ‘annihilation’ mechanism, observed in the current simulations, has also been suggested by Kimura & Moffatt (Reference Kimura and Moffatt2014) in their analysis of a pair of skewed Burgers vortex tubes. In the case of the present LES, the viscous dissipation modelled by the adopted SGS closure makes up for the unresolved dissipation, and, as a result, the LES solution predicts the helicity jump in both Re-case reasonably well. The study of the helicity dynamics associated with strong reconnection events such as this one would also require a careful analysis of the effects of the vortex-core size, as thicker cores prolong the reconnection time, distributing the enstrophy production events, as well as the total helicity change, more gradually in time.

6.2. Centreline and vortex-core-integrated helicity

Scheeler et al. (Reference Scheeler, Kleckner, Proment, Kindlmann and Irvine2014) reported that the centreline helicity, defined as in this manuscript by (6.2), was conserved for the knotted vortex problem, even in the presence of viscous effects. On the other hand, Kimura & Moffatt (Reference Kimura and Moffatt2014), when modelling the same flow as a crossed vortex pair problem, have shown that the helicity integrated within core volumes as in (6.3) decays right after reconnection. In the present study both quantities are examined in figures 22(a) and 22(c) for both Reynolds numbers. Both methods indicate that after the reconnection, the large trailing ring carries the majority of the helicity. However, the tube-integrated helicity agrees well with the global helicity, by showing the helicity production during the reconnection. The centreline helicity, on the other hand, shows discrepancies compared with the reference total helicity results after the reconnection, and appears to not change significantly before and after reconnection.

Figure 22. Evolution of centreline helicity, tube-integrated helicity and circulation, for the entire vortex structure (black), and small leading ring (yellow) and large trailing ring (red) after reconnection. (a,c) Time series of vortex-tube-integrated helicity density $H_V(t)$ $(\square )$, and centreline helicity $H_C(t)$ $(\textrm {------})$ normalized by initial circulation from DNS results at $Re_{\varGamma }=2000$ and $Re_{\varGamma }=6000$, respectively. Green and blue triangles $(\triangleleft )$ show the normalized total helicity $H(t)$ for $Re_{\varGamma }=2000$ and $Re_{\varGamma }=6000$, respectively. (b,d) Evolution of instantaneous circulation $\varGamma _n(t)$ (see (6.8)) from DNS results at $Re_{\varGamma }=2000$ and $Re_{\varGamma }=6000$, respectively, shown by circles $(\circ )$. The grey box is masking the reconnection time when the vortex centreline is not clearly defined.

The helicity obtained by integration along the vortex line via both methods exhibits oscillations after reconnection, which can be explained by colliding Kelvin waves and the simultaneous formation of annular structures.

It is noteworthy that the centreline helicity in the present study does not vary significantly during reconnection. One possible reason for this is that this helicity model takes the axial velocity value from the centreline, and does not take into account the distribution of the axial flow $\tilde {w}$, defined in (5.3), within the vortex core. Figure 23 shows that, caused by the vortex bursting during the reconnection and collisions of Kelvin waves, the axial flow distribution on different cross-sections normal to the vortex tube can be irregular and far from being uniform, which is implicitly assumed by the centreline helicity method (6.2). This explains why in figure 22 the discrepancy between the total helicity $H$ and the centreline helicity $H_c$ is noticeably larger in the $Re_\varGamma =6000$ case, where the axial flow distribution on the cross-sections is more inhomogeneous. Since the radial distribution of the axial flow can contribute significantly to the local helicity, the vortex-tube-integrated helicity $H_V$ should be considered as a more accurate representation of the helicity associated with a vortex tube than the mere centreline helicity $H_C$. This explains what is observed in figure 22, that is before reconnection, when the vortex tube is practically unperturbed, centreline helicity shows a good agreement with global helicity and tube-integrated helicity, whereas after reconnection, there is a visible discrepancy among the three quantities.

Figure 23. Axial flow distribution ($\hat {w}$) on three different cut planes on the large vortex ring after reconnection, based on the $Q$-isosurface of $Q\bar {R}^4/\varGamma ^2=323$ at $t^*=3.97$ and $Re_\varGamma =6000$ from DNS results.

The instantaneous measured circulation $\varGamma _n(t)$ is shown in figures 22(b) and 22(d). The measured circulation $\varGamma _n$ is averaged along the vortex length after obtaining a local value on every normal slice according to

(6.8)\begin{equation} \varGamma_n(t)=\frac{1}{N_{{cut}}} \left(\sum_{i}^{N_{{cut}}} \iint_{A_n(\tilde{z}_i)} \boldsymbol{{\omega}}\cdot\boldsymbol{{n}} \,\mathrm{d}A_n \right), \end{equation}

where $A_n(\tilde {z}_i)$ is the area of a circle lying in a cut plane perpendicular to the vortex tube at a given axial position $\tilde {z}_i$. It is centred at the vortex centreline, and its normal vector aligns with the tangential vector of the local centreline. The radius of the thus obtained sampling circle varies in time and it is heuristically set to $r_s = 1.333 r_e (t)$ to collect all the vortex lines within and around the vortex tube. Changing $r_s$ in the range $1.1 r_e (t)$ to $1.5 r_e (t)$ only introduces a variation of about 5 % in $\varGamma _n$. For all of the vortex loops considered, $N_{{cut}}=100$ evenly spaced cut planes are used. Note the measured circulation $\varGamma _n$ is different from the initial nominal circulation $\varGamma$, because this method does not ensure that the circulation is strictly calculated along a closed material line, as required by Kelvin's theorem.

Nonetheless, the large trailing ring (red symbols and lines) formed after reconnection does preserve the circulation of the initial trefoil reasonably well. A slow decay of the circulation is observed after reconnection, and it is more pronounced at the lower Reynolds numbers. This can be explained by the fact that at lower Reynolds numbers, the thread structures tend to interact and merge with the larger trailing ring, while for higher Reynolds numbers, a larger portion of thread structures are left suspended in between the small and large rings, and are eventually left to dissipate, as seen in the strip of high vorticity region between the small and large ring between $t\varGamma /\bar {R}^2=3.5$ and $t\varGamma /\bar {R}^2=4.0$ in figure 7(b). When there is an interaction between the thread structures and a large coherent ring, the vortex lines in the threads will pass through the normal cut plane $\mathrm {d} A_n$ irregularly; also, such interaction will distort the vortex lines in the large vortex tube. All of these will add to the disturbance in the measured circulation $\varGamma _n$.

The smaller ring (yellow symbols and lines) after reconnection has a lower measured circulation $\varGamma _n(t)$ than the trefoil and the larger trailing ring, and is also characterized by more unsteady behaviour at both Reynolds numbers. This is because the small ring is surrounded by intense secondary vortical activity, as made evident in figure 23.