Velocity gradient analysis of a head-on vortex ring collision

Abstract

We simulate the head-on collision between vortex rings with circulation Reynolds numbers of 4000 using an adaptive, multiresolution solver based on the lattice Green's function. The simulation fidelity is established with integral metrics representing symmetries and discretization errors. Using the velocity gradient tensor and structural features of local streamlines, we characterize the evolution of the flow with a particular focus on its transition and turbulent decay. Transition is excited by the development of the elliptic instability, which grows during the mutual interaction of the rings as they expand radially at the collision plane. The development of antiparallel secondary vortex filaments along the circumference mediates the proliferation of small-scale turbulence. During turbulent decay, the partitioning of the velocity gradients approaches an equilibrium that is dominated by shearing and agrees well with previous results for forced isotropic turbulence. We also introduce new phase spaces for the velocity gradients that reflect the interplay between shearing and rigid rotation and highlight geometric features of local streamlines. In conjunction with our other analyses, these phase spaces suggest that, while the elliptic instability is the predominant mechanism driving the initial transition, its interplay with other mechanisms, e.g. the Crow instability, becomes more important during turbulent decay. Our analysis also suggests that the geometry-based phase space may be promising for identifying the effects of the elliptic instability and other mechanisms using the structure of local streamlines. Moving forward, characterizing the organization of these mechanisms within vortices and universal features of velocity gradients may aid in modelling turbulent flows.

1. Introduction

1.1. Vortex rings

Vortex rings are ubiquitous flow phenomena in both applied and theoretical settings, with applications including sound generation, transport, mixing and vortex interactions (Shariff & Leonard 1992). In geophysical settings, vortex rings can be used to model entrainment and dispersion in particle clouds (Bush, Thurber & Blanchette 2003). They play important roles in the initial jets of volcanic eruptions (Taddeucci et al. 2015) and the transport of contaminated sediments disposed of in open-water settings (Ruggaber 2000). In biomechanical settings, vortex rings have been observed in the motions of blood in the human heart (Arvidsson et al. 2016) and in the propulsive motion of oblate medusan jellyfish (Dabiri 2005). Remarkably, separated vortex rings augment dandelion seed dispersal by prolonging flight through drag enhancement (Cummins et al. 2018). In aerodynamic settings, vortex rings are responsible for the so-called vortex ring state, which negatively impacts lift in helicopters (Johnson 2005) and the performance of offshore wind turbines (Kyle, Lee & Früh 2020). In experimental and numerical settings, the formation and pinch-off of vortex rings are of particular interest in jet flows involving nozzles and orifices (Gharib, Rambod & Shariff 1998; Mohseni, Ran & Colonius 2001; Krueger & Gharib 2003; O'Farrell & Dabiri 2014; Limbourg & Nedić 2021).

Vortex rings are also associated with complex instabilities and dynamics that relate more generally to the sustenance of turbulence. Flow instabilities in vortex rings depend primarily on the core vorticity distribution, the circulation Reynolds number $( Re_\varGamma = \varGamma / \nu )$ and the slenderness ratio $( \delta = a / R )$ (Balakrishna, Mathew & Samanta 2020). Here, $\varGamma$ is the circulation, $\nu$ is the kinematic viscosity, $a$ is the core radius and $R$ is the ring radius. We focus on the evolution of thin-cored vortex rings with Gaussian core vorticity profiles, no swirl and centroids ($Z$) that propagate along the $z$-axis. In cylindrical coordinates ($r, \theta, z$), this initial vorticity profile is written as

(1.1)\begin{equation} \omega_\theta (r, z; t = 0) ={\pm} \frac{\varGamma_0}{{\rm \pi} a_0^2} \exp \left(-\frac{(z - Z_0)^2 + (r - R_0)^2}{a_0^2}\right)\!, \end{equation}

where subscripts $({\cdot })_0$ denote parameter values at $t = 0$ and the sign of $\omega _\theta$ dictates the propagation direction. Since Gaussian vortex rings only satisfy the governing equations with infinitesimal core thickness, they initially undergo a rapid period of equilibration in which vorticity is redistributed throughout the core (Shariff, Verzicco & Orlandi 1994; Archer, Thomas & Coleman 2008; Balakrishna et al. 2020). Following instability growth, transition is often marked by the development of secondary vorticity in a halo around the core vorticity (Dazin, Dupont & Stanislas 2006; Bergdorf, Koumoutsakos & Leonard 2007; Archer et al. 2008). During turbulent decay, the shedding of secondary vortex structures to the wake can result in a stepwise decay in circulation (Weigand & Gharib 1994; Bergdorf et al. 2007).

Stability analyses of thin vortex rings are often (classically) formulated in terms of asymptotic expansions in $\delta$ (Widnall, Bliss & Tsai 1974; Widnall & Tsai 1977; Fukumoto & Hattori 2005). Infinitesimally thin vortex rings ($\delta \to 0$) are neutrally stable (Shariff & Leonard 1992). For rings with finite thickness ($\delta > 0$), the curvature instability occurs at first order in $\delta$ and the elliptic instability occurs at second order in $\delta$. The curvature and elliptic instabilities occur at short wavelengths and arise due to parametric resonance between Kelvin waves with core azimuthal wavenumbers separated by one and two, respectively (Fukumoto & Hattori 2005; Hattori, Blanco-Rodríguez & Le Dizès 2019). The curvature instability is attributed to a dipole field produced by the vortex ring curvature (Fukumoto & Hattori 2005; Blanco-Rodríguez et al. 2015; Blanco-Rodríguez & Le Dizès 2017). By contrast, the elliptic instability is attributed to a quadrupole field generated by straining induced by the ring or some external source (Fukumoto & Hattori 2005; Blanco-Rodríguez et al. 2015; Blanco-Rodríguez & Le Dizès 2016).

This elliptic instability acts to break up elliptic streamlines and is key to the development of three-dimensional transitional and turbulent flows (Kerswell 2002). In the context of vortex rings (or, more generally, strained vortices), it is sometimes called the Moore–Saffman–Tsai–Widnall instability (Fukumoto & Hattori 2005; Chang & Llewellyn Smith 2021) based on the initial investigations of Moore & Saffman (1975) and Tsai & Widnall (1976). The elliptic instability dominates the curvature instability for thin Gaussian vortex rings without swirl. However, the curvature instability becomes increasingly important for vortex rings with increasing $Re_\varGamma$ and decreasing $\delta$, as well as in vortex rings with swirl (Blanco-Rodríguez & Le Dizès 2017; Hattori et al. 2019).

While interesting in their own right, thin vortex rings often form canonical building blocks of more complex turbulent flows. Modified vortex geometries, such as elliptic vortex rings (Cheng, Lou & Lim 2016, 2019) and trefoil knots (Yao, Yang & Hussain 2021; Zhao et al. 2021), provide alternative means of probing vortex dynamics and interactions. Collisions between vortex rings and other vortex rings, walls and free surfaces are also commonly studied to investigate mechanisms underlying the turbulent cascade and the generation of small scales (see Mishra, Pumir & Ostilla-Mónico (2021) for a review). These mechanisms can be characterized using a variety of collision geometries, including head-on collisions (Cheng, Lou & Lim 2018; McKeown et al. 2018, 2020; Mishra et al. 2021), inclined collisions (Kida, Takaoka & Hussain 1991; Yao & Hussain 2020a,c) and axis-offset collisions (Zawadzki & Aref 1991; Smith & Wei 1994; Nguyen et al. 2021), among others. Boundary layers play an important role in vortex–wall interactions (e.g. by causing rebounding events) (Walker et al. 1987) and interactions with free surfaces can often be understood in terms of image vortices (Archer, Thomas & Coleman 2010).

Here, we focus on head-on collisions between identical vortex rings of opposite circulation, which have been noted for their rapid enstrophy production (Lu & Doering 2008; Ayala & Protas 2017; Kang, Yun & Protas 2020). They have been classically studied in the contexts of the formation of smaller rings through vortex reconnection and the formation of turbulent clouds at high $Re_\varGamma$ (Oshima 1978; Lim & Nickels 1992; Chu et al. 1995). Many recent investigations have focused particularly on the mechanisms (e.g. instabilities) underlying these transitional and turbulent processes (McKeown et al. 2018, 2020; Mishra et al. 2021).

For the head-on vortex ring collisions under consideration, the elliptic instability competes and interacts with the longer-wavelength Crow instability. The Crow instability (Crow 1970) is associated with the mutual interaction of perturbed counter-rotating vortices, which, in the linear regime, locally displaces the vortices without modifying their core structures (Leweke, Le Dizès & Williamson 2016). Mishra et al. (2021) provides a focused review of vortex ring collisions in the context of the these instabilities. For collisions at relatively low Reynolds numbers, the Crow instability can lead to the pinch-off of secondary vortex rings via local reconnections. At higher Reynolds numbers, the elliptic instability favours rapid disintegration of the vortex rings into a turbulent cloud.

McKeown et al. (2020) proposed that iterative elliptic instabilities between successive generations of antiparallel vortices can mediate the turbulent cascade in head-on vortex ring collisions. Mishra et al. (2021) also observed that the elliptic instability tends to dominate at high $Re_\varGamma$, although this behaviour is also sensitive to the slenderness ratio and vorticity distribution. In a different configuration involving symmetrically perturbed antiparallel vortices, Yao & Hussain (2020b) attributed the turbulent cascade at high $Re_\varGamma$ to an avalanche of successive vortex reconnections. In general, Ostilla-Mónico et al. (2021) found that collisions between counter-rotating vortices are indeed highly sensitive to the geometry of their configuration. They particularly found that the mechanisms mediating the cascade bear resemblance to the reconnection scenario (Yao & Hussain 2020b) when the vortices are nearly perpendicular, whereas they are more reminiscent of the iterative elliptic instability scenario (McKeown et al. 2020) when the vortices are more acutely aligned. These recent works share two common themes, (i) that the mode of transition and the formation of a cascade are sensitive to the details of the initial flow configuration and (ii) that the interplay between relevant instabilities is simultaneously important to the flow physics and difficult to capture.

1.2. Velocity gradients and vortices

The elliptic instability, which typically dominates head-on collisions between the vortex rings of interest at high $Re_\varGamma$ (McKeown et al. 2020; Mishra et al. 2021), is associated with elliptic streamlines (Kerswell 2002). This generic feature of strained vortical flows can be used to characterize the elliptic instability, which is typically difficult to discern in the complex interactions of multiscale vortices (Mishra et al. 2021; Ostilla-Mónico et al. 2021). Given the inherent complexity of turbulent flows, the geometry of local streamlines provides a relatively simple and interpretable means for characterizing flow features (e.g. vortices).

The instantaneous trajectory of a materially advecting fluid particle follows the streamlines, which are frame dependent. At a critical point, e.g. in a frame advecting with the particle, the velocity gradient tensor (VGT), $\boldsymbol{\mathsf{A}} = \boldsymbol {\nabla } \boldsymbol {u}$, determines, to linear order, the local structure of streamlines (Perry & Fairlie 1975; Perry & Chong 1987; Chong, Perry & Cantwell 1990). The scale-invariant shape of local streamlines is captured by normalizing the VGT as $\tilde{\!\boldsymbol{\mathsf{A}}} = \boldsymbol{\mathsf{A}}/A$ (Girimaji & Speziale 1995; Das & Girimaji 2019), where $A = \left\Vert \boldsymbol{\mathsf{A}} \right\Vert_{\rm F} = {\rm tr}( \boldsymbol{\mathsf{A}}^{\rm T} \boldsymbol{\mathsf{A}})^{1/2}$ is the Frobenius norm of the VGT, $(\cdot)^{\rm T}$ represents the transpose and, unless otherwise stated, non-bold versions of bold tensor quantities represent their Frobenius norms. This normalized VGT has been used to investigate the scalings, forcings and non-local features of the VGT dynamics (Das & Girimaji 2019, 2020a, 2022) and a similar analysis of vorticity gradients has been used to classify the geometry of local vortex lines (Sharma, Das & Girimaji 2021).

The principal invariants of $\tilde {\!\boldsymbol{\mathsf{A}}}$ instantaneously characterize local streamline topologies and geometries (Chong et al. 1990; Das & Girimaji 2019, 2020a,b). They are given by

(1.2a–c)\begin{equation} p_{\rule{0pt}{1.6ex} \mkern-3mu A} =-{\rm tr}(\tilde{\!\boldsymbol{\mathsf{A}}}),\quad q_{\rule{0pt}{1.6ex} \mkern-3mu A} = \tfrac{1}{2}({\rm tr}(\tilde{\!\boldsymbol{\mathsf{A}}})^{2} - {\rm tr}(\tilde{\!\boldsymbol{\mathsf{A}}}^2)),\quad r_{\rule{0pt}{1.6ex} \mkern-4mu A} =-{\rm det}(\tilde{\!\boldsymbol{\mathsf{A}}}), \end{equation}

where ${\rm tr}({\cdot })$ and ${\rm det}( {\cdot } )$ represent the trace and determinant, respectively. For incompressible flows ($p_{\rule {0pt}{1.6ex} \mkern -3mu A} = 0$), four classes of local streamline topologies are separated by degenerate geometries in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane. Using the invariants of $\tilde {\!\boldsymbol{\mathsf{A}}}$ is advantageous compared with using the invariants of $\boldsymbol{\mathsf{A}}$ since the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane is a bounded phase space and it provides a more complete representation of streamline geometries (Das & Girimaji 2019, 2020a,b). For example, the aspect ratio of purely elliptic local streamlines ($q_{\rule {0pt}{1.6ex} \mkern -3mu A} > 0$ and $r_{\rule {0pt}{1.6ex} \mkern -4mu A} = 0$) is completely characterized by $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$, but not by $Q = A^2 q_{\rule {0pt}{1.6ex} \mkern -3mu A}$. However, while the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane efficiently characterizes local streamline geometries at critical points, additional parameters are required to fully describe all geometries (Das & Girimaji 2020a).

Following Das & Girimaji (2020b), we consider the local streamline geometry in the context of the modes of deformation of a fluid parcel: extensional straining, (symmetric and antisymmetric) shearing and rigid rotation. The well-known Cauchy–Stokes decomposition of the VGT, $\tilde {\!\boldsymbol{\mathsf{A}}} = \tilde {\!\boldsymbol{\mathsf{S}}} + \tilde {\!\boldsymbol{\mathsf{W}}}$, disambiguates contributions from the symmetric strain rate tensor, $\tilde {\!\boldsymbol{\mathsf{S}}} = ( \tilde {\!\boldsymbol{\mathsf{A}}} + \tilde {\!\boldsymbol{\mathsf{A}}}^{\rm T})/2$, and the antisymmetric vorticity tensor, $\tilde {\!\boldsymbol{\mathsf{W}}} = ( \tilde {\!\boldsymbol{\mathsf{A}}} - \tilde {\!\boldsymbol{\mathsf{A}}}^{\rm T})/2$. It has enabled insightful characterizations of the VGT dynamics from the perspective of the strain rate eigenframe (Tom, Carbone & Bragg 2021). However, it does not disambiguate symmetric shearing from extensional straining in $\tilde {\!\boldsymbol{\mathsf{S}}}$ or antisymmetric shearing from rigid rotation in $\tilde {\!\boldsymbol{\mathsf{W}}}$. This limitation motivated the development of the triple decomposition of the VGT (Kolář 2007), which disambiguates all three fundamental modes of deformation.

Kolář (2004, 2007) originally formulated the triple decomposition of the VGT by identifying a ‘basic’ reference frame in which motions associated with elongation, rigid rotation and pure shearing can be isolated. However, identifying a basic reference frame requires a challenging pointwise optimization problem, the solution of which is typically approximated over a finite number of frames (Kolář 2007; Nagata et al. 2020). More recently, Gao & Liu (2018, 2019) introduced a unique triple decomposition, based on a related vorticity tensor decomposition (Liu et al. 2018; Gao et al. 2019), that is more computationally practical than that of Kolář (2004, 2007). This triple decomposition is formally performed in a local ‘principal’ coordinate system ($x^*, y^*, z^*$), which is related to the global coordinate system ($x, y, z$) by an orthogonal transformation. In this principal frame, denoted by $( {\cdot } )^*$, the triple decomposition is given in normalized form by

(1.3)\begin{equation} \tilde{\!\boldsymbol{\mathsf{A}}}^* = \underbrace{\begin{bmatrix} \dot{\epsilon}_{x^*} & 0 & 0 \\ 0 & \dot{\epsilon}_{y^*} & 0 \\ 0 & 0 & \dot{\epsilon}_{z^*} \end{bmatrix}}_{\boldsymbol{\dot{\epsilon}}^*} + \underbrace{\begin{bmatrix} 0 & 0 & 0 \\ \dot{\gamma}_{z^*} & 0 & 0 \\ \dot{\gamma}_{y^*} & \dot{\gamma}_{x^*} & 0 \end{bmatrix}}_{\boldsymbol{\dot{\gamma}}^*} + \underbrace{\begin{bmatrix} 0 & -\dot{\varphi_{z^*}} & 0 \\ \dot{\varphi_{z^*}} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}}_{\boldsymbol{\dot{\varphi}}^*}. \end{equation}

