2 Governing equations

The flow is governed by the incompressible Navier–Stokes equations:

Decomposing the velocity and pressure of the flow into components corresponding to the base flow and to a small perturbation, $\boldsymbol{u}=\boldsymbol{U}+{\it\epsilon}\boldsymbol{u}^{\prime }$ , $p=P+{\it\epsilon}p^{\prime }$ with ${\it\epsilon}\ll 1$ and linearising the equations around the base flow yields the linearised Navier–Stokes equations given here at first order of ${\it\epsilon}$ :

The base flow of the simulation is composed of a pair of counter-rotating Lamb–Oseen vortices, with initial azimuthal velocity $U_{{\it\theta}}$ defined as

with ${\it\Gamma}$ the circulation of the vortex.

All results are rendered non-dimensional by the initial distance separating the vortices $b$ and the drift velocity of the vortex pair through mutual induction $U_{drift}={\it\Gamma}/2{\rm\pi}b$ . The characteristic evolution time of the dipole is the time it takes to descend a distance equal to $b$ through mutual induction, ${\it\tau}=b/U_{drift}=2{\rm\pi}b^{2}/{\it\Gamma}$ . The resulting Reynolds number is $\mathit{Re}={\it\Gamma}/2{\rm\pi}{\it\nu}$ with ${\it\nu}$ the kinematic viscosity of the fluid. In this case, the Reynolds number is taken to be $\mathit{Re}=1000$ and the initial aspect ratio of the vortex pair is $a/b=0.18$ . In order to avoid the transient phase due to the adaptation of the Lamb–Oseen vortex pair to the Navier–Stokes equations, the base flow is time-stepped forward by several reference times ( $t\approx 8$ ) before embarking on the perturbation analysis. The details of this method are described by Sipp (Reference Sipp1999). During the linear simulations to obtain the linear optimal perturbation the base flow is ‘frozen’. This is valid as long as the viscous timescale is large compared with the timescale of the perturbation dynamics. The timescale of the viscous diffusion of a vortex being $t_{{\it\nu}}=2{\rm\pi}a^{2}/{\it\nu}$ , the ratio of these two timescales $t_{{\it\nu}}/{\it\tau}=(a/b)^{2}\mathit{Re}$ is indeed large.

As the wavelength of the most unstable Crow mode is around ${\it\lambda}_{z}=7b$ , the simulations are carried out on a mesh of extent $[-L_{x}/2,L_{x}/2]\times [-L_{y}/2,L_{y}/2]\times [-L_{z}/2,L_{z}/2]$ with $L_{x}=9.6b$ , $L_{y}=10.8b$ , $L_{z}=7b$ in order to calculate one complete wavelength of the instability. Cartesian coordinates $(x,y,z)$ are used throughout the study. The largest divisions in each direction are ${\rm\Delta}x_{max}/b\approx 0.025$ , ${\rm\Delta}y_{max}/b\approx 0.03$ , ${\rm\Delta}z_{max}/b\approx 0.02$ . The DNS is carried out using the incompressible Navier–Stokes solver Nek5000 (Fischer, Lottes & Kerkemeier Reference Fischer, Lottes and Kerkemeier2008), which is based on a spectral elements method. The computational domain is split into $N_{e}^{3}$ uniformly distributed elements with $N_{e}=70$ , containing $N_{GLL}=8$ Gauss–Lobatto–Legendre (GLL) points in each direction. The largest mesh interval is given by ${\rm\Delta}d_{max}={\rm\pi}L_{d}/2(N_{GLL}-1)N_{e}$ with $d\in x,y,z$ .

Although the flow studied in this paper is initially laminar, the nature of the resulting dynamics is unknown and potentially turbulent. Therefore, it is judicious to create a mesh that can handle a turbulent flow calculation were it to materialise. The mesh must be fine enough to capture the Kolmogorov dissipation scale ${\it\lambda}_{Kolm}=\mathit{Re}^{-3/4}b\approx 6\times 10^{-3}$ and the shear layer scale at the leading hyperbolic point ${\it\delta}=(2{\rm\pi}b^{2}{\it\nu}/3^{1/2}{\it\Gamma})^{1/2}\approx 2.5\times 10^{-2}$ determined here for a point vortex model. The contribution of eddies of scale $l\leqslant 5{\it\lambda}_{Kolm}$ is negligible (Pope Reference Pope2000), and therefore the current mesh, designed on this threshold, is sufficient. As confirmation a test was carried out using half the number of elements for the same domain size, which gave close to identical results, with less than $3\,\%$ difference between the values of the circulation at $t=10$ . Similarly the adequacy of the domain size was validated using a mesh twice as large, keeping the same number of elements and GLL points. Boundary conditions in the axial direction $(z)$ are periodic. In order to compensate for the descent of the vortex pair through mutual induction a vertical velocity is imposed in the domain. Therefore, the boundary conditions on the negative and positive $y$ -borders are inflow and outflow, respectively. Finally in the spanwise direction symmetrical conditions are imposed.

