1. Introduction
Vortex reconnection is described by Yao & Hussain (Reference Yao and Hussain2020b ) as a ‘fundamental topology-transforming dynamical event’. Studying its mechanism is important for predicting the behaviour of trailing vortices of an aircraft, understanding the turbulence cascade, and more importantly, the occurrence of a finite-time singularity (FTS) in Euler or Navier–Stokes equations (Yao & Hussain Reference Yao and Hussain2022). The starting configuration ranges from simple vortex rings and vortex tubes to complex knotted and linked vortices such as the Hopf link and the trefoil knot. Numerical simulations are usually carried out with either direct numerical simulations (DNS) or simplified models.
With a Biot–Savart (B–S) model, Kimura & Moffatt (Reference Kimura and Moffatt2017) (see also references therein) found that anti-parallel vortex filaments stretch out as they approach closely, forming a tent-like (or pyramid-like) structure before reconnection regardless of the initial configuration, thereby suggesting a universal route. To desingularise the B–S integral, a ‘cut-off’ parameter is added to the denominator. This implies that the integral will produce spurious results for length scales smaller than the cut-off, such as at the time of reconnection. Moffatt & Kimura (Reference Moffatt and Kimura2019a
) argued that the evolution at reconnection time depends only on the curvature, the core radius and the separation distance. They used the B–S law to obtain analytical expressions for the rate of change of these variables, resulting in a nonlinear dynamical system, hereafter referred as the Moffatt–Kimura (MK) model. With this and subsequent work (Moffatt & Kimura Reference Moffatt and Kimura2019b
), they suggested the possible occurrence of ‘physical’ singularity for both inviscid (Euler) flow and viscous flows when the vortex Reynolds number is
$Re_{\varGamma } = \varGamma /\nu \geqslant 4000$
. Here,
$\varGamma$
is the circulation strength, and
$\nu$
is the kinematic viscosity. In particular, they note that the vorticity amplification, i.e. the ratio of maximum vorticity at some time
$t_c$
to the initial vorticity
$\omega _{{max}}/\omega _{0}$
, takes very large values with increase in
$Re_{\varGamma }$
. This is important as it was shown with the theorem of Beale, Kato & Majda (Reference Beale, Kato and Majda1984) that if a FTS occurs at a critical time
$t_c$
, then the maximum norm of vorticity becomes unbounded as a consequence of
$\int_{0}^{t_c} \|\boldsymbol{\omega}(\boldsymbol{x}, t)\|_{\infty} dt = \infty$
. Employing a similar set-up with
$Re_{\varGamma }$
up to 4000, Yao & Hussain (Reference Yao and Hussain2020a
) (hereafter YH) used DNS to show that the vorticity amplification is, in fact, much smaller than that reported by Moffatt & Kimura (Reference Moffatt and Kimura2019b
). In particular, YH attributed the suppression of vortex growth to flattening of the vortex cores as they approach closely, and the braking effect of the bridges, both of which were ignored in the MK model.
Early work with DNS at low
$Re_{\varGamma }$
, e.g. Kida, Takaoka & Hussain (Reference Kida, Takaoka and Hussain1991), was instrumental in establishing a clear understanding of the physical mechanism of reconnection. The process is usually divided into three phases. In the first phase, called inviscid advection, the vortex rings approach each other due to self-induction and mutual induction, and they undergo strong vortex stretching. The cores also flatten, resulting in a dipole structure. Next, during the bridging phase, the anti-parallel vortex lines annihilate each other, resulting in the formation of bridge- or hairpin-like structures in a direction orthogonal to the interaction. While the bridge structures rapidly recede, they remain connected by remnant threads from incomplete reconnection, resulting in the threading phase. High
$Re_{\varGamma } \leqslant 40\,000$
simulations with vortex tubes (Yao & Hussain Reference Yao and Hussain2020b
, Reference Yao and Hussain2022) have revealed more intricate details with increasing
$Re_{\varGamma }$
, such as the formation of the
$k^{-5/3}$
spectrum as a result of successive reconnections, and the core flattening becoming more pronounced. Since DNS may become prohibitively expensive for larger
$Re_{\varGamma }$
, it is useful to study the problem with simplified toy models that preserve some of the physical processes elucidated above while being computationally cheaper.
In this work, we use a projection of Navier–Stokes equations on a set of logarithmically spaced discrete Fourier modes to reach very large Reynolds numbers, at moderate numerical cost. Such a projection was invented by Campolina & Mailybaev (Reference Campolina and Mailybaev2021) and termed logarithmic-lattice (or log-lattice, for short). Although it may superficially look like a mere three-dimensional (3-D) generalisation of well-known shell models of turbulence (see the review by Biferale Reference Biferale2003), the projection on log-lattice is a mathematically well-defined procedure, allowing us to preserve the main symmetries and conservation laws of the original equations, without the need of adjustable parameters. The interest in this technique in the context of a FTS has already been demonstrated by Campolina & Mailybaev (Reference Campolina and Mailybaev2018), who showed that the Euler equation on log-lattice develops a blow-up at finite time. This blow-up is characterised by a chaotic attractor that spans a wide range of scales, out of reach of present DNS (see figure 2 of Campolina & Mailybaev Reference Campolina and Mailybaev2018). Log-lattice is also gaining considerable interest in geophysical applications where large Reynolds numbers are common and one can perform numerous simulations with a wide range of parameters to obtain general scaling laws that can later be turned into meaningful parametrisations for climate models (Pikeroen et al. Reference Pikeroen, Barral, Costa and Dubrulle2023). Both Pikeroen et al. (Reference Pikeroen, Barral, Costa and Dubrulle2023) and Barral & Dubrulle (Reference Barral and Dubrulle2023) note that the scaling laws obtained for homogeneous Rayleigh–Bénard convection appear to be mostly compatible with existing experimental and numerical work despite the decimation of modes. Furthermore, with log-lattice simulations of reversible Navier–Stokes equations, Costa, Barral & Dubrulle (Reference Costa, Barral and Dubrulle2023) were able to observe the existence of a second-order phase transition previously seen with DNS at a higher resolution by reaching scales smaller than usual DNS. The unique combination of using the equations of motion in their original form and the ability to span a wide range of scales with few modes makes log-lattice an attractive tool for a wide range of problems. The goal of the current work is to assess the similarity of log-lattice vortex reconnections to those seen with DNS and the B–S law by comparing with the results of MK and YH. In particular, the following questions are addressed.
-
(i) Since the initial conditions for a log-lattice simulation need to be defined in Fourier space, is it possible to construct vortex rings and set up the reconnection problem similar to Moffatt & Kimura (Reference Moffatt and Kimura2019a ) and Kida et al. (Reference Kida, Takaoka and Hussain1991)?
-
(ii) The reconstruction of physical space from lattice variables is currently an open question. Similar to Bohr et al. (Reference Bohr, Jensen, Paladin and Vulpiani1998) and Gürcan (Reference Gürcan2017), can one use the standard discrete Fourier transform to reconstruct the velocity field in physical space? If so, what are the possible issues with this approach?
-
(iii) Qualitatively, how does the log-lattice vortex reconnection compare with vortex reconnections seen with DNS and the B–S law? Which physical processes, if any, are preserved?
-
(iv) Finally, do log-lattice vortex reconnections show a blow-up of global quantities such as enstrophy
$\mathcal{{E}}$
and maximum norm of vorticity
$\|\boldsymbol{\omega }\|_{\infty }$
at finite
$Re_{\varGamma }$
up to
$10^8$
? How does the Euler simulation differ?
2. Numerical framework and initial conditions
2.1. Log-lattice framework
Only a brief description of the log-lattice framework is given here. Additional details can be found in Campolina & Mailybaev (Reference Campolina and Mailybaev2021). Starting with the Fourier transform of the incompressible Navier–Stokes equation,
\begin{align} \partial _t u_i + \textrm{i} k_j u_i * u_j &= -\textrm{i} k_i p - \nu k^2 u_i + f_i, \end{align}
\begin{align} \textrm{i} k_j u_j &= 0, \\[-6pt] \nonumber \end{align}
where
$\textrm{i} = \sqrt {-1}$
,
$k_i$
is the
$i$
th component of the
$d$
-dimensional wave vector
$\boldsymbol{k} = (k_1, \ldots , k_d)$
,
$p$
is the complex pressure field,
$f_i$
denotes the forcing, and
$\nu$
is the kinematic viscosity. In this work,
$f_i = 0$
for all simulations. When
$\nu = 0$
, the flow is inviscid, and the system reduces to the incompressible Euler equations. Here, (2.3) describes the main convolution operation that couples Fourier modes in triadic interactions such that
$k = q + r$
and
$q, r, k$
are any three nodes on a logarithmic lattice 
. The logarithmic lattice is the set
where
$\lambda \gt 1$
is the spacing factor. As shown in Campolina & Mailybaev (Reference Campolina and Mailybaev2021), non-trivial triad interactions exist only when the equation
$\lambda ^n = \lambda ^q + \lambda ^r$
has integer solutions for any
$(n, q, r) \in \mathbb{Z}^3$
. The following three families of solutions are known to exist, each with
$z$
interactions in
$d$
dimensions.
-
(i)
$\lambda = 2$
with
$z = 3^d$
. -
(ii)
$\lambda = \sigma \approx 1.325$
is the plastic number with
$z = 12^d$
. -
(iii) Any
$\lambda$
that satisfies
$1 = \lambda ^b - \lambda ^a$
where
$(a, b)$
are mutual prime integers and
$(a, b) \notin (1, 3), (4, 5)$
as they correspond to the plastic number. In this case,
$z = 6^d$
. For
$(a, b) = (1, 2)$
, the spacing factor
$\lambda = \phi \approx 1.618$
is the golden mean.
While decreasing the spacing factor
$\lambda$
can increase the node density and the number of interactions per node, which is desirable for simulating turbulent flows, this increases the computational cost. In this sense, simulations with
$\lambda = 2$
are computationally cheaper, but as pointed out in Barral & Dubrulle (Reference Barral and Dubrulle2023), this should be avoided for divergence-free flows because of the lack of a forward cascade. Furthermore, the more general absence of backscatter prohibits any inverse cascade as well. For instance, if one starts with an initial condition where
$\hat {u}(k) = 0$
for
$k \lt K$
, then this remains true for all subsequent time. This problem is not present for other values of
$\lambda$
.
In this work,
$\lambda = \phi$
is chosen for all simulations. The minimum wave vector
$k_{\textit{min}}$
is set to
$2\pi$
to match a simulation on a box of size
$L = 1$
. The initial grid size is set as
$(2N)^3 = 40^3$
. The Python code used in the simulations is adaptive, and new grid points are added or removed based on the fraction of energy contained in the outermost shells, i.e.
where
$E = ( {1}/{2}) (\boldsymbol{u},\boldsymbol{u})$
is the total energy,
$k_{\textit{min}}$
is the minimum wavenumber,
$\mathcal{M}$
is a spectral mask, and
$r_{\textit{tol}}$
is a very small threshold that is set to
$10^{-200}$
. When the above condition is satisfied, larger wavenumbers corresponding to smaller scales are added. The computational cost is reduced further by taking advantage of the grid symmetry along the initial axis
$f(-k) = \overline {f(k)}$
where
$\overline {({\cdot })}$
denotes complex conjugation. This means that the actual simulation is performed for
$N \times 2N \times 2N$
nodes instead of
$(2N)^3$
nodes for each velocity component.
2.2. Flow visualisation with log-lattice
To enable qualitative studies, 3-D velocity fields need to be reconstructed in physical space from the lattice variables. As suggested in Bohr et al. (Reference Bohr, Jensen, Paladin and Vulpiani1998) and Gürcan (Reference Gürcan2017), a simple algorithm would involve the discrete Fourier transform (DFT)
where
$\hat {\boldsymbol{u}}(\boldsymbol{k})$
is the velocity vector in Fourier space,
$\boldsymbol{k}$
is the non-uniformly spaced wave vector, and
$\boldsymbol{x}$
are evenly spaced sample points. If is an evenly spaced lattice, then (2.6) reduces to the regular DFT. While the complexity of the regular one-dimensional DFT is
$\mathcal{O}(N^2)$
, the complexity for log-lattices is
$\mathcal{O}(N\log _{\lambda }N)$
, similar to a fast Fourier transform, which reduces computation time. However, for a large number of sample points, this is still computationally expensive, therefore the Python DFT code was GPU-parallelised with the CuPy library (Okuta et al. Reference Okuta, Unno, Nishino, Hido and Loomis2017).
2.3. Initial conditions and simulation details
Two sets of initial conditions are chosen to study the oblique interaction of two vortex rings as shown in figure 1(a). The first set is similar to Moffatt & Kimura (Reference Moffatt and Kimura2019a
), with initially circular vortex rings with radii
$R$
, thickness
$\delta _0$
, and inclination angle
$\alpha = 45^{\circ }$
. They have equal and opposite circulation
$\varGamma$
, and are separated at their centrelines by a distance
$2s_0 = 0.4$
. The core size is
$\delta _0/R = 0.01$
. The only difference with the second set is the core size, which is
$\delta _0/R = 0.2$
, which is similar to that studied by Kida et al. (Reference Kida, Takaoka and Hussain1991). The larger core size enables better visualisation of the various reconnection processes.
Figure 1. (a) Schematic of the initial configuration. Figure adapted from YH. (b) Visualisation of vortex rings from log-lattices for the Kida-type ring with
$\delta _0/R = 0.2$
with the spacing factor
$\lambda = \phi$
. (c) As in (b) for
$\lambda \approx 1.237$
. The Q-criterion isosurfaces on grid size
$128^3$
are plotted at
$\text{Q} \gt 0.35\, \text{Q}_{{max}}$
and
$\text{Q} \gt 0.3\, \text{Q}_{{max}}$
, respectively.
Due to the spectral nature of log-lattices, the initial conditions need to be described in Fourier space as well. This immediately leads to the question of how to represent a vortex ring structure in Fourier space. An early idea exploited the fact that the Fourier transform of a 3-D Gaussian function is another 3-D Gaussian, and a rudimentary ring can be constructed by imposing the divergence-free condition. However, this type of ring has essentially one length scale to control its thickness, but the radius of the ring cannot be controlled. The current idea makes use of the Dirac delta function to represent a two-dimensional circle in 3-D space. For a vortex ring of radius
$R$
, the vorticity is given by
\begin{equation} \boldsymbol{\omega} (\boldsymbol{x}) = \varGamma \int \delta ^{(3)}(\boldsymbol{x} - \boldsymbol{R}(\phi ))\, \boldsymbol{T}(\phi )\,\textrm{d}\phi, \end{equation}
where
$\boldsymbol{R}(\phi ) = (R\cos \phi , R\sin \phi , 0)$
, and the unit tangent vector is
$\boldsymbol{T}(\phi ) = (-\sin \phi , \cos \phi , 0)$
. Taking the Fourier transform of (2.7) gives
\begin{equation} \hat {\boldsymbol{\omega }}(\boldsymbol{k}) = \varGamma \int \textrm{e}^{-\boldsymbol{i} \boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}(\phi )}\, \boldsymbol{T}(\phi )\,\textrm{d}\phi. \end{equation}
Using polar coordinates with
$x = R\cos \phi$
,
$y = R\sin \phi$
, and correspondingly
$k_x = k_{\perp }\cos \alpha$
,
$k_y = k_{\perp }\sin \alpha$
, such that
$R = \sqrt {x^2 + y^2}$
and
$k_{\perp } = \sqrt {k_x^2 + k_y^2}$
, the dot product
$\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{R}(\phi )$
reduces to
$Rk_{\perp } \cos (\phi - \alpha )$
, and (2.8) can be written as
\begin{equation} \hat{\boldsymbol{\omega }}(\boldsymbol{k}) = \varGamma \int \textrm{e}^{-\textit{i} Rk_{\perp}\cos (\phi - \alpha )} [{-}\sin \phi , \cos \phi , 0]\,\textrm{d}\phi. \end{equation}
Let
$\theta = \phi - \alpha$
such that
$\textrm{d}\theta = \textrm{d}\phi$
, and after some trigonometric manipulation, we note that the result
$\int \textrm{e}^{\textrm{i}Rk_{\perp }\cos \theta } \cos \theta\, \textrm{d}\theta$
can be conveniently expressed as a Bessel function
$-\textrm{i}\,J_{1}(k_{\perp }R)$
. Further, assuming
$\varGamma = 1$
and multiplying (corresponding to a convolution in real space) with a 3-D Gaussian of width
$\delta _0$
to represent the thickness of the vortex ring, the following closed-form expressions are obtained:
\begin{align} \hat {\omega }_{x}(\boldsymbol{k}) &= -\textrm{i}\,J_1(k_{\perp }R) \frac {k_y}{k_{\perp }} \textrm{e}^{\frac{-(\| \boldsymbol{k}\| \,\delta _0)^2}{2}}, \end{align}
\begin{align} \hat {\omega }_{y}(\boldsymbol{k}) &= \textrm{i}\,J_1(k_{\perp }R) \frac {k_x}{k_{\perp }} \textrm{e}^{\frac{-(\| \boldsymbol{k}\| \,\delta _0)^2}{2}}, \end{align}
\begin{align} \hat {\omega }_{z}(\boldsymbol{k}) &= 0. \end{align}
This describes a single vortex ring in Fourier space having radius
$R$
and thickness
$\delta _0$
oriented along the
$z$
direction. Straightforward extensions of the formula allow for translation and rotation of the ring. For instance, to translate the ring in the
$-z$
direction by a distance
$s_0$
, one can multiply (2.10)–(2.12) with
$e^{-\textrm{i}k_zs_0}$
(see (A1) for the complete equation). Multiple rings can be added by superposition. Finally, the initial velocity field can be obtained by applying the B–S law,
\begin{equation} \hat{\boldsymbol{u}}(\boldsymbol{k}) = \frac {\textrm{i}\boldsymbol{k} \times \hat{\boldsymbol{\omega}}(\boldsymbol{k})}{\|\boldsymbol{k}\|^2} \quad \text{for } \boldsymbol{k}\neq 0. \end{equation}
To study the qualitative aspects of reconnection, one needs to define a vortex. In this work, the second invariant of the velocity gradient tensor
$\boldsymbol{\nabla }\boldsymbol{u}$
, i.e. the Q-criterion (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988), is chosen, which identifies vortices as regions where rotation dominates over the strain:
\begin{equation} \text{Q} = \frac {1}{2} \big(\|\boldsymbol{\varOmega }\|^2 - \|\boldsymbol{S}\|^2\big) \gt 0, \end{equation}
where
$\boldsymbol{\varOmega } = ( {1}/{2}) [\boldsymbol{\nabla }\boldsymbol{u} - (\boldsymbol{\nabla }\boldsymbol{u})^{\textrm{T}}]$
is the spin tensor, and
$\boldsymbol{S} = ({1}/{2}) [\boldsymbol{\nabla }\boldsymbol{u} + (\boldsymbol{\nabla }\boldsymbol{u})^{\textrm{T}}]$
is the strain-rate tensor. This can be seen as an immediate improvement over vorticity magnitude
$\|\boldsymbol{\omega }\| = \|\boldsymbol{\nabla }\times \boldsymbol{u}\|$
, which is known to misidentify shearing motions as vortices as both regions possess non-zero vorticity (Lugt Reference Lugt1979). While numerous works in the vortex reconnection literature appear to favour the
$\lambda _2$
criterion, which identifies regions as vortices when the middle eigenvalue
$\lambda _2$
of the symmetric tensor
$(\boldsymbol{S}^2 + \boldsymbol{\varOmega }^2)$
possesses negative values
$\lambda _2 \lt 0$
, Chakraborty, Balachandar & Adrian (Reference Chakraborty, Balachandar and Adrian2005) showed that it is possible to obtain equivalent thresholds among popular methods such as the Q,
$\lambda _2$
and
$\varDelta$
criteria so that a particular choice among these criteria should not affect the qualitative results presented here.
The initial conditions for the Kida-type ring with
$\delta _{0}/R = 0.2$
, after applying the inverse DFT, are visualised with the Q-criterion in figure 1(b). Applying a threshold
$\text{Q} \gt 0.35\text{Q}_{{max}}$
reveals the vortex rings. The visualisation seems to suggest that the vortex rings are deformed, which if true, may invoke some unintended instability during the simulation. By comparing an additional log-lattice simulation of head-on collision of vortex rings with DNS, it is shown in Appendix A that the apparent deformation can be attributed to the approach taken for the real space representation, and the vortex rings are actually circular. One can also take advantage of the fact that large-scale modes are grouped closer on a logarithmic grid to ensure better physical space representation. This means that using rings with large radii (e.g.
$R = 0.1L$
, where
$L$
is the integral scale) and large thicknesses leads to a better result. Further improvements can be made at the expense of increased computational cost, e.g. including the modes containing zero components
$\boldsymbol{k}_i = 0$
, and reducing the grid spacing factor
$\lambda$
. All visualisations in this paper are shown with simulations performed with the
$\boldsymbol{k}_i = 0$
modes. The zero modes where
$\boldsymbol{k}_i = 0$
would couple to every other mode
$p$
such that
$\boldsymbol{k}_i = p_i - p_i$
in triadic interactions. These interactions are non-local, and they may enable, for instance, sweeping effects of small-scale structures by the large-scale structures. A qualitative comparison of the full-time evolution of the vortex rings with and without
$\boldsymbol{k}_i = 0$
modes is shown in the supplementary movie is available at https://doi.org/10.1017/jfm.2025.10900. The effect of reducing
$\lambda$
can be seen from figure 1(c), where the rings appear more circular than in figure 1(b). Due to the irregular spacing of the Fourier modes, using the inverse DFT approach as is will include artefacts/images with lower amplitude (similar to ringing artefacts in the signal processing literature) that are generally not visible at larger thresholds but become more apparent as the simulation progresses and turbulence decays. All visualisations and plots produced with real space data were made early times where the impact of these artefacts was minimal. To track the behaviour of global quantities such as enstrophy and maximum norm of vorticity over long times, they were calculated directly in Fourier space, which is unaffected by the artefacts.
3. Results and discussion
3.1. Core flattening, bridging and threading
Figure 2(a–d) show the temporal evolution of log-lattice vortex reconnection within the smaller box indicated in figure 1(b) for the Kida-type ring at
$Re_{\varGamma } = 10^4$
. The Q-criterion fields are thresholded at
$10\,\%$
of the maximum at each time step, and shaded with axial vorticity
$\omega _y$
to delineate the opposite rotation of the vortices. As explained in Kida et al. (Reference Kida, Takaoka and Hussain1991), the first phase of reconnection is inviscid advection, where the anti-parallel vortex structures approach each other due to curvature-driven self-induction, and collide. A similar process is evident in the log-lattice simulation as shown in figure 2(a,b), where the vortex rings are advected towards each other.
Figure 2. Q-criterion isosurfaces shaded with contours of axial vorticity at
$t = 0.01, 0.12, 0.2, 0.225$
. The shaded region with the colour map is the small box indicated in figure 1(b) thresholded at
$\text{Q} \gt 0.1\text{Q}_{{max}}$
. Other regions are uniformly shaded with
$30\,\%$
opacity and thresholded at
$\text{Q} \gt 0.08\text{Q}_{{max}}$
.
Next are the bridging and threading phases. Kida et al. (Reference Kida, Takaoka and Hussain1991) further explain that during the collision, the outermost vortex lines cancel each other at the point of contact due to viscous cross-diffusion until they form hairpin-like bridge structures in a direction orthogonal to the initial approach of the two vortices. Since the hairpin-bridge structures are strongly curved at the tip, this generates a large self-induced velocity that pushes the structures out of the plane of the paper, and backwards away from each other, which effectively stops the cancellation of vortex lines, and the remnant thread structures remain connected. At later times of the log-lattice simulation, i.e. figure 2(c,d) and 3(b), the regions shaded with axial vorticity clearly reveal the formation of both bridge and thread structures. This suggests that the three phases of vortex reconnection can be captured with a log-lattice simulation.
The entire process including core compression is better illustrated in figure 4(a–d), where contours of axial vorticity along the symmetric plane
$(x,z)$
are plotted at various times. The core becomes flatter upon approach, resulting in a head–tail dipole structure similar to that of Yao & Hussain (Reference Yao and Hussain2020a
). Interestingly, the vortex core flattens further and almost smears out with increasing
$Re_{\varGamma }$
, even when the core is thin, as evidenced in figure 4(e–h), where contours of axial vorticity are plotted for the MK-type ring at
$t = 0.115$
.
Figure 4. (a–d) Temporal evolution of the vortex core shape is shown with contours of axial vorticity with level sets
$\omega _y = [0.75{-}0.99]\,\omega _{y, {max}}$
in the symmetric
$(x,z)$
plane for the Kida-type ring at
$Re_{\varGamma } = 10^4$
. (e–h) Contours of axial vorticity at
$t = 0.115$
and level sets
$\omega _y = [0.85{-}0.99]\,\omega _{y, {max}}$
are shown in the symmetric
$(x,z)$
plane for the MK-type ring for increasing
$Re_{\varGamma }$
.
Reconnection is generally accompanied by the formation of vortex sheets (cf. figure 6 of Yao & Hussain Reference Yao and Hussain2022), which along with vortex tubes are considered to be the geometric building blocks of turbulence (Bermejo-Moreno, Pullin & Horiuti Reference Bermejo-Moreno, Pullin and Horiuti2009). Since
$Q\gt 0$
identifies regions where spin dominates strain, it is not clear from the visualisations shown in figure 2 if vortex sheets are formed in log-lattice simulations. Following Bermejo-Moreno et al. (Reference Bermejo-Moreno, Pullin and Horiuti2009) and Yao & Hussain (Reference Yao and Hussain2020b
), one can use the positive eigenvalues of symmetric second-order velocity gradient tensor
$A_{\textit{ij}} \equiv S_{ik} \varOmega _{kj} + S_{jk} \varOmega _{ki}$
proposed by Horiuti & Takagi (Reference Horiuti and Takagi2005) to educe vortex sheets. This takes advantage of the fact that both vorticity and strain rate in a vortex sheet are large and correlated. Examining closely at a large threshold of
$[A_{\textit{ij}}]_+ \gt 0.5[A_{\textit{ij}}]_{+, {max}}$
, it can be seen from figure 5(a) that these structures are at the edge of the Q-criterion vortex structure where the threads are formed. This is to be expected since it was already noted in Bermejo-Moreno et al. (Reference Bermejo-Moreno, Pullin and Horiuti2009) that the Q-criterion structures tend be surrounded, overlapped or intersected with structures of
$[A_{\textit{ij}}]_+$
. Interestingly, the vortex sheets seem to persist after the bridge structures are formed (see figure 5
b), and start thinning out and move away from the centre of the reconnection zone.
Figure 5. The dark green shaded regions are structures of
$[A_{\textit{ij}}]_+$
at times (a)
$t = 0.2$
, (b)
$t = 0.225$
, thresholded at
$[A_{\textit{ij}}]_+ \gt 0.5[A_{\textit{ij}}]_{+, {max}}$
. Other regions shaded uniformly at
$30\,\%$
opacity are Q-criterion structures at
$\text{Q} \gt 0.08\text{Q}_{{max}}$
.
At large
$Re_{\varGamma }$
, numerous small-scale structures are generated as a result of successive reconnections where threads of sufficient intensity undergo further reconnections, as explained with the cascade picture by Melander & Hussain (Reference Melander and Hussain1988). These are distinctly absent in the log-lattice visualisations shown in figures 2 and 3. While one can reach very large wavenumbers (small scales) with few modes on a log-lattice, the small scales are sparse and may not be adequate to visualise the associated structures. A smaller spacing factor could help in ascertaining more details. This effect can be clearly seen in figure 6, where a lowpass filter is applied to the velocity field in Fourier space, i.e.
\begin{equation} \hat {\boldsymbol{u}}^{(c)}(\boldsymbol{k}) = \begin{cases} \hat {\boldsymbol{u}}(\boldsymbol{k}) & \text{if}\ \|\boldsymbol{k}\|_\infty \lt k_h, \\ 0 & \text{otherwise}. \end{cases} \end{equation}
The wavenumbers are filtered up to
$k_h$
in the inertial range
$\|\boldsymbol{k}\|_\infty \lt 30$
for
$\lambda = \phi$
, and
$\|\boldsymbol{k}\|_\infty \lt 40$
for reduced
$\lambda \approx 1.237$
. Comparing with figure 3(b), one can no longer see the thread structures for
$\lambda = \phi$
, but they are visible in figure 6(b), where more modes are available in the inertial range. However, the inertial range is much larger for the reduced
$\lambda$
, and the structures appear to be different. While qualitatively the dynamics displays the usual reconnection between the rings (not shown) seen at
$\lambda = \phi$
, further simulations with lower
$\lambda$
are necessary before it can be conclusively stated that the reconnection dynamics is not different. Computing the energy spectrum, defined as
\begin{equation} E(\boldsymbol{k}) = \langle |\hat {\boldsymbol{u}}|^2 \rangle _{S_k}, \end{equation}
where the average
$\langle {\cdot }\rangle _{S_k}$
is taken over wave vectors in the shell
$S_k$
, its time evolution as shown in figure 7 clearly illustrates the flow of energy to finer scales (larger wavenumbers) as time progresses for the Kida-type ring with
$\lambda = \phi$
. Interestingly, there is a marked increase in the flow of energy at time
$t = 0.2$
when the threads were visualised (cf. figures 2(c) and 3(b)), suggesting that the reconnection triggers the proliferation of small-scale structures.
Figure 6. Q-criterion structures similar to figure 3(b) filtered up to
$k_h$
in the inertial range (a)
$\|\boldsymbol{k}\|_\infty \lt 30$
for
$\lambda = \phi$
, (b)
$\|\boldsymbol{k}\|_\infty \lt 40$
for
$\lambda \approx 1.237$
.
Figure 7. Time evolution of the energy spectrum for the Kida case at
$Re = 10^4$
with
$\boldsymbol{k}_i = 0$
modes and
$\lambda = \phi$
. The
$y$
direction shows wavenumbers up to
$k_{{max}}$
for this simulation. The contour lines represent level sets of the function
$\log E(|\boldsymbol{k}|, t)$
.
It has been suggested by MK that there is a counter-intuitive increase of curvature at the tipping points induced by strong vortex stretching, which may be responsible for the large growth of vorticity. This was examined by YH, who concluded that the radius of curvature decreases slower than the MK model, and subsequently increases once the bridges are formed. Prompted by one of the reviewers, who suggested that the limited resolution in the DNS of YH could have curbed the growth of curvature, this is examined for the Kida-type ring with the log-lattice simulation. Similar to YH, the curvature is calculated with vortex lines by integrating,
\begin{equation} \frac {\textrm{d}\boldsymbol{x}(s)}{\textrm{d}s} = \frac {\boldsymbol{\omega }}{|\boldsymbol{\omega }|}, \end{equation}
and the seed points are chosen as the peak of
$\omega _y$
along the symmetry
$(x, z)$
plane. The vortex lines at different times show an increase in curvature at the tip region as seen in figure 8(a). While calculating the temporal evolution of
$\kappa _0^2/\kappa ^2(t)$
, a smooth plot could not be obtained (possibly due to artefacts), and a plot was made with every fourth time step so that one could get a sense of the trend. As illustrated in figure 8(b), the radius of curvature decreases strongly, which is consistent with the results of MK and YH. However, unlike YH, the radius of curvature decreases very strongly and almost vanishes just after
$t = 0.2$
. Whether this is due to the availability of access to much smaller scales with log-lattices than that achievable with DNS, or the naive approach taken to real space representation hindering the analysis, requires further investigation. Subsequently,
$\kappa _0^2/\kappa ^2(t)$
does start to increase, similar to YH, which could be attributed to the bridges receding.
Figure 8. (a) The black and red lines show the vortex lines at
$t = 0, 0.167$
, respectively. (b) Temporal evolution of
$\kappa _0^2/\kappa ^2(t)$
for the Kida case at
$Re = 10^4$
and
$\lambda = \phi$
.
3.2. Circulation transfer and separation distance scaling
Key evidence of reconnection is the transfer of circulation from the symmetry
$(x, z)$
plane to the collision
$(x,y)$
plane. This is calculated for one half of each plane with
$\varGamma _y = \int \boldsymbol{\omega }\boldsymbol{\cdot}\boldsymbol{n}_s\, \textrm{d}S$
and
$\varGamma _z = \int \boldsymbol{\omega }\boldsymbol{\cdot}\boldsymbol{n}_c\, \textrm{d}S$
, where
$\boldsymbol{n}_s$
and
$\boldsymbol{n}_c$
are unit vectors normal to the symmetric and collision planes, plotted in figure 9(a) for both Kida- and MK-type rings at
$Re_{\varGamma } = 10^4$
. The inviscid advection phase is not evident as in YH due to the choice of the box (cf. figure 1
b) where the complete vortex core is not visible during the initial time (see the shaded region in figure 2
a). Therefore, the circulation along the symmetry plane initially increases as a larger portion of the vortex core comes into view. Later, a pattern similar to that in YH is visible where the circulation continuously drops along the symmetry plane while consequently increasing along the collision plane.
An important question in reconnection studies is the rate of approach and separation of the vortices. If
$\varGamma$
is the only relevant dimensional quantity, then it was shown with dimensional analysis that the separation distance
$s_0$
would scale as
$s_0 \sim (\varGamma t)^{1/2}$
. Indeed, for superfluid vortex reconnections with the Gross–Pitaevskii model, the approach and separation rates were found to follow the same scaling, suggesting that this scaling may be universal (Villois, Proment & Krstulovic Reference Villois, Proment and Krstulovic2017). However, Yao & Hussain (Reference Yao and Hussain2020b
) note that for viscous vortex reconnection, this scaling depends on the core size, with slender cores maintaining the local assumption required for the
$1/2$
scaling. To examine the temporal evolution of separation distance during approach, an approach similar to that of YH is taken where the centroid of
$\omega _y$
at
$0.75\omega _{y, {max}}$
is calculated for one half of the symmetry plane, and is taken to be the centre of the tube. As shown in figure 9(b), the separation distance before reconnection is found to scale as
$s_0 \sim (t_0 - t)^{0.95}$
for the Kida-type ring, far from the
$1/2$
scaling obtained from dimensional analysis. However, as found in Yao & Hussain (Reference Yao and Hussain2020b
), reducing the core size better maintains the local assumption, and for the MK-type ring,
$s_0$
scales as
$(t_0 - t)^{0.55}$
.
Figure 9. Temporal evolution of (a) circulation
$\varGamma _y, \varGamma _z$
normalised by the total circulation
$\varGamma = \varGamma _y + \varGamma _z$
along the symmetric
$(x, z)$
plane shaded black, and collision
$(x,y)$
plane shaded cyan, respectively. (b) The separation distance
$s_0$
is shown for both the Kida-type (solid line) and MK-type (dashed line) rings at
$Re_{\varGamma } = 10^4$
.
Figure 10. Time evolution of (a,d) total enstrophy
$\mathcal{E}$
, (b,e) maximum norm of vorticity
$\| \omega \| _{\infty }$
, and (c,f) kinetic energy spectrum at peak enstrophy for increasing Reynolds numbers, for (a–c) Kida-type rings and (d–f) MK-type rings. Here,
$Re_{\varGamma } = 5\times10^3, 10^4, 10^5, 10^6, 10^7, 10^8, \infty $
are shaded red, maroon, green, cyan, blue, violet and black, respectively.
3.3. Behaviour of enstrophy and maximum norm of vorticity
The temporal evolution of enstrophy and maximum norm of vorticity are plotted in figures 10(a,b) and 10(d,e) for the Kida- and MK-type rings, respectively. As explained in § 2.3, these quantities are calculated directly from lattice variables in Fourier space. This means that they are free from artefacts, allowing us to examine them for long times. It is clear that both enstrophy and vorticity, while increasing with
$Re_{\varGamma }$
, remain finite even for very large
$Re_{\varGamma }$
. Only the Euler simulation (indicated with a solid black line) shows a blow-up for both cases, albeit at different times. This difference is to be expected, as it is already shown by Pikeroen et al. (Reference Pikeroen, Barral, Costa, Campolina, Mailybaev and Dubrulle2024) that the blow-up times for log-lattice simulations are sensitive to initial conditions. Figure 10(c,f) show the energy spectrum at peak enstrophy, where one can observe a Kolmogorov-type
$k^{-5/3}$
spectrum. At low
$Re_{\varGamma }$
, the range of the
$-5/3$
slope is confined to the small wavenumbers, but this expands with increasing
$Re_{\varGamma }$
, suggesting a large number of small-scale structures being generated as a result of successive reconnections.
Figure 11. Time evolution of enstrophy for the Kida-type ring (a) with and without the
$\boldsymbol{k}_i = 0$
modes and
$\lambda = \phi$
, (b) with different
$\lambda$
along with the
$\boldsymbol{k}_i = 0$
modes.
Influence of the
$\boldsymbol{k}_{\boldsymbol{i}} \boldsymbol{=} \boldsymbol{0}$
modes and spacing factor
$\boldsymbol{\lambda}$
While some comments have already been directed towards the influence of the zero modes and using different spacing factors, here it is addressed with the behaviour of total enstrophy for the Kida-type ring. In figure 11(a), the total enstrophy over its initial value is examined for early times until
$t = 4$
. The simulation without
$\boldsymbol{k}_i = 0$
modes shows an enstrophy growth as seen in the DNS reconnection of Yao & Hussain (Reference Yao and Hussain2020b
). However, the addition of the
$\boldsymbol{k}_i = 0$
modes almost doubles this peak. As explained previously in § 2.3, this is to be expected since the zero modes would couple other modes in triadic interactions, thereby resulting in an increase of enstrophy. The addition of
$\boldsymbol{k}_i = 0$
becomes necessary when non-local interactions are important. In this paper, it is being used for visualisation purposes. Whether this is sufficient to include all non-local interactions is an important question. An early work by Plunian & Stepanov (Reference Plunian and Stepanov2007) introduced non-local interactions in the sabra shell model. This could be straightforwardly implemented for log-lattices, and the effect of non-local interactions can be tested by tweaking the non-locality strength.
Now the effect of changing
$\lambda$
is examined. Lowering
$\lambda$
from
$2$
(dyadic) to
$\phi \approx 1.618$
(golden mean) and
$\sigma \approx 1.325$
(plastic number) increases the enstrophy each time as the number of interactions per mode is doubled (cf. figure 11
b). However, as pointed out in Campolina & Mailybaev (Reference Campolina and Mailybaev2021) and § 2.1, one can choose other values of
$(a{,}b)$
such that the lattice is non-degenerate and the number of interactions are same as that of the golden mean. With
$(a{,}b) = (2, 5)$
, this corresponds to
$\lambda \approx 1.237$
(visualised previously in figures 1(c) and 6(b)). With this spacing factor, the peak of enstrophy is similar to the golden mean (cf. figure 11
b), but the peak is reached at an earlier time. These results suggest that the addition of
$\boldsymbol{k}_i = 0$
modes and changing the spacing factor may have an impact on the log-lattice reconnection dynamics, and their examination is left for future work.
4. Conclusion
With the log-lattice technique, numerical simulations of two inclined vortex rings with core sizes
$\delta _0/R = 0.2, 0.01$
are performed for increasing
$Re_{\varGamma }$
up to
$10^8$
, along with an inviscid (Euler) simulation. It is shown that the log-lattice vortex reconnections appear to retain key physical processes seen in DNS, including core flattening and the formation of hairpin-like bridge structures that suppress vorticity amplification. In line with the DNS results of YH, the peak of the maximum norm of vorticity increases with
$Re_{\varGamma }$
but remains finite even at
$Re_{\varGamma } = 10^8$
, and a blow-up is seen only for the inviscid case. Other qualitative results observed with DNS studies – such as the presence of vortex sheets, proliferation of small-scale structures (by examining the temporal evolution of energy), behaviour of curvature, circulation transfer and separation distance scaling – are also captured quite well by log-lattices, making them suitable toy models to study vortex reconnections and their links to an FTS at a much lower computational cost than DNS, as shown in table 1. Future studies can be directed towards testing the accuracy of the log-lattice vortex reconnections by establishing quantitative comparisons with their DNS counterparts.
Table 1. Initial and final grid sizes, and time taken, for some cases of the log-lattice simulation. Bold font indicates simulations performed with the
$\boldsymbol{k}_i = 0$
modes. The CPU time indicates the time taken for the convolution operation (2.3) at the last time step of the simulation. The simulation at
$Re_{\varGamma } = 10^4$
with
$\boldsymbol{k}_i = 0$
is examined in detail in § 3.1.
While it is possible to use the inverse DFT to visualise the reconnection process, this is limited to early times due to the presence of artefacts that are a result of using non-uniform, non-integer Fourier modes. Further improvements can be focused on interpolating the Fourier modes to linearly spaced integer values to not only suppress artefacts but also to employ the standard fast Fourier transform algorithms. Another direction is to hide the artefacts with a clever choice of sampling points. This is currently being explored and will be a part of a future publication. Research can also be directed towards understanding the reconnection dynamics at different spacing factors and also by introducing further non-local interactions with the method described in Plunian & Stepanov (Reference Plunian and Stepanov2007). It is also currently unknown if log-lattice simulations converge to DNS when the spacing factor
$\lambda$
tends to 1. The delta function can be used to study other initial conditions, such as vortex tubes, by replacing the equation of a circle with that of a line, and convolving with a 3-D Gaussian to give it some thickness. This is at the core of our current efforts to study the interaction of vortex tubes with increasing complexity – by varying the circulation strength and core sizes, and introducing axial flow in the vortex cores at very high
$Re_{\varGamma }$
. Such studies of asymmetrical reconnection could be useful in devising a method for directly identifying vortex reconnection in turbulent flow data, which is a major challenge.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2025.10900.
Acknowledgements
We thank G. Krstulovic for suggesting the Dirac delta function to represent vortex rings in Fourier space. We are also grateful the editor D. Lohse and the three anonymous reviewers whose constructive feedback and criticism have helped us to improve the paper.
Funding
This research has been funded though the Agence Nationale pour la Recherche, via the grants ANR TILT grant agreement no. ANR-20-CE30-0035 and ANR BANG grant agreement no. ANR-22-CE30-0025.
Declaration of interests
The authors report no conflict of interest.
Data and code availability statement
The driver scripts and data generated from log-lattice simulations can be requested directly from the authors. The log-lattice Python code used in the paper was developed by A. Barral as a part of his PhD work, and is freely available along with a detailed documentation. See Barral et al. (Reference Barral, Dubrulle, Costa, Pikeroen and Lopez2024) for more details.
Appendix A. On the circularity of the vortex ring with
$\boldsymbol{\lambda = \phi}$
The real-space visualisation of the initial condition in figure 1(b) seems to suggest that the vortex rings may not be actually circular when the spacing factor is
$\lambda = \phi$
. If true, then this deformation of the rings could invoke an unintended instability during their temporal evolution. By comparing with a well-known result from DNS, it is shown here that the rings are actually circular, and the apparent deformation can be linked to the naive approach taken for the real-space representation.
With DNS, Kang, Yun & Protas (Reference Kang, Yun and Protas2020) (hereafter KYP) conducted a systematic search for potentially singular behaviour of the 3-D Navier–Stokes equation by solving a family of partial differential equation optimisation problems where the initial conditions with a given enstrophy are sought in such a way that the enstrophy within a time window is maximised. This resulted in two sets of initial conditions, which are termed symmetric and asymmetric (cf. figure 3 of KYP), whose subsequent flow evolution always shows a growth of enstrophy. Our interest is in the symmetric initial condition as shown in figure 8(b) of KYP, where two vortex rings oriented along the diagonal of the computational domain can be seen. The subsequent evolution of total enstrophy and component-wise enstrophy corresponding to each coordinate direction (cf. figure 3(a) of KYP) clearly show that the total enstrophy is always equally distributed among its components.
To orient the vortex ring along the diagonal of the computational domain in Fourier space, the unit vector is set to
$\boldsymbol{n} = (\pm 1, \pm 1, \pm 1)$
(with positive values being used for one ring, and negative values for the other) in the equation
\begin{equation} \hat {\boldsymbol{\omega }}(\boldsymbol{k}) = \varGamma J_1(\|\boldsymbol{n} \times \boldsymbol{k}\|\,R) \frac {\boldsymbol{n} \times \boldsymbol{k}}{\|\boldsymbol{n} \times \boldsymbol{k}\|}\, \textrm{e}^{\frac{-(\|\boldsymbol{k}\|\,\delta _0)^2}{2}} \textrm{e}^{- \textrm{i} k_z s_0}, \end{equation}
and
$\varGamma = \pm 1$
. When the unit vector is
$\boldsymbol{n} = (0, 0, -1)$
,
$\boldsymbol{n} \times \boldsymbol{k} = (-k_y, k_x, 0)$
,
$\|\boldsymbol{n} \times \boldsymbol{k}\| = k_{\perp } = \sqrt {k_x^2 + k_y^2}$
and
$s_0 = 0$
, the above reduces to (2.10)–(2.12).
Figure 12. (a) Temporal evolution of total enstrophy
$\mathcal{E}$
(maroon thick solid line) and component-wise enstrophy corresponding to each coordinate direction
$\mathcal{E}_1, \mathcal{E}_2, \mathcal{E}_3$
(black solid, dotted and dashed lines, respectively). (b,c) Initial condition for the head-on collision of vortex rings oriented along the diagonal of the computational domain visualised at different angles. (d–f) Evolution of head-on collision at
$t = 0.127, 0.204, 0.324$
. All visualisations show Q-criterion isosurfaces thresholded at
$\text{Q} \gt 0.4\text{Q}_{{max}}$
.
With a Kida-type ring, i.e.
$\delta _0/R = 0.2$
and all other parameters being the same as described in § 2.3, the vortex rings oriented along the diagonal are initialised on log-lattice. Similar to oblique reconnection, the real-space representation with Q-criterion in figure 12(b) shows an apparent deformation of the ring. However, it can be easily checked that the total enstrophy remains equipartitioned among the three coordinate directions for the initial condition. This equipartitioning is immediately disturbed even if the vortex rings are tilted gently from this position. Running the simulation with
$Re_{\varGamma } = 10^4$
, the temporal evolution of total enstrophy and component-wise enstrophy are plotted in figure 12(a). Similar to the results of KYP, the total enstrophy remains equipartitioned among its components throughout the simulation. This suggests that the initialised vortex rings are more or less circular, even if this is not apparent from the visualisation.
Open access