Here, $\boldsymbol {\dot {\epsilon }}^*$, $\boldsymbol {\dot {\gamma }}^*$ and $\boldsymbol {\dot {\varphi }}^*$ represent the normal straining, pure shearing and rigid body rotation tensors, respectively. Their constituents can be directly identified from the components of the VGT in the principal frame (Gao & Liu 2018, 2019; Das & Girimaji 2020b). Their representations in the global coordinates ($\boldsymbol {\dot {\epsilon }}$, $\boldsymbol {\dot {\gamma }}$ and $\boldsymbol {\dot {\varphi }}$) can subsequently be recovered by inverting (i.e. transposing) the original orthogonal transformation (Gao et al. 2019).

The components of the normal straining tensor represent the real parts of the eigenvalues of $\tilde {\!\boldsymbol{\mathsf{A}}}$, which are identical to those of $\tilde {\!\boldsymbol{\mathsf{A}}}^*$. For points with rotational local streamlines, $\tilde {\!\boldsymbol{\mathsf{A}}}$ has a pair of complex eigenvalues and the real eigenvector defines the local rotation axis. In this case, the transformation to the principal frame is identified by (i) using a real Schur decomposition to align the $z^*$-axis with the real eigenvector of the VGT and (ii) orienting the $x^*$–$y^*$ plane to minimize the local rotational speed (Liu et al. 2018; Das & Girimaji 2020b). One advantage of (1.3) is that it provides representations of the strength ($2 \dot {\varphi }_{z^*}$) and the axis ($z^*$) of rigid rotation that are Galilean invariant (Wang, Gao & Liu 2018). Unlike the rotational case, the VGT has only real eigenvalues when the local streamline geometry is non-rotational ($\boldsymbol {\dot {\varphi }}^* = \boldsymbol {0}$). In this case, the principal frame is identified by using a Schur decomposition to transform the VGT into a triangular tensor. The modes of deformation are then isolated by decomposing this transformed tensor into a normal, diagonal tensor representing normal straining and a non-normal, strictly triangular tensor representing pure shearing (Keylock 2018; Das & Girimaji 2020b).