Farrell (Reference Farrell1988) demonstrated that the linear optimal perturbation of a flow is the adjoint of the most unstable mode of the flow. The adjoint of the most unstable Crow mode is determined by solving the adjoint linearised Navier–Stokes equations. The adjoint $\tilde{A}$ of an operator $A$ is defined as $(\boldsymbol{x},A\boldsymbol{y})=(\tilde{A}\boldsymbol{x},\boldsymbol{y})+\text{b.c.}$ with $\boldsymbol{x}$ and $\boldsymbol{y}$ vectors, $(\cdot ,\cdot )$ the scalar product and b.c. meaning boundary conditions. Applying this definition to the operator corresponding to the linearised Navier–Stokes equations gives the adjoint equations:

The initial velocity field for each DNS simulation described in this paper is expressed as $\boldsymbol{u}_{{\it\epsilon}}(0)=\boldsymbol{U}+{\it\epsilon}\Vert \boldsymbol{U}\Vert \boldsymbol{u}_{0}^{\prime }$ with $\boldsymbol{u}_{0}^{\prime }$ the normalised perturbation velocity field (the adjoint Crow mode is taken here). Once started the simulation yields the flow field $\boldsymbol{u}_{{\bf\epsilon}}(t)$ . In contrast to the linear simulations throughout which the base flow was frozen, the nonlinear simulation implicates an evolving base flow. In the linear stage of the evolution the perturbation $\boldsymbol{u}^{\prime }(t)$ can be retrieved by the following approximation $\boldsymbol{u}_{{\bf\epsilon}}(t)\simeq \boldsymbol{u}_{{\bf\epsilon}=0}(t)+{\it\epsilon}\Vert \boldsymbol{U}\Vert \boldsymbol{u}^{\prime }(t)$ . An unperturbed flow $({\it\epsilon}=0)$ is thus required in order to follow the evolution of the perturbation $\boldsymbol{u}^{\prime }(t)$ and the perturbation kinetic energy $E^{\prime }(t)=\Vert \boldsymbol{u}^{\prime }(t)\Vert ^{2}$ can then be computed as

3 Reference case: response to infinitesimal forcing of the linear optimal perturbation

The reference case for the nonlinear dynamics of the Crow instability is established for an infinitesimal initial perturbation amplitude to get as close to the linear limit as possible ( ${\it\epsilon}\rightarrow 0$ ). This case can be assimilated to the far-field wake of an aircraft in very calm and homogeneous atmospheric conditions. The linear optimal perturbation is superimposed on the base flow with an initial amplitude ${\it\epsilon}=10^{-3}$ . It is important to note that at such a small initial amplitude, the overall dynamics will be the same as if the flow was initialised with the most unstable Crow mode. The use of the linear optimal perturbation here has only one objective: accelerating the development of the instability to limit simulation time costs.

The dynamics of the vortex pair is concurrent with that described by Crow (Reference Crow1970). Isovorticity contours illustrate the dynamics in figure 1. As the simulations are carried out over one wavelength of the Crow instability, only one wavelength is shown. The vortices undergo a gradual symmetrical sinusoidal deformation with the cores remaining in a plane of angle approximately equal to $45^{\circ }$ with respect to the spanwise direction $(x)$ . Eventually at time $t=9$ the closest parts of the vortices meet and are pressed together. The descent of the dipole and the continued deformation of the vortices creates a distinctive head–tail dipole shape in the connection $z$ -planes. Vorticity cancellation between the opposite signed vortices becomes very strong in these areas, and the linking phenomenon is observed. This well-known event is described in more detail in the following paragraph. Once linking is complete, the vortices form a chain of vortex rings connected by thin vorticity threads. Due to the original symmetrical deformation of the vortices, the sides of the rings remain curved upwards, as can be observed at $t=11$ .

Figure 1. Isovorticity contours at $20\,\%$ of the initial maximum of vorticity ${\it\omega}=0.2|{\it\omega}_{max}(t=0)|$ showing the initial development of the Crow instability and the linking of the vortices: (a)  $t=0$ ; (b)  $t=7$ ; (c)  $t=8$ ; (d)  $t=9$ ; (e)  $t=10$ ; (f)  $t=11$ .

Figure 2. Evolutions of the circulation ${\it\gamma}_{z}$ in the connection plane ( $z=L_{z}/2$ ) and the circulation ${\it\gamma}_{x}$ in the dividing plane ( $x=0$ ). Evidence of the transfer of axial vorticity to spanwise vorticity throughout the linking process.