The triple decomposition enables refined analyses of the influences of fundamental constituents of the VGT. For example, the original triple decomposition (Kolář 2007) has been used to show that lifetimes of fundamental flow structures at macroscopic scales (where viscosity can be neglected) can be related to stability of rigid rotation, linear instability of pure shearing and exponential instability of irrotational straining (Hoffman 2021). At small scales, the more recent triple decomposition (Gao & Liu 2018, 2019) has been used to show that pure shearing is typically the dominant contributor to energy dissipation (Wu et al. 2020) and intermittency (Das & Girimaji 2020b) in turbulent flows. Further, the symmetric and antisymmetric components of $\boldsymbol {\dot {\gamma }}^*$ are given by $\boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu S} = ( \boldsymbol {\dot {\gamma }}^* + \boldsymbol {\dot {\gamma }}^*{}^{\rm T} )/2$ and $\boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu W} = ( \boldsymbol {\dot {\gamma }}^* - \boldsymbol {\dot {\gamma }}^*{}^{\rm T} )/2$, respectively. In this manner, the triple decomposition is more refined than the Cauchy–Stokes decomposition since  $\tilde {\!\boldsymbol{\mathsf{S}}}^* = \boldsymbol {\dot {\epsilon }}^* + \boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu S}$  and  ${\tilde {\!\boldsymbol{\mathsf{W}}}^* = \boldsymbol {\dot {\varphi }}^* + \boldsymbol {\dot {\gamma }}^*_{\rule {0pt}{-2ex} \mkern -2mu W}}$ (Gao & Liu 2019; Das & Girimaji 2020b). As described in detail by Das & Girimaji (2020b), the triple decomposition also enables a natural characterization of local streamline topologies and geometries. Similar topological analyses of vortical flow features have also been proposed (Nakayama 2017), but we focus on the triple decomposition for the advantages outlined herein.

The ability of the triple decomposition to capture local streamline structure in terms of fundamental modes of deformation has guided efforts to define improved vortex criteria. There are an abundance of criteria to identify vortices that are based on various features (e.g. eigenvalues) of the VGT and that adopt various philosophies of what constitutes a vortex (Chakraborty, Balachandar & Adrian 2005; Epps 2017; Günther & Theisel 2018; Liu et al. 2019a; Haller 2021). Debates surrounding these criteria primarily involve their (i) philosophical underpinnings, (ii) threshold sensitivities and (iii) observational invariances.

Regarding (i), the Cauchy–Stokes decomposition underlies many common symmetry- based vortex criteria, including the $Q$ (Hunt, Wray & Moin 1988) and $\lambda _2$ (Jeong & Hussain 1995) criteria. Local streamline topology underlies many common geometry-based vortex criteria, including the $\varDelta$ (Chong et al. 1990) and $\lambda _{ci}$ (Chakraborty et al. 2005) criteria. Like the geometry-based methods, and unlike the symmetry-based methods, the rigid vorticity criterion ($\dot {\varphi }_{z^*} > 0$) (Tian et al. 2018) captures all rotational local streamline geometries under the assumption that rigid rotation is an essential ingredient of a vortex (Liu et al. 2019a; Das & Girimaji 2020b). The philosophical distinction between symmetry-based and geometry-based criteria also underlies the so-called ‘disappearing vortex problem’ in which, fixing the VGT configuration and strain rate, increasing only the vorticity magnitude can remove a geometry-based vortex from the flow (Chakraborty et al. 2005; Kolář & Šístek 2020, 2022). However, we here adopt the geometry-based viewpoint since, unlike vorticity, rigid rotation persistently underlies rotational local streamline topologies in all inertial frames. This interpretation in terms of local streamline topology has the potential to elucidate connections to related (e.g. elliptic) instabilities.

Regarding (ii), the Omega ($\varOmega$) class of vortex criteria (Liu et al. 2016; Dong et al. 2018; Dong, Gao & Liu 2019; Liu & Liu 2019) is advantageous since it uses quantities that are bounded and less threshold sensitive than the aforementioned methods. Regarding (iii), whereas most common vortex criteria are Galilean invariant, they are typically not objective since they are not preserved in rotating reference frames (Epps 2017; Günther & Theisel 2018). However, the objectivized $\dot {\varphi }_{z^*}$ (Liu, Gao & Liu 2019b) and objectivized $\varOmega$ (Liu et al. 2019c) criteria, which are formulated by replacing $\tilde {\!\boldsymbol{\mathsf{W}}}$ with its deviation from its global spatial mean, remain invariant in these reference frames. Moreover, they are among the only compatible (i.e. self-consistent) objectivized vortex criteria out of the modifications commonly associated with the vortex criteria we have discussed (Haller 2021). This advantage enhances the experimental verifiability and clarifies the physical significance of visualizations of the corresponding vortex structures.

Synthesizing the advantages of the geometry-based vortex definitions and the $\varOmega$ class of vortex criteria, we identify vortices using the $\varOmega _r$ method in the present investigation. This criterion is formulated in terms of the quantity

(1.4)\begin{equation} \varOmega_r = \frac{(\boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{z^*})^2}{2 (\boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{z^*})^2 - 4 \lambda_{ci}{\vphantom{\left( \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{z^*} \right)}}^{\mkern-14mu 2} + 4\varepsilon_{vort}}, \end{equation}

where $\lambda _{ci}$ is the imaginary part of the complex eigenvalues of $\boldsymbol{\mathsf{A}}$, $\boldsymbol {e}_{z^*}$ is the unit vector along the $z^*$-axis and $\varepsilon _{vort}$ is a numerical threshold used to prevent division by zero. Vortices are theoretically identified as spatially connected regions satisfying $\varOmega _r > 0.5$ when $\varepsilon _{vort} = 0$. In practice, however, vortices are identified using a small $\varepsilon _{vort} > 0$ and $\varOmega _r \geq 0.52$ (Liu et al. 2019a) to e.g. remove weak vortices.

1.3. Contributions

In this paper, we utilize the advantageous properties of geometry-based analyses of the VGT to efficiently characterize turbulence initiated by a vortex ring collision. We use the adaptive, multiresolution computational techniques discussed in § 2 to perform a direct numerical simulation of this flow at $Re_{\varGamma _0} = 4000$. In § 3, we establish the fidelity of our simulation and we visualize and discuss the various regimes of its evolution. In § 4, we analyse the partitioning of the velocity gradients to characterize these regimes in terms of the modes of deformation. In § 5, we introduce a geometry-based phase space that characterizes the action of the elliptic instability and its interplay with other mechanisms driving the turbulent flow. Our analyses reveal statistical features of the VGT that are similar to those of previous simulations. They also provide tools with the potential to help disentangle mechanisms underlying vortex interactions during transition and turbulent decay. Finally, we summarize our results in the context of previous works and highlight promising future research prospects in § 6.

2. Methods

2.1. Computational method

To efficiently simulate a turbulent vortex ring collision, we adopt a recently developed multiresolution solver for viscous, incompressible flows on unbounded domains (Liska & Colonius 2016; Dorschner et al. 2020; Yu 2021; Yu, Dorschner & Colonius 2022). Yu et al. (2022) provide a detailed discussion of the formulation, properties and performance of the method. We summarize the key advantages of the solver here and expound the computational formulation in Appendix A. The advantages we discuss allow us to simulate a relatively high Reynolds number vortex ring collision at a relatively low computational cost.

The Navier–Stokes equations (NSE) are spatially discretized onto a staggered Cartesian grid using a second-order-accurate finite-volume scheme that endows discrete operators with useful properties (i) (Liska & Colonius 2016). Discrete differential operators are constructed to mimic the symmetry, orthogonality and integration properties of their continuous counterparts. They also commute with the Laplacian and integrating factor operators, as defined in Appendix A. Furthermore, the discretization of the nonlinear term in the momentum equations preserves relevant (e.g. energy) conservation properties in the absence of viscosity. Together, the mimesis, commutativity and conservation properties of the discretization scheme facilitate fast, stable, high-fidelity simulations of turbulent flows.

The computational methods we employ also have high parallel efficiency (ii) and linear algorithmic complexity (iii). The computational efficiency of the flow solver is primarily centred around solving the discrete pressure Poisson equation on a formally unbounded grid using the lattice Green's function (LGF) (Liska & Colonius 2014, 2016; Liska 2016). The Poisson equation is obtained by taking the divergence of the NSE in rotational form, such that the source term is $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {r}$, where $\boldsymbol {r} = \boldsymbol {u} \times \boldsymbol {\omega }$ is the Lamb vector. By considering flows with at least exponentially decaying far-field vorticity, the approximate support of this source field can be captured using a finite computational domain. Given a source cutoff threshold, the finite domain is adaptively truncated to capture only the regions relevant to the Poisson problem. Solving the Poisson problem over this domain involves the convolution of the LGF with the source field. The flow solver achieves (ii) and (iii) by efficiently evaluating this convolution via a fast multipole method (Liska & Colonius 2014) that compresses the kernel using polynomial interpolation. This method is accelerated by exploiting the efficiency of fast Fourier transforms on a block-structured Cartesian grid.

In addition to spatially adapting the extent of the computational domain, adaptive multiresolution discretization (iv) is achieved by using adaptive mesh refinement (AMR) to reduce the number of degrees of freedom required for solutions. As discussed previously (Dorschner et al. 2020; Yu et al. 2022), the present AMR framework is carefully constructed to preserve the desirable operator properties (i) and augment the efficiencies (ii, iii) associated with the uniform-grid framework (Liska & Colonius 2014, 2016). In the AMR framework, the computational grid is partitioned into multiple levels, each with double the resolution in each direction as the previous level. The spatial regions associated with each level are non-overlapping, except for extended regions that are used to compute a combined source term that includes a correction induced by the difference between the coarse-grid and fine-grid partial solutions. As formulated in Appendix A, a region is refined when its combined source exceeds a threshold and it is coarsened when its combined source falls below a smaller threshold. As shown in § 3.1 (see table 1), this AMR formulation drastically reduces the number of computational cells required to capture a head-on vortex ring collision compared with a fixed-resolution scheme.

Table 1. Reference times used to analyse the initial (red), transitional (green) and turbulent (blue) regimes of the present vortex ring collision, where $t^* = 14.77$ is the time of maximum dissipation. Here, $\bar {R}_p$ represents the mean vortex ring radius (see Appendix B), $V$ is the volume of the computational domain and $\phi _k$ is the fraction of $V$ occupied by level $k$ of the AMR grid, which has grid spacing $\Delta x_k = \Delta x_{base}/2^{k}$ with $\Delta x_{base} = 0.04$.

2.2. Vortex ring collision simulation

As depicted in figure 1, we consider a flow configuration in which the vortex rings are initialized with opposing circulations such that they propagate toward one another along the $z$-axis and meet at the collision plane at $z = 0$. The rings are initialized a distance $L_z = 2.5 R_0$ apart, which is sufficiently large to mitigate their mutual influence during the most vigorous period of equilibration. Both rings are initialized with Gaussian vorticity distributions (1.1) such that $Re_{\varGamma _0} = 4000$ and $\delta _0 = 0.2$. Unless otherwise stated, we use the initial circulation, $\varGamma _0 = 1$, and radius, $R_0 = 1$, of each ring to non-dimensionalize all variables. To excite transition, we randomly perturb the radii of the vortex rings using the first 32 Fourier modes in $\theta$, which are prescribed random phases and uniform magnitudes, $R_{pert} = 5 \times 10^{-4}$. Consistent with previous tests (Yu et al. 2022), these initial perturbations are sufficiently large to dominate perturbations incurred by discretization errors.

Figure 1. Initial geometry of the flow configuration used to simulate the head-on collision between vortex rings. The shading of the vortex cores reflects their Gaussian vorticity profiles.

The computational mesh we use has $N_{level} = 2$ levels of refinement beyond the base level such that the ratio of the coarsest-grid spacing to the finest-grid spacing is $\Delta x_{base} / \Delta x_{fine} = 4$. Based on preliminary simulations of turbulent vortex rings (Liska & Colonius 2016) and vortex ring collisions (Yu et al. 2022), we select $a_0/\Delta x_{base} = 5$ and $\Delta t/ \Delta x_{fine} = 0.35$ to ensure the flow is well resolved throughout the simulation. Finally, parameters controlling the spatial and mesh refinement thresholds are chosen as discussed in Appendix A.

2.3. Simulation integral metrics

We track the evolution and fidelity of the simulation using integral metrics associated with incompressible flows (Liska & Colonius 2016). Particularly, we compute the hydrodynamic impulse, helicity, vortical kinetic energy and enstrophy of the flow, which are denoted by $\boldsymbol {I}_{V}$, $H$, $K_V$ and $E$, respectively. These integrals are formally evaluated on an unbounded domain, but we evaluate them using the finite AMR grid as

(2.1)\begin{equation} \left.\begin{gathered} \boldsymbol{I}_{V}(t) =\int_{V(t)} (\boldsymbol{x} \times \boldsymbol{\omega})\,{\rm d} V,\quad H(t) = \int_{V(t)} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\omega} )\,{\rm d} V, \\ K_{V}(t) = \int_{V(t)} \boldsymbol{u} \boldsymbol{\cdot} (\boldsymbol{x} \times \boldsymbol{\omega})\,{\rm d} V,\quad E(t) = \tfrac{1}{2} \int_{V(t)} \lvert \boldsymbol{\omega} \rvert^2 \,{\rm d} V, \end{gathered}\right\} \end{equation}

where $V(t)$ is the time-varying AMR grid. The impulse is the appropriate measure of momentum since it converges for flows on unbounded domains with compact vorticity. The vortical kinetic energy can also be expressed as $K_V = K + K_{\partial V}$, where $K$ represents the kinetic energy and $K_{\partial V}$ is a correction term based on the flow at the boundary of the grid, $\partial V(t)$. These metrics can be expressed as

(2.2a,b)\begin{equation} K(t) = \tfrac{1}{2} \int_{V(t)} \lvert \boldsymbol{u} \rvert^2\,{\rm d} V,\quad K_{\partial V}(t) = \int_{\partial V(t)} \boldsymbol{x} \boldsymbol{\cdot} \big((\boldsymbol{u} \boldsymbol{u} ) \boldsymbol{\cdot} \boldsymbol{n} - \tfrac{1}{2} \lvert \boldsymbol{u} \rvert^2 \boldsymbol{n}\big)\,{\rm d} S, \end{equation}

where $\boldsymbol {n}$ is the normal vector of $\partial V$ (Wu, Ma & Zhou 2015). For vanishing far-field velocity, $K_{\partial V}$ vanishes on unbounded domains and, for the present grid, we make use of the smallness of $K_{\partial V}$ when analysing dissipation.

In the absence of non-conservative external body forces, the hydrodynamic impulse is conserved for incompressible flows on unbounded domains (Saffman 1993). The helicity would also be conserved in the absence of viscosity, and it is useful for assessing simulation fidelity as the vortex rings initially approach the collision plane since the evolution of the flow is dominated by inviscid effects. These integral metrics initially evaluate to $\boldsymbol{I}_{V}(0) = \boldsymbol{0}$ and $H(0) = 0$ due to the spatial symmetries of the initial flow configuration. These initial symmetries hold to the extent that the vorticity is well captured and the contributions of the random perturbations used to excite instability growth are negligible. As the flow evolves, subsequent deviations from these initial values reflect the degree to which the corresponding initial symmetries are broken.

The enstrophy and the kinetic energy provide a more detailed picture of simulation fidelity during transition and turbulent decay, when viscous dissipation at small scales becomes relevant. For unsteady, incompressible flows on unbounded domains, the dissipation governs the decay rate of kinetic energy and can be expressed in terms of the enstrophy. Comparing these integrals is useful for characterizing (i) the degree to which small-scale features are resolved during peak dissipation and (ii) the flux of kinetic energy out of the finite computational domain (Archer et al. 2008). We therefore introduce effective Reynolds numbers, which are given by

(2.3a,b)\begin{equation} \frac{Re^{eff}_S (t)}{Re_{\varGamma_0}} =-\frac{\varPhi_S (t)}{{\rm d} K/{\rm d} t},\quad \frac{Re^{eff}_W (t)}{Re_{\varGamma_0}} =-\frac{\varPhi_W (t)}{{\rm d} K/{\rm d} t}, \end{equation}

where $\varPhi _S$ is the volume-integrated dissipation and $\varPhi _W = 2 E/Re_{\varGamma _0}$ is its enstrophy-based counterpart (Serrin 1959). Here, we differentiate $K$ instead of $K_V$ when computing the effective Reynolds numbers to prevent amplification of the noise associated with adaptations in the computational domain, to which $K_V$ is more sensitive. This is justified since $K$ and $K_V$ are nearly identical throughout the present simulations (see figure 2). The ratio $Re^{eff}_S$ is useful for assessing spatial resolution since the dissipation can vary significantly during transition and turbulent decay. The corresponding Kolmogorov scale, $\eta = (\nu ^3/\varPhi _S)^{1/4}$, can also be used to validate the selected grid spacings. The difference between $Re^{eff}_S$ and $Re^{eff}_W$ reflects the relative significance of the acceleration of the flow on $\partial V$ through the boundary integral in the Bobyleff–Forsyth formula (Serrin 1959). Together, the error metrics defined in this section comprehensively characterize the fidelity of the simulation as its flow structures evolve and the computational domain adapts accordingly.

Figure 2. Temporal evolution of the integral metrics defined in § 2.3 over the course of the simulation. The vertical lines correspond to the reference times in table 1 and they are coloured accordingly. The horizontal lines in the $Re^{eff}$ panel represent $Re_{\varGamma _0}$ (solid) with 10 % margins (dashed). The impulse magnitude is normalized by that of each vortex ring in isolation, $\lvert \boldsymbol {I}_{V1} \rvert \approx 1.02{\rm \pi} \approx 3.204$. The enstrophies $E_{\mathcal {E}}$ and $E_{\mathcal {C}}$ are computed using vorticities located at the edges and centres, respectively, of the computational cells.

3. Evolution of integral metrics and vortical structures

3.1. Evolution of integral metrics

Figure 2 shows the evolution of the integral metrics from § 2.3 over the course of the simulation. In the subsequent analysis, we reference the various regimes of flow development with respect to the time, $t^* = 14.77$, at which maximum dissipation is attained. Table 1 qualitatively characterizes the state of the simulation at each reference time we consider for the initial, transitional and turbulent regimes of the simulation.

The initial evolution of the flow involves a rapid period of equilibration ($t \lesssim 0.25 t^*$) and the propagation of the equilibrated rings towards the collision plane ($0.25 t^* \lesssim t \lesssim 0.50 t^*$). The interaction of the rings accelerates their radial expansion ($0.50 t^* \lesssim t \lesssim 0.75 t^*$) and the elliptic instability eventually emerges along the expanding rings ($0.75 t^* \lesssim t \lesssim 0.90 t^*$). Appendix B supports the importance of the elliptic instability during the early stages of transition. Subsequently, the flow transitions to turbulence ($0.90 t^* \lesssim t \lesssim t^*$) and rapidly produces small-scale flow structures. Following transition, the flow undergoes turbulent decay for the remainder of the simulation (i.e. for $t \gtrsim t^*$). See § 3.2 for visualizations of the flow at the reference times from table 1 associated with each of these regimes of evolution.

As the vortex rings initially propagate towards the collision plane ($t \lesssim 0.50 t^*$), the kinetic energy decays slowly and the enstrophy and dissipation are relatively small. The effective Reynolds numbers rapidly adjust to the value of $Re_{\varGamma _0}$ during the initial equilibration period ($t \lesssim 0.25 t^*$) and remain roughly constant as the equilibrated rings approach the collision plane ($0.25 t^* \lesssim t \lesssim 0.50 t^*$). The helicity is well conserved in this regime since the flow evolves in a nearly inviscid fashion. Further, the impulse is initially small and grows relatively slowly during this period. These results suggest that the symmetries associated with the handedness and momentum distribution of the flow are well preserved in the initial regime of evolution.