Over recent years a great deal of attention has been paid to vortex linking because of its supposed involvement in turbulence cascade and turbulent noise generation. In spite of a great number of analytical, numerical and experimental studies, some aspects of this complicated process are still unclear. Saffman (Reference Saffman1990) developed a model describing the linking of two counter-rotating vortices and Moriconi (Reference Moriconi2000) extended the model into high Reynolds numbers. Detailed descriptions of the physical mechanisms at play are given by many authors (Melander & Hussain Reference Melander and Hussain1989; Marshall, Brancher & Giovannini Reference Marshall, Brancher and Giovannini2001; Hussain & Duraisamy Reference Hussain and Duraisamy2011). To summarise quickly, when the vortices touch, viscous diffusion causes the touching parts to cancel each other out. The remaining strands on the outer sides of the vortex cores then reconnect to form bridges joining the two vortices. These bridges become strong and move away from each other through self-induction, leaving behind the remnants of the original vortices in threads. The linking of the vortex pair into vortex rings can be observed particularly well when studying the evolution of the axial circulation ${\it\gamma}_{z}$ in the connection $z$ -plane ( $z=\pm L_{z}/2$ ) and the spanwise circulation ${\it\gamma}_{x}$ in the symmetry plane ( $x=0$ ), as can be seen in figure 2. These circulations are defined by the following expressions:

with ${\it\omega}_{x}$ and ${\it\omega}_{z}$ the $x$ - and $z$ -components of the vorticity vector, respectively. Vorticity cancellation between the colliding counter-rotating vortices leads to a dramatic drop of the axial circulation ${\it\gamma}_{z}$ from $t=8$ to $t=10$ . At this point the threads linking the vortex rings are still present and so the decrease of the axial circulation slows before finally reaching almost zero at around $t=13.5$ . From this point the axial circulation is nil, demonstrating that once the vortex rings are formed and the threads have disappeared, no axial vorticity exists in the ( $z=L_{z}/2$ )-plane. The spanwise circulation ${\it\gamma}_{x}$ is negligible initially but from $t=8$ increases rapidly, eventually attaining $98\,\%$ of the initial value of ${\it\gamma}_{z}$ . The evolution of these quantities demonstrates the transfer of vorticity from the axial direction to the spanwise direction. Of the numerous ways to characterise this phenomenon, one of the most striking is the connection time $t_{rc}$ . This corresponds to the duration between the vortex cores touching and the bridges separating. From pure observation of the connection progress, it is virtually impossible to determine a definitive connection time value. However Melander & Hussain (Reference Melander and Hussain1989) and Leweke & Williamson (Reference Leweke and Williamson2011) noted that it is possible to deduce a value by inspecting the evolution of the circulation in the connection plane. Using the same method, illustrated in figure 2, the linking time in this case is found to be equal to $t_{rc}=1.5\pm 0.2$ compared with $t_{rc}=0.7\pm 0.1$ for Melander & Hussain (Reference Melander and Hussain1989) at $Re={\it\Gamma}/{\it\nu}=1000$ and $t_{rc}=0.9\pm 0.2$ for Leweke & Williamson (Reference Leweke and Williamson2011) for $Re={\it\Gamma}/{\it\nu}=2340$ . The error estimation of this type of method is around 15–20 %. The disparity between the connection time found here and those established in other papers is most certainly due to the difference in Reynolds number and the difference in axial perturbation wavelength. Although the aspect ratios used are very similar ( $a/b=0.18$ for Leweke and Williamson and $a/b=0.22$ for Melander and Hussain) the axial wavelengths chosen for the deformation of the vortex pair are quite different: ${\it\lambda}_{z}/b=5.13$ for Leweke and Williamson and ${\it\lambda}_{z}/b=3.85$ for Melander and Hussain (as opposed to ${\it\lambda}_{z}/b=7$ in the present study).

The subsequent evolution of the flow is analysed in the light of several quantities, namely the length-averaged axial circulation ${\it\Gamma}_{z}(t)$ and spanwise circulation ${\it\Gamma}_{x}(t)$ which indicate the strength of the vortices, the enstrophy ${\it\xi}(t)$ and the total kinetic energy $E(t)$ . These quantities are expressed as integrals over the computational volume  $V$ :

Figure 3(a,b) shows the evolution of these quantities in time. The evolution of the length-averaged circulations highlights the exchange of circulation between the axial and spanwise directions at the linking stage, and then again at $t=15$ and $t=24$ . After linking, the maximum of length-averaged circulation $\max ({\it\Gamma}_{x},{\it\Gamma}_{z})$ remains at around $60\,\%$ of the initial value of ${\it\Gamma}_{z}$ and is attained at $t\approx 13$ by ${\it\Gamma}_{x}$ and $t\approx 19$ by ${\it\Gamma}_{z}$ . This demonstrates the existence of a periodic behaviour that will be analysed further in the following paragraphs by observing the dynamics of the vortex rings. The nature of the flow can be examined by considering the evolution of the enstrophy ${\it\xi}$ . The progression in time of the enstrophy is governed by