As the rings expand radially at the collision plane ($0.50 t^* \lesssim t \lesssim 0.75 t^*$) and the elliptic instability emerges ($0.75 t^* \lesssim t \lesssim 0.90 t^*$), the kinetic energy decays more rapidly and the dissipation grows. During these periods, the helicity symmetry remains well preserved and the effective Reynolds numbers remain relatively constant near $Re_{\varGamma _0}$, suggesting that the flow is well resolved. However, the impulse magnitude varies more rapidly in time due to the rapid radial expansion of the rings. In following the expanding vortical flow at the collision plane, the adaptations of the domain break the symmetry associated with impulse integral more significantly than during the initial evolution of the rings. The resulting growth in $\lvert \boldsymbol {I}_V \rvert$ is primarily attributed to its component in the $z$ direction, along which the domain is compressed as the flow concentrates about the collision plane.

As the flow transitions to turbulence ($0.90 t^* \lesssim t \lesssim t^*$), the kinetic energy decays even more rapidly and the dissipation approaches its maximum value. Due to the proliferation of small-scale flow structures during this period, the effective Reynolds numbers drop to their minimum values at $t \approx t^*$, when the flow is most difficult to resolve. The increased difference between $Re^{eff}_S$ and $Re^{eff}_W$ reflects that the acceleration of the flow near $\partial V$ is more relevant at this time. The rapid generation of small-scale flow structures also implies that viscosity plays a more important role in this regime. Correspondingly, the helicity begins to vary in time in this regime, reaching its maximum rate of change at the time of peak dissipation. Its variations reflect that vortex lines in the flow undergo rapid topological changes during transition. By contrast, the impulse magnitude decays to a roughly constant value as the radial expansion of the rings slows in this transitional regime.

During the turbulent decay of the flow ($t \gtrsim t^*$), the kinetic energy becomes small and the dissipation decays rapidly, eventually falling below its initial value (at $t \approx 1.62 t^*$). The dissipation matches the kinetic energy decay rate more closely for this regime than for transition. Further, as the turbulence develops, the effective Reynolds numbers agree well with one another and, to a lesser extent, with $Re_{\varGamma _0}$. These features reflect, respectively, that the acceleration near $\partial V$ is less significant and that the small scales are relatively well resolved, especially with respect to the transitional period. The helicity variations in this turbulent regime also eventually slow relative to those observed during transition. Similarly, the impulse remains roughly constant around its value at $t = t^*$. Whereas the $z$ component dominates the impulse magnitude during the radial expansion of the rings, all impulse components have similar magnitudes in this turbulent regime.

The evolution of the integral metrics characterizes the various regimes of the flow and supports the fidelity of our simulation. For example, the helicity is well conserved during the nearly inviscid evolution of the vortex rings and its subsequent variations are relatively small in magnitude. Further, the variations in the impulse magnitude throughout the simulation remain less than 5 % of the impulse associated with each vortex ring in isolation. These results indicate that the symmetries associated with the handedness and momentum distribution of the flow also remain well preserved in the appropriate regimes of the simulation.

Our dissipation analysis also suggests that the small-scale flow structures remain reasonably well resolved throughout the entire simulation. After equilibration, the maximum relative deviation in $Re^{eff}_S$ from $Re_{\varGamma _0}$ is roughly 6.5 % and it occurs around $t \approx t^*$. This relative deviation is similar to that of a previous simulation of a single vortex ring at $Re_{\varGamma _0} = 7500$ using a finite computational grid (Archer et al. 2008). Moreover, it is considerably smaller during the approach and radial expansion of the rings and, to a lesser extent, during turbulent decay. Even during peak dissipation, when the Kolmogorov scale is smallest, the finest grid has acceptable resolution since $\Delta x_{fine} / \eta \approx 3.42$. Altogether, these results suggest that our simulation is well resolved and support our analysis of the mechanisms underlying transition and turbulent decay.

3.2. Evolution of vortical flow structures

For the present simulation, we identify vortices using the $\varOmega _r$ criterion (Liu et al. 2019a) with a numerical threshold of $\varepsilon _{vort}$ = 0.04. This criterion provides connections to the triple decomposition of the VGT and the structure of local streamlines. Due to the well-preserved symmetries of the flow, the global spatial mean of the vorticity tensor is nearly zero and, hence, the $\varOmega _r$ criterion is nearly objective (Liu et al. 2019c) for the present simulation. We specifically visualize the flow using $\varOmega _r = 0.52$ and $\varOmega _r = 0.93$ to investigate the structures of the vortex boundaries and the vortex cores, respectively. We colour these structures using $\cos \theta ^*$, where $\theta ^*$ is the angle between the $z^*$-axis and the $z$-axis. This colour scheme enables the identification of antiparallel vortices along the $z$-axis, which play an important role in mediating transition and generating small-scale flow structures in the present vortex ring collision. In figure 3, we visualize the vortical structures in the flow at each reference time from table 1.

Figure 3. Visualizations of the vortex boundaries ($\varOmega _r = 0.52$, left side) and vortex cores ($\varOmega _r = 0.93$, right side), coloured by $\cos \theta ^*$, for each reference time from table 1. A movie depicting the evolution of the vortex boundaries from the auxiliary viewpoint (leftmost column) is provided as supplementary material available at https://doi.org/10.1017/jfm.2024.90.

During the initial evolution of the flow, the equilibrated vortex rings approach the collision plane and expand radially due to their mutual interaction. In this regime, the thinning of the vortex boundaries and cores illustrates the mechanisms driving the shift from a rigid-rotation-dominated regime to a shearing-dominated regime. Further, the visualizations at $t = 0.75 t^*$ depict the emergence of the short-wave elliptic instability, which is consistent with previous vortex ring collision simulations in similar parameter regimes (McKeown et al. 2020; Mishra et al. 2021).

The transitional regime of the flow is marked by the development of secondary vortex filaments and the subsequent generation of small-scale vortical flow structures. At $t = 0.90 t^*$, the visualizations show the development of secondary vortical structures around the circumference of the collision. These structures consist of antiparallel vortex filament pairs that arise in regions where the elliptic instability drives local interactions between the rings. This behaviour supports the notion that the elliptic instability mediates the initial transition of the rings, leading to the development of secondary vortical structures. The antiparallel secondary filaments become increasingly densely packed as transition progresses and they mediate the proliferation of small-scale vortical flow structures, e.g. as observed at $t = t^*$.

In Appendix B, we decompose the flow into azimuthal Fourier modes to characterize the wavenumbers of the perturbations that dominate transition. Our analysis confirms that the short-wave elliptic instability, with wavelength of the order of the core radius, mediates the initial stages of transition. We further show that, at $t = t^*$, the most prominent antiparallel vortices occur at the second harmonic of an originally dominant perturbation. Taken together, our visualizations and perturbation analysis are consistent with the initial stages of the iterative elliptic instability pathway, which is driven by subsequent generations of antiparallel vortex filaments (McKeown et al. 2020).

During the turbulent decay of the flow, the geometric features of the vortex boundaries remain similar at each reference time. However, as energy is dissipated, the smallest-scale vortices are progressively destroyed and the vortical flow structures grow larger in time. The structures of the vortex cores and boundaries reinforce the importance of the interactions between the secondary vortex filaments in mediating the evolution of the turbulent flow. The vortex boundaries also show the formation and ejection of vortex rings from the turbulent cloud resulting from the collision. These ejections, which are a hallmark of the Crow instability, often occur in regions where antiparallel vortex filaments interact and are of similar size to those filaments. This observation provides further evidence of the interplay between the elliptic and Crow instabilities driven by interacting vortex filaments around the turbulent cloud (Mishra et al. 2021; Ostilla-Mónico et al. 2021). In what follows, we develop machinery to probe these mechanisms in the context of features of the velocity gradients, with a particular emphasis on characterizing the action of the elliptic instability among other mechanisms.

4. Partitioning of velocity gradients

Here, we investigate the partitioning of the velocity gradients to characterize the evolution of the flow in the context of the fundamental modes of deformation. We first consider volumetric weighted averages of the relative contributions of various constituents of $\tilde {\!\boldsymbol{\mathsf{A}}}$ to the strength of the velocity gradients. These averages may be expressed as

(4.1)\begin{equation} \langle {\xi}\rangle_{\mkern-3mu A^2} = \frac{\displaystyle\int_{V(t)} A^2 \xi \,{\rm d} V}{\displaystyle\int_{V(t)} A^2\,{\rm d} V}, \end{equation}

where $\xi \in \{\tilde {S}^2, \tilde {W}^2\}$ for the Cauchy–Stokes decomposition and $\xi \in \{\dot {\epsilon }^2, \dot {\gamma }^2, \dot {\varphi }^2, 2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} \}$ for the triple decomposition. Here, we have used that the Frobenius inner product, denoted by :, of a symmetric tensor with an antisymmetric tensor is zero and that $\dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W} = \dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu S} = \dot {\gamma }^2/2$. The shear–rotation correlation term, $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} = 2 \dot {\varphi }_{z^*} \dot {\gamma }_{z^*} > 0$, reflects the presence of shearing in the plane of rigid rotation. All of the relative contributions discussed are unitarily invariant and, thus, they apply to both the principal coordinates and the global coordinates. Figure 4 shows how these relative contributions evolve during the simulation for both decompositions of $\tilde {\!\boldsymbol{\mathsf{A}}}$.

Figure 4. Relative contributions of the constituents of the Cauchy–Stokes decomposition and the triple decomposition to $A^2$. The vertical lines correspond to the reference times in table 1 and they are coloured accordingly.