with $S_{{\it\omega}{\it\omega}}=2[(\boldsymbol{n}_{{\it\omega}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}})\boldsymbol{u}]\boldsymbol{\cdot }\boldsymbol{n}_{{\it\omega}}$ the component of strain in the vorticity direction and $\boldsymbol{n}_{{\it\omega}}={\bf\omega}/\Vert {\bf\omega}\Vert$ . The left-hand side of the equation constitutes the advection–diffusion term, the first term of the right-hand side is the production term and the last term ensures viscous dissipation. Production of enstrophy occurs through stretching of vorticity. During the first stage of the evolution of the vortices, before the connection, the enstrophy decreases rapidly which can be expected in the case of a laminar flow. At $t=10$ the linking process entails massive vorticity stretching which generates a large spike of enstrophy. The subsequent enstrophy peaks are due to the stretching of smaller structures formed around the vortex ring. Finally the kinetic energy of the flow decreases steadily over time as a consequence of viscous dissipation. This behaviour can be expected from a laminar viscous flow.

Figure 3. (a) Evolution of length-averaged axial and spanwise circulations ${\it\Gamma}_{z}$ and ${\it\Gamma}_{x}$ over time. (b) Total kinetic energy $E$ and enstrophy ${\it\xi}$ evolution over time. Note the kinetic energy and enstrophy were integrated over the whole domain whereas the circulations were integrated over half of the domain.

At the end of the linking process $t=11$ , the vortex rings produced are circular and the sides corresponding to the original straight vortices are curved upwards as a result of the Crow instability (see figure 1). Self-induced vorticity due to the curvature of the vortex draws these sides outwards (in the spanwise direction) until the ring is flat. The resulting rings are then elliptic with the major axis aligned in the spanwise direction. The dynamics of an elliptic vortex ring have been studied analytically by Arms & Hama (Reference Arms and Hama1965). Using the Biot–Savart law and adopting a localised-induction theory Arms & Hama (Reference Arms and Hama1965) established several principles that regulate the dynamics of elliptic vortex rings. They concluded that elliptic vortex rings experience periodic deformations, where the major and minor axes periodically exchange their orientations. The period $T$ of this behaviour depends on the eccentricity $e=\sqrt{1-l^{2}/L^{2}}$ of the original ellipse, where $L$ is the length of the semi-major axis and $l$ the length of the semi-minor axis. The greater the eccentricity, the more complex is the deformation. These results were validated numerically and experimentally by Dhanak & Bernardinis (Reference Dhanak and De Bernardinis1981).

An estimation for the initial eccentricity in this case is $e\approx 0.97$ with an axis ratio $l/L\approx 0.26$ . The different stages of evolution of the vortex rings, illustrated in figure 4 with vorticity norm isocontours, are very similar to those described by Arms & Hama (Reference Arms and Hama1965) and Dhanak & Bernardinis (Reference Dhanak and De Bernardinis1981), in particular the ‘figure-of-eight’ shape attained at the half-period $T/2$ (here at $t=19$ ). The cores are very close in the centre of the domain at this time, but they do not touch enough for linking to occur and the deformation continues. Had the cores touched at this point, viscous diffusion would certainly have lead to the annihilation of vorticity between the parts of opposing sign and the vortices would have split once more into two smaller vortex rings. It is conceivable that another Crow wavelength, generating a different eccentricity, could lead to such a scenario. At time $T$ the resulting ring is an ellipse of similar size and shape to the original ellipse.

Figure 4. Vorticity isocontours at $14\,\%$ of the initial vorticity norm maximum ${\it\omega}=0.14|{\it\omega}_{max}(t=0)|$ for ${\it\epsilon}=10^{-3}$ , $t=13$ , 15, 16, 17, 19, 22, 23, 24, 26 (a,b) and $20\,\%$ of the initial vorticity norm maximum ${\it\omega}=0.2|{\it\omega}_{max}(t=0)|$ for ${\it\epsilon}=10^{-2}$ , $t=7$ , 9, 10, 11, 12, 13, 14, 15, 16.4 (c) and ${\it\epsilon}=3\times 10^{-2}$ , $t=7$ , 8.2, 9, 10, 12 (d) evolution of the vortex rings during the first period $T$ . View from above (a) and the front (b) for each stage of the ${\it\epsilon}=10^{-3}$ case, view from above only for ${\it\epsilon}=10^{-2}$ (c) and ${\it\epsilon}=3\times 10^{-2}$ (d).

By observing the movement of the vortex ring, the deformation period is estimated to be $T=13$ . Two complete periods are captured by the present simulation although only one is shown in figure 3. Using the expression of Arms & Hama (Reference Arms and Hama1965) derived by linearised analysis of a circular vortex ring and adapted to account for the eccentricity of the elliptic ring, the theoretical period of oscillation is $T_{th}=10.4$ , providing reasonable agreement with the simulated value. Such discrepancies between theory and simulation were also observed by Arms & Hama (Reference Arms and Hama1965): typically a $6\,\%$ difference for $e=0.94$ and larger for increased eccentricity.

This behaviour of vortex rings issued from the linking of a vortex pair was observed experimentally by Leweke & Williamson (Reference Leweke and Williamson2011). They observed the successive deformations characteristic of elliptic vortex rings over less than a period, and in some cases a second reconnection of the vortex rings into longer ellipses resembling the initial straight vortices. The simulated results of this paper confirm the initial development described by Leweke & Williamson (Reference Leweke and Williamson2011) (stages A to E in their figure 14) but the second reconnection of the vortices is not retrieved.

Figure 5. Perturbation energy growth rates for the linear evolution of the Crow instability eigenmode (dashed line) and the linear optimal perturbation (solid line), as well as those obtained for the evolution of the linear optimal perturbation initialised with ${\it\epsilon}=10^{-3}$ (triangles), ${\it\epsilon}=10^{-2}$ (squares) and ${\it\epsilon}=3\times 10^{-2}$ (circles).

These results indicate that when subjected to infinitesimal forcing, a pair of counter-rotating vortices go through multiple evolutions without losing much coherence over more than 25 characteristic times. Loss of coherence could not be reached in the simulated time frame meaning that in calm atmospheric conditions aircraft wakes may be present for durations longer than those required for vortex linking.

4 Nonlinear response to the linear optimal perturbation

In this section the nonlinear response to the linear optimal perturbation is analysed by a parametric study based on the initial perturbation amplitude ${\it\epsilon}$ . Figure 5 shows the perturbation energy growth $G(t)$ over time for two example values of ${\it\epsilon}$ as well as the infinitesimal case ${\it\epsilon}=10^{-3}$ and the equivalent growths for the Crow eigenmode and linear optimal perturbation as outlined by Brion et al. (Reference Brion, Sipp and Jacquin2007). The perturbation energy growth is defined as follows:

with $E^{\prime }$ the perturbation kinetic energy defined at (2.8). First note that all the perturbation energy growths initially follow the linear behaviour and then branch off, signalling at that moment the beginning of nonlinear behaviour. The recording of the perturbation energy growth is halted soon after the divergence from the linear behaviour as from that time the definition of the perturbation kinetic energy (2.8) is no longer valid. In the case of infinitesimal initial amplitude ${\it\epsilon}=10^{-3}$ this occurs at $t\approx 8$ which coincides with the start of vortex linking, the first nonlinear event to occur. Increasing the initial amplitude simply shifts this scenario in time and generates the vortex linking sooner. However, for ${\it\epsilon}=3\times 10^{-2}$ the deviation from the linear behaviour occurs before the exponential perturbation energy growth characteristic of the Crow instability is even established. Therefore, a critical amplitude ${\it\epsilon}_{c}$ which distinguishes two behaviours can be defined such that $10^{-2}<{\it\epsilon}_{c}<3\times 10^{-2}$ . For ${\it\epsilon}<{\it\epsilon}_{c}$ the Crow instability is reached before nonlinear effects come into play, whereas for ${\it\epsilon}>{\it\epsilon}_{c}$ nonlinear dynamics start during the initial transient growth due to the linear optimal perturbation.

Figure 6. (a) Kinetic energy spectrum with ${\it\epsilon}=10^{-3}$ (solid line) and ${\it\epsilon}=3\times 10^{-2}$ (dashed line) at $t=3$ ; (b) spectrum width at $t=1$ (light grey), $t=2$ (grey) and $t=3$ (black) for different initial amplitudes ${\it\epsilon}$ of the linear optimal perturbation.

An efficient way to measure the impact of nonlinear phenomena in an initially purely harmonic flow is to look at its wavelength composition. In the case of an infinitesimal perturbation, only the wavelength corresponding to the perturbation will experience growth over time. This is demonstrated in figure 6(a) which gives the volume-averaged energy spectrum of the ${\it\epsilon}=10^{-3}$ and ${\it\epsilon}=3\times 10^{-2}$ cases at $t=3$ . For ${\it\epsilon}=10^{-3}$ the majority of the energy is concentrated at wavenumber $kb=0.9$ (the first peak in the spectrum) which corresponds to the wavelength of the initial perturbation (and that of the computational domain). However, as soon as the perturbation reaches sufficient amplitude, interactions of the perturbation with itself produce smaller scales, thus widening the spectrum as can be seen for ${\it\epsilon}=3\times 10^{-2}$ . Here the Crow wavenumber $kb=0.9$ is less amplified compared with the infinitesimal case. A gain of around $100$ is obtained for ${\it\epsilon}=3\times 10^{-2}$ as opposed to about $2700$ for ${\it\epsilon}=10^{-3}$ as can also be seen in figure 5. Note that such important gains, which are much larger than those permitted by the Crow amplification, are made possible by the use of the adjoint mode as the initial perturbation.