Consistent with the equivalence of $\varPhi _S$ and $\varPhi _W$ for incompressible flows on unbounded domains (Serrin 1959), $\langle {\tilde {S}^2}\rangle _{\mkern -3mu A^2} \approx \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2} \approx 0.50$ for the present simulations. The largest deviations from this balance occur during equilibration ($t \lesssim 0.25 t^*$) and around the time of peak dissipation ($t \approx t^*$). These deviations are consistent with the behaviour of the effective Reynolds numbers in figure 2 and their smallness further validates the ability of the finite computational grid to approximate a formally unbounded flow. However, since $\langle {\tilde {S}^2}\rangle _{\mkern -3mu A^2}$ and $\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$ remain relatively constant throughout the simulation, they provide limited information about the nature of the velocity gradients as the flow progresses through its initial, transitional and turbulent regimes.

Compared with the constituents of the Cauchy–Stokes decomposition, the constituents of the triple decomposition show more pronounced variations associated with the different regimes of evolution. For the initial Gaussian vorticity profiles, the contribution of rigid rotation to the enstrophy dominates the contribution of antisymmetric shearing. During equilibration ($t \lesssim 0.25 t^*$), the fluctuations in all contributions of the triple decomposition constituents reflect the redistribution of velocity gradients in the cores of the vortex rings. As the equilibrated rings approach the collision plane and spread ($0.25 t^* \lesssim t \lesssim 0.75 t^*$), $\langle {\dot {\epsilon }^2}\rangle _{\mkern -3mu A^2}$ and $\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2}$ decrease and $\langle {\dot {\gamma }^2}\rangle _{\mkern -3mu A^2}$ and $2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ increase. As the elliptic instability emerges, these contributions level off in a regime where antisymmetric shearing dominates rigid rotation and shear–rotation correlations are enhanced. The subsequent development of the elliptic instability ($0.75 t^* \lesssim t \lesssim 0.90t^*$) is marked by slight rebounds in the contributions of $\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2}$ and $\langle {\dot {\gamma }^2}\rangle _{\mkern -3mu A^2}$. These rebounds are associated with the emergence of secondary vortex filaments and, hence, the nonlinear evolution of the elliptic instability. The transition to turbulence ($0.90 t^* \lesssim t \lesssim t^*$), which is associated with the generation of small scales and enhanced dissipation, is marked by a decrease in the contribution of $2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$.

Remarkably, even though the flow is not stationary during turbulent decay, the relative contributions of the constituents of the triple decomposition to the strength of the velocity gradients remain roughly constant after transition. Further, as summarized in table 2, these ‘equilibrium’ relative contributions are similar to those computed by Das & Girimaji (2020b) for forced isotropic turbulence at high Taylor-scale Reynolds numbers. This agreement suggests that the velocity gradient partitioning may encode a relatively common balance in unbounded, incompressible turbulence with appropriate symmetries. In this balance, shearing makes the largest contribution to the velocity gradients and rigid rotation makes the smallest contribution.

Table 2. Comparison of the equilibrium partitioning of the velocity gradients for the present vortex ring collision with the partitioning computed for forced isotropic turbulence (Das & Girimaji 2020b). Here, the equilibrium partitioning is computed as the mean over the turbulent decay regime ($t \gtrsim t^*$) and it is insensitive to the length of the averaging interval.

Beyond the strength of velocity gradients, it is also useful to examine the interplay between the modes of deformation in the context of vortical flow structures. Here, we introduce a new phase space defined by the relative contributions of $\langle {\dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W}}\rangle _{\mkern -3mu A^2}$, $\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2}$ and (implicitly) $2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ to $\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. The upper bound of this phase space is found by maximizing $\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$, which occurs when $2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2} = 0$. The bottom boundary is found by minimizing $\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$, which occurs when $2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2} / \langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2} = 2 \sqrt{\langle {\dot {\varphi }^2}\rangle _{\mkern -3mu A^2} \langle {\dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W}}\rangle _{\mkern -3mu A^2}}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. Correspondingly, these boundaries may be expressed as

(4.2) \begin{align} &\left(1-\sqrt{\frac{\langle {\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}\rangle_{\mkern-3mu A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}}\right)^2 \leq \frac{\langle {\dot{\varphi}^2}\rangle_{\mkern-3mu A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}\ \, \, \leq 1 - \frac{\langle {\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}\rangle_{\mkern-3mu A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}, \end{align}

(4.3) \begin{align} &\Updownarrow \qquad\qquad \Updownarrow \nonumber\\ &-2 \left(\frac{\langle {\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}\rangle_{\mkern-3mu A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}} - \sqrt{\frac{\langle {\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}\rangle_{\mkern-3mu A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}}\right) \geq \frac{2 \langle {\boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}}}\rangle_{\mkern-3mu A^2}}{\langle {\tilde{W}^2}\rangle_{\mkern-3mu A^2}}\!\geq 0, \end{align}

highlighting that lower and upper bounds of the rigid rotation contribution correspond to the upper and lower bounds, respectively, of the shear–rotation correlation contribution. The maximum value of $2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$ varies along the lower boundary of this phase space to ensure that the relative contributions sum to unity. The global maximum occurs when $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2 = 0.5$ or, equivalently, $\dot {\varphi }^2 / \tilde {W}^2 = \dot {\gamma }^2_{\rule {0pt}{-2ex} \mkern -2mu W} / \tilde {W}^2 = 0.25$ at all points, which corresponds to the (pointwise) maximum of $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} = ( \sqrt {2} + 1 )^{-1}$ reported by Das & Girimaji (2020b). This maximum corresponds to local streamlines in the principal frame for which shearing occurs exclusively in the plane of rigid rotation.

Figure 5 depicts the trajectory of the flow in this shear–rotation phase space and elucidates how the relationships between the constituents of enstrophy associated with the triple decomposition evolve in time. Following rapid variations during equilibration, the trajectory returns to a position in phase space similar to that of the initial condition at $t \approx 0.25 t^*$. The trajectory then undergoes a shift across the phase space as the equilibrated rings approach the collision plane and spread radially ($0.25 t^* \lesssim t \lesssim 0.75 t^*$). This shift from a rigid-rotation-dominated regime to a shearing-dominated regime is associated with an enhanced contribution from $2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. The development of the elliptic instability ($0.75 t^* \lesssim t \lesssim 0.90 t^*$) and transition ($0.90 t^* \lesssim t \lesssim t^*$) are associated with a shift in the direction of the trajectory towards smaller contributions of $2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}/\langle {\tilde {W}^2}\rangle _{\mkern -3mu A^2}$. During turbulent decay ($t \gtrsim t^*$), the trajectory remains roughly fixed in the phase space and very close to its position at $t = t^*$.

Figure 5. Shear–rotation phase space trajectory of the flow, highlighting the evolution during equilibration (left inset) and transition and turbulent decay (right inset). The contours represent increments of 0.05 in the colour scale. The white circle marks the initial condition and the fill colours of the black circles correspond to the reference times they represent from table 1.

Considering this phase space trajectory in the context of the dissipation (see figure 2) reveals that, while the initial growth in dissipation is associated with enhanced shear–rotation correlations, its subsequent enhancement during transition is associated with a reduction in shear–rotation correlations. In Appendix C, we reexamine the visualizations from figure 3 in terms of the shear–rotation correlations in the flow to highlight their relationship to the vortical flow structures. In § 5, we interpret the effects of these shear–rotation correlations using a new, related phase space, based on local streamline geometry, to characterize the elliptic instability and other mechanisms.

5. Statistical geometry of local streamlines

5.1. Phase space transformations

The elliptic instability is associated with the resonance of the vortical flow with the underlying strain field and acts to break up elliptic streamlines. Consistent with this picture, we introduce a new geometry-based phase space that captures local flow features that (i) are conducive to the elliptic instability and (ii) characterize its action.

To address (i), we consider the angle, $\theta _\omega$, between the vorticity vectors associated with antisymmetric shearing and rigid rotation, which is given by

(5.1)\begin{equation} \theta_\omega = \cos^{-1}\left(\sqrt{\frac{(\boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}})^2}{\dot{\varphi}^2 \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}}\right)\!. \end{equation}

Our focus on shear straining is consistent with the classical models of strained vortices used to characterize the elliptic instability (Kerswell 2002). Since decreasing $\theta _\omega$ corresponds to increasing the alignment between shearing and rigid rotation, it can be associated with conditions conducive to the elliptic instability.

To address (ii), we consider the aspect ratio, $\zeta$, of the elliptic component of rotational local streamlines in the plane of rigid rotation, which is given by

(5.2)\begin{equation} \zeta = \sqrt{1-e^2} = \sqrt{\frac{\dot{\varphi}^2}{\dot{\varphi}^2 + 2 \boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}}}}, \end{equation}

where $e$ represents the eccentricity of an ellipse with aspect ratio $\zeta$. This aspect ratio characterizes the scale-invariant geometry of the local streamlines in the plane of rigid rotation. As such, it can be used alongside $\theta _\omega$ to characterize the action of the elliptic instability by identifying how alignment between shearing and rigid rotation affects local streamline geometry.

The present shear–rotation and geometry-based phase spaces can be understood through nonlinear transformations of the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. The transformations we derive express the relative contributions of the triple decomposition to the velocity gradients using $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$ and $r_{\rule {0pt}{1.6ex} \mkern -4mu A}$, and they represent the inverse transformations to those presented by Das & Girimaji (2020b). However, an additional parameter, $\theta _\omega$, is generally required to evaluate our transformations. This extra parameter demonstrates that the invariants of $\tilde {\!\boldsymbol{\mathsf{A}}}$ alone cannot generally characterize the relative contributions of the constituents of the triple decomposition to the velocity gradients.

For rotational local streamlines, the transformations are given by