The impact of the initial amplitude ${\it\epsilon}$ is analysed with the width of the energy spectrum ${\rm\Delta}k$ in the axial ( $z$ ) direction, defined as

with $E(k)$ the total kinetic energy of the flow at the wavenumber $k$ . Figure 6(b) shows the resulting spectrum widths for varying ${\it\epsilon}$ for $t=1$ –3. Only a few values of ${\it\epsilon}$ have been tested but they are sufficient to draw the conclusion that over a certain limit the nonlinear behaviour induced by the linear optimal perturbation is drastically different and leads to a rapid development of smaller-scale structures. The instant $t=3$ is particularly interesting because at that time (see figure 5) either exponential Crow growth has been reached ( ${\it\epsilon}<{\it\epsilon}_{c}$ ) or nonlinear behaviour has taken place ( ${\it\epsilon}>{\it\epsilon}_{c}$ ).

Next, the effect of the initial amplitude on the length-averaged axial circulation, the enstrophy and kinetic energy is evaluated. Figure 7 gives the evolutions of these quantities through time for ${\it\epsilon}=10^{-2}$ and ${\it\epsilon}=3\times 10^{-2}$ as well as ${\it\epsilon}=10^{-3}$ for comparison.

The evolution of the length-averaged axial circulation ${\it\Gamma}_{z}$ of the two simulations (see figure 7 a) shows that raising the initial perturbation amplitude accelerates the displacement of the vortices and anticipates linking. For ${\it\epsilon}=10^{-2}$ the same periodic evolution of the circulation as in the infinitesimal case can be observed after the drop due to connection. The vorticity isocontours of figure 4 confirm this behaviour. The elliptic ring shape is attained at around $t=7$ and the deformation of the ring is the same as in the infinitesimal case. However there is a major difference at later times: after the ‘figure-of-eight’ pattern the ring slowly separates into two smaller rings. This case was also described by Dhanak & Bernardinis (Reference Dhanak and De Bernardinis1981) for rings of axis ratio $l/L$ under $0.2$ . The axis ratio of the elliptic ring at $t=7$ in this case is $l/L=0.195$ . The separation of the ring into two smaller rings occurs in a very similar manner to the first connection of the vortices. Such a doubling of the dominant longitudinal wavelength has been observed in previous works (Misaka et al. Reference Misaka, Holzpfel, Henningson, Gerz, Manhart and Schwertfirm2012, for example). For ${\it\epsilon}=3\times 10^{-2}$ the periodic behaviour stops soon after the first exchange between ${\it\gamma}_{z}$ and ${\it\gamma}_{x}$ . This occurs after the first growth of ${\it\Gamma}_{z}$ at $t=3$ in figure 7(a). ${\it\Gamma}_{z}$ then slowly drops to around zero at $t\approx 9$ and eventually becomes negative, implying that the periodic ring state is halted and a completely different state is in progress. In fact, analysis of figure 4 shows that once the first elliptic ring shape is reached (at $t=7$ ) the deformation of the ring is stopped and it slowly loses coherence and breaks up under the effect of the many small-scale structures that have developed around the ring over time.

Figure 7. (a) Length-averaged axial circulation ${\it\Gamma}_{z}(t)$ , (b) enstrophy ${\it\xi}(t)$ and (c) total kinetic energy $E(t)$ for ${\it\epsilon}=10^{-2}$ (dashed lines), ${\it\epsilon}=3\times 10^{-2}$ (full lines) and ${\it\epsilon}=10^{-3}$ (dotted lines) for comparison.