(5.3)\begin{align} \left.\begin{gathered} \dot{\epsilon}^2 = \frac{(-3 \times 2^{2/3} q_{\rule{0pt}{1.6ex} \mkern-3mu A} + 3^{2/3} (2 \sqrt{3 \varDelta} + 9 \lvert r_{\rule{0pt}{1.6ex} \mkern-4mu A} \rvert)^{2/3})^2}{6^{5/3} (2 \sqrt{3 \varDelta} + 9 \lvert r_{\rule{0pt}{1.6ex} \mkern-4mu A} \rvert)^{2/3}},\quad \dot{\gamma}^2 = 2 \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W} = 1 - 2 \dot{\epsilon}^2 - 2 q_{\rule{0pt}{1.6ex} \mkern-3mu A}, \\ \dot{\varphi}^2 = \left(-\sqrt{\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}} \cos \theta_{\omega} + \sqrt{1 - \dot{\epsilon}^2 - \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W} (2 - \cos^2 \theta_{\omega})}\right)^2,\quad 2 \boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}} = 2 \sqrt{\dot{\varphi}^2 \dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}} \cos \theta_{\omega}, \end{gathered}\right\} \end{align}

where $\varDelta = q_{\rule {0pt}{1.6ex} \mkern -3mu A}{\mkern -6mu}^3 + \frac {27}{4} r_{\rule {0pt}{1.6ex} \mkern -4mu A}{\mkern -6mu}^2$ is proportional to the discriminant of the characteristic equation of  $\tilde {\!\boldsymbol{\mathsf{A}}}$. In this case, $\dot {\epsilon }^2$ and $\dot {\gamma }^2$ can be determined directly from $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$ and $r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ and $\dot {\varphi }^2$, $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and $\zeta$ can be determined if $\theta _\omega$ is known. For non-rotational local streamlines, the transformations are given by

(5.4a–c)\begin{equation} \dot{\epsilon}^2 =-2 q_{\rule{0pt}{1.6ex} \mkern-3mu A},\quad \dot{\gamma}^2 = 1- \dot{\epsilon}^2,\quad \dot{\varphi}^2 = 0, \end{equation}

which can be determined directly from a single parameter, $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$. The rotational and non-rotational transformations are continuous with one another at their boundary ($\varDelta = 0$) and they are both symmetric about the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis. However, the aspect ratio is only well defined for $\varDelta > 0$, consistent with our focus on rotational local streamlines.

In figure 6, we illustrate how $\dot {\epsilon }^2$, $\dot {\gamma }^2$, $\dot {\varphi }^2$, $2\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and $\zeta$ vary within the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. In this phase space, the rotational geometries are externally bounded by $3\sqrt {3}\lvert r_{\rule {0pt}{1.6ex} \mkern -4mu A} \rvert = ( 1 + q_{\rule {0pt}{1.6ex} \mkern -3mu A} ) ( 1 - 2 q_{\rule {0pt}{1.6ex} \mkern -3mu A} )^{1/2}$ and the non-rotational geometries are externally bounded by $q_{\rule {0pt}{1.6ex} \mkern -3mu A} = -\frac {1}{2}$. For the rotational geometries, we display $\dot {\varphi }^2$, $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and $\zeta$ for $\theta _\omega = 43.57^\circ$, which corresponds to the mean value at $t = 1.25 t^*$ (see figure 8) and approximates the equilibrium value during turbulent decay. We document how each of these quantities varies with $\theta _\omega$ in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space in Appendix D.

Figure 6. The $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (a) and the variations in $\dot {\epsilon }^2$, $\dot {\gamma }^2$, $\dot {\varphi }^2$, $2\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$ and $\zeta$ in this space when $\theta _\omega = 43.57^\circ$ is held constant. Panel (a) shows the boundaries corresponding to the symmetry-based (dotted green) and the geometry-based (solid blue) vortex criteria (see § 1.2). It also labels the four classes of non-degenerate local streamline topologies: stable-focus–stretching (SFS), unstable-focus–compression (UFC), stable-node–saddle-saddle (SNSS) and unstable-node–saddle-saddle (UNSS).

Figure 6 shows that, generically, the contribution of pure shearing ($\dot {\gamma }^2$) tends to dominate the velocity gradients near the origin of the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane. The contribution of normal straining ($\dot {\epsilon }^2$) grows large when the velocity gradients are dominated by the strain rate tensor. By contrast, for rotational geometries, the contribution of rigid rotation ($\dot {\varphi }^2$) grows large near the external boundary in regions where the vorticity tensor dominates. The contribution of shear–rotation correlations ($2\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}$) grows largest in the intermediate region of the phase space and it decays near the discriminant line and the external boundary. The aspect ratio $\zeta$ is unity along the external boundary and it decays to zero at the discriminant line.

The elliptic instability, which is relevant in strained vortical flows, is expected to be most active for an intermediate range of aspect ratios. Further, vortex stretching and squeezing, which are known to play important roles in turbulent flows, can most readily be associated with the SFS and UFC streamline topologies, respectively (Lopez & Bulbeck 1993). In the corresponding regions in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ plane, the transformations in figure 6 suggest that the interplay between shearing and rotation is pertinent to these fundamental turbulent processes. Within this context, the transformations have the potential to characterize the role of the elliptic instability, among other mechanisms, in mediating such processes. In what follows, we investigate this premise by analysing the evolution of the velocity gradient distributions in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$, shear–rotation and $\zeta$–$\theta _\omega$ phase spaces.

5.2. Phase space distributions

The joint probability density functions (p.d.f.s) of the normalized velocity gradients in the phase spaces we investigate encode information about the local streamline geometries. Although these phase spaces are related to one another, the choice of phase space plays an important role in interpreting the statistical distributions of the flow. In the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, incompressible turbulent flows often follow a near-universal teardrop-like distribution about the origin (Das & Girimaji 2019, 2020a). The shear–rotation phase space highlights the distribution of rotational streamline geometries and it characterizes the interplay of rigid rotation and antisymmetric shearing. The $\zeta$–$\theta _\omega$ phase space also considers rotational streamline geometries and it characterizes the flow in terms of geometric features of local streamlines that are associated with the elliptic instability.

Figure 7 shows the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ and the shear–rotation phase space distributions at reference times, selected from table 1, that pertain to the development of the elliptic instability, transition and turbulent decay. As the elliptic instability emerges ($t \approx 0.75 t^*$), the velocity gradients are concentrated near the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis. Since $2\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ is relatively large around this time, the rotational regions of the flow are concentrated near the bottom boundary of the shear–rotation phase space, particularly near the location where $2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ is maximized.

Figure 7. Joint p.d.f.s of the velocity gradients satisfying $A^2/\textrm {max}(A^2) \geq 0.1\,\%$ in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (a) and $W^2/\textrm {max}(W^2) \geq 0.1\,\%$ in the shear–rotation phase space (b) at times $t = 0.75 t^*$, $0.90 t^*$, $1.00 t^*$ and $1.25 t^*$ (from top to bottom). The blue triangles represent the centroids of the distributions and the magenta contours represent the p.d.f. levels for which 90 % of the flow (by volume) resides at higher p.d.f. levels. These contours are smoothed by using coarser p.d.f. bins to ensure that they roughly enclose the regions with higher p.d.f. levels.

As the flow transitions ($t \approx 0.90t^*$), the distribution remains centred about the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis but it expands towards larger $\lvert r_{\rule {0pt}{1.6ex} \mkern -4mu A} \rvert$ as it begins to fill the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space. Similarly, the p.d.f. in the shear–rotation phase space begins extending away from the bottom boundary of the phase space, although its bulk remains concentrated at the boundary. The broadening p.d.f.s represent the generation of more diverse local streamline topologies (see figure 6), which is consistent with the formation of more complex vortical structures (see figure 3).

At the time of peak dissipation ($t \approx t^*$), the distributions populate nearly all of the area in both the phase spaces. During the subsequent turbulent decay (as shown for $t \approx 1.25 t^*$), the flow approaches its equilibrium distributions in both phase spaces, which remain similar to the distributions for $t \approx t^*$. In the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, the equilibrium p.d.f. above the $r_{\rule {0pt}{1.6ex} \mkern -4mu A}$-axis is concentrated slightly left of the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$-axis. This slight preference of the SFS topology is consistent with the typical presence of positive vortex stretching in regions of turbulent flows with rotational geometries. Below the $r_{\rule {0pt}{1.6ex} \mkern -4mu A}$-axis, the p.d.f. is concentrated along the discriminant line for $r_{\rule {0pt}{1.6ex} \mkern -4mu A} > 0$. The equilibrium distribution of our vortex ring collision in this phase space is similar to the near-universal teardrop-like shapes reported previously for forced isotropic turbulence (Das & Girimaji 2019, 2020a). This similarity suggests that, in addition to the velocity gradient partitioning (see table 2), the teardrop-like distribution may be more broadly applicable to incompressible flows with appropriate symmetries.

The p.d.f.s in the shear–rotation phase space evolve similarly to those in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (e.g. by broadening) and specifically highlight rotational geometries. However, the difference between vortex stretching and squeezing, which is encoded in the sign of the real eigenvalue of the VGT, cannot be distinguished in this phase space since its constituents are non-negative. The high concentration of the p.d.f.s near the lower corner of this phase space, representing the origin of the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, shifts the centroids accordingly and highlights the importance of shearing in the generation and evolution of turbulent flows. The equilibrium distribution remains relatively concentrated about the lower boundary of this phase space, including regions where 2$\langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ is relatively large. This behaviour highlights the potential for the elliptic instability to be active during turbulent decay, but the breadth of the distribution suggests it may not be a completely dominant mechanism driving turbulent flow in rotational regions.

To further investigate the evolution of the flow in the context of the elliptic instability, we again consider the relationship between the $\zeta$–$\theta _\omega$ and shear–rotation phase spaces. Using $\theta _\omega$ and $\zeta$, the constituents of the shear–rotation phase space and the corresponding shear–rotation correlation term and vortex identification criterion are given by