Initially, the enstrophy in both simulations decreases. As before, considerable growth (see mark $P_{0}$ in figure 7 b) occurs at times coinciding with the collision of the vortices, which involves massive vorticity stretching on the centreplane separating the two cores. Subsequently for ${\it\epsilon}=10^{-2}$ the enstrophy decreases much more slowly and finally reaches the reference state ${\it\epsilon}=10^{-3}$ , meaning that no spectacular event occurs. However for ${\it\epsilon}=3\times 10^{-2}$ several peaks in enstrophy ( $P_{0}$ – $P_{3}$ ) occur after linking. Moreover, the final mean level of enstrophy remains much higher. This can be related to the persistence of small scales in the flow. The first peak ( $P_{1}$ ), at $t\approx 3.8$ corresponds to the stretching of vorticity in the later stages of the linking phase, notably within the bridges and induced by the bridges on the remaining threads. This strongly nonlinear amplification has been observed by Melander & Hussain (Reference Melander and Hussain1989), Marshall et al. (Reference Marshall, Brancher and Giovannini2001) and Hussain & Duraisamy (Reference Hussain and Duraisamy2011). Two other bumps ( $P_{2}$ and $P_{3}$ ) occur at $t\approx 6.7$ and $t\approx 10$ which were absent in the low ${\it\epsilon}$ cases and which correspond to the encounter of opposite parts of the ring. Figure 7(c) gives the evolution of the kinetic energy of the flow. Up to $t\approx 2$ both cases see the kinetic energy decreasing at the same rate, which is about that of the linear perturbation ${\it\epsilon}=10^{-3}$ . However, linking in the case ${\it\epsilon}=3\times 10^{-2}$ promotes an increased decay which is not observed in the case ${\it\epsilon}=10^{-2}$ . This increased dissipation stems from the larger production of smaller scales.

5 Accelerated decay

The results presented in § 4 have demonstrated that for initial amplitudes greater than ${\it\epsilon}_{c}$ , the evolution of the flow is significantly changed compared with the reference evolution described in § 3. Most importantly the flow rapidly loses most of its coherence. The underlying physical mechanisms are described in the following paragraphs. Note that the formation of a pair of vortex rings from the initially single elliptic ring observed in the case ${\it\epsilon}=10^{-2}$ and the subsequent evolution which, from our simulation results, remains strongly coherent, is not discussed in the present study.

Figure 8. Initial vorticity contour evolution for (a–c) ${\it\epsilon}=10^{-3}$ at $0.2\,\%\,|{\it\omega}_{max}(t=0)|$ (§ 3) and (d–f) ${\it\epsilon}=3\times 10^{-2}$ at $10\,\%\,|{\it\omega}_{max}(t=0)|$ : (a,d) $t=0$ , (b,e) $t=0.4$ , (c,f) $t=0.85$ .

According to figure 5 the evolution over the first characteristic time is linear, regardless of the perturbation amplitude. This is confirmed by comparing the development of the vorticity contours of figure 8 in the linear (a) and nonlinear (b) cases up to $t\approx 0.4$ . Note that the contour level is chosen to be much lower for the ${\it\epsilon}=10^{-3}$ case so that the perturbation can be seen. The perturbation sheets are amplified at the leading hyperbolic point and stretched horizontally beneath the cores, leading to the beginning of the vortex displacement. Details of this process are given by Brion et al. (Reference Brion, Sipp and Jacquin2007). At $t=0.85$ in the ${\it\epsilon}=10^{-3}$ case (§ 3), the perturbation rapidly disappears whereas for ${\it\epsilon}=3\times 10^{-2}$ the initial perturbation energy is sufficient to resist viscous dissipation and becomes further involved in the vortex pair dynamics.

By definition, the perturbation is modulated at the Crow wavelength and therefore the vorticity sheets change sign every half-wavelength. The subsequent progression can be split into two domains: one in which the perturbation vorticity is of the same sign as the vortex in the corresponding half-domain ( $x>0$ or $x<0$ ), and one in which the signs are opposite. When stretched at the leading hyperbolic point, the sheets become more compact and form vortex-like structures (see figure 8, $t=0.85$ ) below the primary vortices. From this point, the flow can be viewed locally as a four-vortex type configuration, counter-rotating in $z\in [-L_{z}/4;L_{z}/4]$ and co-rotating in $z\in [-L_{z}/2;-L_{z}/4[\,\cup \,]L_{z}/4;L_{z}/2]$ .

Figure 9. An $(x,y)$ -plane slice of the axial vorticity at $z=0$ at times $t=0$ (a), $0.85$ (b), $1.5$ (c), $1.9$ (d) and $2.4$ (e). Here $+$ gives the position of the leading hyperbolic point.