(5.5)\begin{equation} \left.\begin{gathered} \frac{\dot{\varphi}^2}{\tilde{W}^2} = \frac{4 \zeta^4 \cos^2 \theta_\omega}{1 + (-2 + 4\cos^2 \theta_\omega)\zeta^2 + \zeta^4},\quad \frac{\dot{\gamma}^2_{\rule{0pt}{-2ex} \mkern-2mu W}}{\tilde{W}^2} = \frac{(1 - \zeta^2)^2}{1 + (-2 + 4\cos^2 \theta_\omega)\zeta^2 + \zeta^4},\\ \frac{2 \boldsymbol{\dot{\varphi}} \mkern-5mu : \mkern-5mu \boldsymbol{\dot{\gamma}}}{\tilde{W}^2} = \frac{4 \zeta^2 \cos^2 \theta_\omega (1 - \zeta^2)}{1 + (-2 + 4\cos^2 \theta_\omega) \zeta^2 + \zeta^4}, \quad \varOmega_r = \left(1 + \left( 1 - \sqrt{\frac{\dot{\varphi}^2}{\tilde{W}^2}} \right)^2\right)^{-1}. \end{gathered}\right\} \end{equation}

Figure 8 shows the p.d.f.s of the velocity gradients alongside the distributions of $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}/\tilde {W}^2$ and $\varOmega _r$ in the geometry-based phase space at the same times as those in figure 7.

Figure 8. Joint p.d.f.s (a–d) of the velocity gradients in the $\zeta$–$\theta _\omega$ phase space at the same times, at the same points and in the same style as those in figure 7. The bottom plots superimpose the 90 % contours and centroids for $t = 0.75 t^*$ (red) and $t = 1.25 t^*$ (black) on $2 \boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }} / \tilde {W}^2$ (e) and on $\varOmega _r$ (f), as given by (5.5).

As the elliptic instability emerges ($t \approx 0.75 t^*$), the p.d.f. is highly concentrated at small values of $\theta _\omega$ over a broad range of $\zeta$. This distribution is consistent with the enhancement of $2 \langle {\boldsymbol {\dot {\varphi }} \mkern -5mu : \mkern -5mu \boldsymbol {\dot {\gamma }}}\rangle _{\mkern -3mu A^2}$ around this time (see figure 4) and reflects conditions conducive to the elliptic instability. The centroid of the distribution is located around $\zeta \approx 0.5$, which suggests that rotational streamlines with this aspect ratio may be particularly susceptible to these conditions.

During transition ($t \approx 0.90 t^*$), the distribution broadens to a much larger range of $\theta _\omega$ and, at the high end of this range, $\zeta$ becomes increasingly correlated with $\theta _\omega$. This broadening reflects the diversification of rotational local streamline topologies to include those for which shearing and rigid rotation are not well aligned. However, despite this broadening, the p.d.f. is still concentrated at small $\theta _\omega$. This behaviour is consistent with the importance of the elliptic instability during transition (see Appendix B).

During turbulent decay ($t \gtrsim t^*$), the features of the p.d.f. are relatively invariant in time. In the equilibrium distribution (approximated at $t \approx 1.25 t^*$), the centroid is located around $\theta _\omega \approx 44^\circ$, which is considerably larger than its value ($\theta _\omega \approx 3^\circ$) at $t = 0.75 t^*$. Further, the 90 % contour of the p.d.f. spans nearly the entire range of $\theta _\omega$ and highlights a well-defined sharpening in the correlation of $\zeta$ with $\theta _\omega$ with increasing $\theta _\omega$. The upper limit $(\theta _\omega,\zeta ) = (90^\circ,1)$ corresponds to circular local streamlines in the plane of rigid rotation subject to out-of-plane shearing.

The shifts in the centroid and the 90 % contour support the hypothesis that, while the elliptic instability still plays a role in the turbulent regime, other mechanisms also contribute significantly to the local flow structure. Specifically, they suggest that, in addition to the breakup of elliptic streamlines via the elliptic instability, the deformation of vortices, e.g. through the action of the Crow instability, may play an important role during turbulent decay. As discussed in Appendix C, the partitioning of shear–rotation correlations between the cores and boundaries of vortices in the flow (see figure 12) provides an interesting opportunity for analysing these mechanisms. Given current challenges in disentangling the elliptic and Crow instabilities in turbulent flows (Mishra et al. 2021; Ostilla-Mónico et al. 2021), the present geometry-based ($\zeta$–$\theta _\omega$) phase space has the potential to help distinguish flow features associated with these ubiquitous mechanisms.

Altogether, the results in this section reinforce the notion that the elliptic instability is the dominant mechanism mediating the transition of the present vortex ring collision. They also support the notion that, after transition, the elliptic instability is no longer a strictly dominant mechanism underlying the turbulent decay of the flow. The results point to increased contributions from rotational geometries with out-of-plane shearing, which may reflect interactions associated with mechanisms like the Crow instability. The results also highlight the ability of the new shear–rotation and geometry-based phases spaces to characterize the relative contributions of different modes of deformation and rotational features of local streamlines, respectively.

6. Conclusions

We use a recently developed adaptive, multiresolution numerical scheme based on the LGF to efficiently simulate the head-on collision between two vortex rings at a relatively high Reynolds number ($Re_{\varGamma _0} = 4000$). The fidelity of this simulation is confirmed using various integral metrics that reflect the symmetries, conservation properties and discretization errors of the flow. We provide a detailed analysis of the initial evolution, transition and turbulent decay of the flow to elucidate flow features that are pertinent to the mechanisms driving its evolution, e.g. the elliptic instability.

Our visualizations of vortex structures enable qualitative characterizations of the various regimes through which the flow evolves. They depict the short-wave elliptic instability as the mechanism driving the initial transition of the rings as they merge at the collision plane. Consistent with previous studies (McKeown et al. 2020; Mishra et al. 2021), late transition and (to a lesser extent) turbulent decay are mediated by antiparallel secondary vortex filaments that arise from local interactions associated with the elliptic instability. We confirm that the elliptic instability dominates transition by analysing the scales of dominant wave-like perturbations in that regime. During turbulent decay, we observe local ejections of vortex rings in regions where antiparallel vortex filaments interact. This observation supports the notion that interplay between the elliptic and Crow instabilities can impact vortex interactions, consistent with previous findings (Mishra et al. 2021; Ostilla-Mónico et al. 2021).

Our analysis of the flow centres around using the triple decomposition of the VGT to characterize the contributions of axial straining, shearing, rigid rotation and shear–rotation correlations to the velocity gradients. The mutual interaction of the rings is marked by the development of shearing-dominated vorticity and enhanced shear–rotation correlations, reflecting conditions conducive to the elliptic instability. These conditions are consistent with the initial elliptic instability observed in our visualizations and previously in similar configurations (McKeown et al. 2020; Mishra et al. 2021). During turbulent decay, the relative contributions of the different modes of deformation to the velocity gradient strength (which is not stationary) are roughly invariant in time, suggesting an equilibrium partitioning of the VGT. This equilibrium partitioning is remarkably similar to the partitioning observed for forced isotropic turbulence (Das & Girimaji 2020b), suggesting that it may provide a broadly applicable avenue for modelling incompressible flows with appropriate symmetries.

During the transition and turbulent decay of the flow, we also consider instantaneous distributions of the velocity gradients in various phase spaces. The broadening of the phase space distributions in these regimes reflects the generation of more diverse local streamline topologies. The distributions in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space show that the present vortex ring collision produces velocity gradients that follow the near-universal teardrop-like distribution observed previously for forced isotropic turbulence (Das & Girimaji 2019, 2020a). In addition to the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space, we introduce the shear–rotation phase space to characterize the interplay of shearing and rigid rotation in rotational settings and highlight the role of their correlations during transition and turbulent decay.

Finally, we introduce a geometry-based ($\zeta$–$\theta _\omega$) phase space to further characterize the action of the elliptic instability (and other mechanisms) during transition and turbulent decay. As the rings interact, the emergence of the elliptic instability spurring transition is associated with the alignment of shearing and rigid rotation ($\theta _\omega \approx 3^\circ$). In this regime, the elliptic local streamlines in the plane of rigid rotation have aspect ratios centred about $\zeta \approx 0.5$. During late transition and turbulent decay, the generation and interaction of secondary vortical structures broadens the distribution to include larger $\theta _\omega$, and the equilibrium distribution is ultimately centred near $\theta \approx 44^\circ$. In this regime, regions with high $\theta _\omega$ and high $\zeta$ become increasingly correlated as they approach ($\theta _\omega,\zeta$) = ($90^\circ,1$). In conjunction with our visualizations, these results suggest that proximity to vortex cores and boundaries may be a useful tool for modelling the interplay between mechanisms such as the elliptic and Crow instabilities. As a whole, the geometry-based phase space we introduce has the potential to help distinguish effects associated with the elliptic instability (small $\theta _\omega$) and other mechanisms, which is an ongoing challenge for turbulent flows driven by interacting vortex filaments (Mishra et al. 2021; Ostilla-Mónico et al. 2021).

Moving forward, the VGT phase spaces we introduce may provide a useful setting for analysing a broad class of turbulent flows. For vortex ring collisions, analysing regimes where the Crow instability dominates the elliptic instability would clarify the extent to which the phase spaces can disambiguate these mechanisms. More generally, it would also be useful to identify the conditions under which (i) the equilibrium partitioning of the VGT (Das & Girimaji 2020b) and (ii) the teardrop-like distribution in the $q_{\rule {0pt}{1.6ex} \mkern -3mu A}$–$r_{\rule {0pt}{1.6ex} \mkern -4mu A}$ phase space (Das & Girimaji 2019, 2020a) are applicable. The present VGT analyses are limited by the local, instantaneous nature of the streamline geometries under consideration. It would be interesting to generalize these analyses to capture features that are non-local and that persist in time.