Figure 9 shows how axial vorticity evolves with time in a transverse $(x,y)$ -plane at $z=0$ , the centre of the counter-rotating section. The secondary vortices are advected along the lower oval streamline and eventually tear away from the sheets. They are then drawn closer to the original vortices by 3D self-induction due to their horseshoe shape as shown in figure 10 at $t=1.5$ . The proximity between each original vortex and its opposite-signed secondary vortex encourages local dynamics on each side, that are added to the interaction between the original vortices. Mutual-induction between the newly formed dipoles pulls the original vortices upwards and outwards, exaggerating the sinusoidal deformation of the vortex cores as can be seen at $t=1.5$ . When the secondary vortices reach the top of the dipole, they are dragged into the downwash between the cores, and recommence their advection around the vortices ( $t=2.4$ ). Let ${\it\Gamma}_{2}$ and $b_{2}$ be the circulation and separation of the secondary vortices and ${\it\Gamma}_{1}$ and $b_{1}$ those of the primary vortices at $z=0$ . The ratio ${\it\Gamma}_{2}/{\it\Gamma}_{1}$ remains approximately constant in time such that ${\it\Gamma}_{2}/{\it\Gamma}_{1}\approx -0.3$ , and the distance ratio $b_{2}/b_{1}\approx 1$ . This matches the observation of Fabre, Jacquin & Loof (Reference Fabre, Jacquin and Loof2002) for a counter-rotating vortex pair: the system is periodic, the vortices orbit around the vorticity centroids.

Figure 10. Vorticity isocontours at a value of $20\,\%\,|{\it\omega}_{max}(t=0)|$ at (a,d)  $t=1.5$ , (b,e)  $t=2.4$ and (c,f)  $t=3.4$ for ${\it\epsilon}=3\times 10^{-2}$ .

As can be observed in figure 10, at $t=3.4$ the vortex cores present visible short-wave deformations. Those in the centre of the domain $(z/{\it\lambda}_{z}\in [-0.2,0.2])$ result from the induction caused by the nearby secondary structures. To evaluate the scale of these deformations the axial evolution of the dispersion radius ${\it\delta}a(z)=(a(z)-\bar{a})/\bar{a}$ in this central region has been evaluated on each $(x,y)$ -plane $S(z)$ , with $a(z)$ calculated by

with $r=((x-x_{c})^{2}+(y-y_{c})^{2})^{1/2}$ and $x_{c}$ and $y_{c}$ the coordinates of the vortex core centre given by

$\bar{a}$ is the average dispersion radius over the section studied:

External noise is filtered out by imposing a minimum vorticity of 10 % of the maximum vorticity in each section $S(z)$ in the calculations. Figure 11 gives the Fourier transform $\hat{{\it\delta}a}(k)$ of the variation of dispersion radius ${\it\delta}a(z)$ as a function of the non-dimensional wavenumber $k\bar{a}$ . Several spikes can be seen, in particular for $k\bar{a}\approx 1.3$ which lies in the range of the optimal perturbation of an isolated vortex (Antkowiak & Brancher Reference Antkowiak and Brancher2004). This suggests that the accelerated decay is additionally promoted by an optimal perturbation of the vortex cores.

Figure 11. Fourier transform of the dispersion radius variation $\hat{{\it\delta}a}$ as a function of non-dimensional axial wavenumber $k\bar{a}$ at $t=3.4$ .

Figure 12 gives the axial vorticity contours in the $z=L_{z}/2$ -plane, the centre of the co-rotating section, which corresponds with the foreground plane of figure 10. The secondary vortices form a second counter-rotating dipole which moves downwards under the combined effects of the velocity field imposed by the primary vortices and of the self-advection between the perturbation vorticity. The separation distance ratio is $b_{2}/b_{1}\ll 1$ and the behaviour is divergent which concurs with the aforementioned study of Fabre et al. (Reference Fabre, Jacquin and Loof2002) except that here the perturbation vorticity diffuses before escaping the primary vortices. In return the secondary vortices induce a displacement of the primary structures towards one another, especially at instigation. This results in an increased 3D deformation of the vortices: the part of the vortices at $z=\pm L_{z}/2$ is drawn to the symmetry plane and downwards while, as mentioned earlier, the part at $z=0$ is pulled away from the symmetry plane and upwards. This, in addition to the natural deformation of the vortices under the $m=1$ Kelvin mode at the Crow wavelength, provokes the accelerated linking of the vortices. At $z=L_{z}/2$ the original cores are then pressed tightly together in a head–tail shape, ready for the connection of the vortices into rings at $t=2.4$ . The evolution of the location of the leading hyperbolic point, shown in figures 9 and 12, also shows that the entire oval surrounding the recirculating flow around the vortices is modified in this process.

Figure 12. An $(x,y)$ -plane slice of the axial vorticity at $z=L_{z}/2$ at times $t=0$ (a), $0.85$ (b), $1.5$ (c), $1.9$ (d) and $2.4$ (e). Here $+$ gives the position of the leading hyperbolic point.

In addition to this analysis of the axial vorticity in two sectional planes, the map of vorticity contours in figure 10 demonstrates the appearance, in the long term, of smaller-scale variations in the flow. They are the remnants of the initial perturbation which were not present for lower values of ${\it\epsilon}$ . Observation of figures 10 and 4 shows that this added noise deteriorates into smaller scales through tilting and stretching, continuously forcing the primary vorticity, until the vortex ring finally breaks, at $t=12$ .