2.1. Outline

Figure 1. Geometry of the hub vortex and the three satellite vortices.

As an initial state, we consider a central hub Lamb–Oseen vortex surrounded symmetrically by three satellite Lamb–Oseen vortices in an incompressible, viscous fluid, where the satellite vortices provide triangular straining to the hub vortex, similar to the illustration in figure 1. The satellite vortices are positioned at the vertices of an equilateral triangle, with the hub vortex at its centre. If $\Gamma$ denotes the circulation of the hub vortex, then the circulation of each satellite vortex is chosen to be $-\Gamma$ . This choice ensures that the velocity induced by the other vortices at the centre of any vortex vanishes. As a result, the vortices are not expected to move. In this manner, the system composed of hub and satellite vortices can form a quasi-steady state, the stability of which can be analysed.

2.2. Numerical methods

The numerical analysis will involve two main steps. First, we obtain the quasi-steady base flow by solving the two-dimensional (2-D) Navier–Stokes equations; then we integrate the linearised Navier–Stokes equations around this base flow over a sufficiently long duration to determine the most unstable mode, if it exists, for a linear disturbance with a specified axial wavenumber, as detailed in the following paragraphs.

The analysis is performed in cylindrical coordinates. The base flow is computed in two dimensions $(r, \theta )$ , while the linear stability analysis extends to three dimensions $(r, \theta , z)$ . The linearised Navier–Stokes equations for the disturbances $u, v, w$ and $p$ around a given 2-D base flow $(U, V)$ can be expressed as

(2.1) \begin{align} {\frac {\partial {u}}{\partial {t}}} &= -U{\frac {\partial {u}}{\partial {r}}} -\frac {V}{r}{\frac {\partial {u}}{\partial {\theta }}} + \frac {2Vv}{r} -u{\frac {\partial {U}}{\partial {r}}} -\frac {v}{r}{\frac {\partial {U}}{\partial {\theta }}} \nonumber \\[6pt] & \quad-{\frac {\partial {p}}{\partial {r}}} + \frac {1}{{Re}} \left [\left (\nabla ^2 -\frac {1}{r^2}\right ) u - \frac {2}{r^2} {\frac {\partial {v}}{\partial {\theta }}}\right ], \end{align}

(2.2) \begin{align} {\frac {\partial {v}}{\partial {t}}} &= -U{\frac {\partial {v}}{\partial {r}}} -\frac {V}{r}{\frac {\partial {v}}{\partial {\theta }}} -u{\frac {\partial {V}}{\partial {r}}} -\frac {v}{r}{\frac {\partial {V}}{\partial {\theta }}} - \frac {Uv +uV}{r} \nonumber \\[6pt] & \quad-\frac {1}{r}{\frac {\partial {p}}{\partial {\theta }}} + \frac {1}{{Re}} \left [\left (\nabla ^2 -\frac {1}{r^2}\right ) v + \frac {2}{r^2} {\frac {\partial {u}}{\partial {\theta }}}\right ], \end{align}

(2.3) \begin{align} {\frac {\partial {w}}{\partial {t}}} &= -U{\frac {\partial {w}}{\partial {r}}} -\frac {V}{r}{\frac {\partial {w}}{\partial {\theta }}} -{\frac {\partial {p}}{\partial {z}}} + \frac {1}{{Re}}\, \nabla ^2 w, \end{align}

where $\nabla ^2 = \partial _r^2 + (1/r)\,\partial _r + (1/r^2)\, \partial ^2_{\theta } + \partial ^2_{z}$ . These equations are discretised and solved numerically on the same 2-D grid as the base flow. By assuming periodic boundary conditions in the $z$ -direction, the above equations are separable in the $z$ -direction, which allows us to consider a single wave of the form

(2.4) \begin{eqnarray} F(t, r, \theta , z) = {\textrm e}^{{\textrm {i}} k z} \sum _{m} \hat {F}_{m} (t,r)\,{\textrm e}^{{\textrm {i}} m\theta } \end{eqnarray}

for $u, v, w$ and $p$ . The 2-D numerical domain spans $0 \leqslant r \leqslant L_r$ and $0 \leqslant \theta \leqslant 2\pi$ , with $L_r = 1000$ . For spatial discretisation, we use a sixth-order-accurate compact scheme (Lele Reference Lele1992) in $r$ , and a Fourier spectral method in $\theta$ , with periodic boundary conditions. To avoid the singularity at $r = 0$ , we expand the $r$ axis to $-L_r \leqslant r \leqslant L_r$ , and omit a grid point at $r = 0$ . The Poisson equation for the pressure is decomposed into a set of ordinary differential equations for individual Fourier modes, which are also solved with the sixth-order compact scheme. For temporal discretisation, we apply the Crank–Nicolson scheme to the viscous terms, and the second-order Adams–Bashforth method to the other terms. Further details are provided in Appendix A of Hattori et al. (Reference Hattori, Blanco-Rodríguez and Le Dizès2019).

The base flow is obtained once the vortices have equilibrated in the field of the other vortices. This rapid relaxation process has been described in the literature for two vortices (Sipp, Jacquin & Cosssu Reference Sipp, Jacquin and Cossu2000; Le Dizès & Verga Reference Le Dizès and Verga2002). After equilibrium is reached, the 2-D solution evolves on a slow viscous time scale. This evolution will not be considered further, as we assume a frozen base flow and linearise the Navier–Stokes equations around it. The same numerical methods and discretisation used for obtaining the base flow are applied to integrate the frozen linearised Navier–Stokes equations.

2.3. Initial conditions and simulation parameters

The vortices are initially positioned as shown in figure 1. The hub vortex is located at the centre of the frame, and the three satellite vortices are placed at distance of $R$ from the hub vortex centre, with angular separation $2\pi /3$ . Initially, all vortices are assumed to have a Gaussian vorticity profile with core radius $a_i$ . After the relaxation process, the core size of the vortices has slightly increased to a larger value $a$ . This new core size will be used to non-dimensionalise all spatial variables in the theory. The vortex circulation, however, is conserved. Thus the Reynolds number ${Re}$ defined by

(2.5) \begin{equation} {Re} =\frac {\Gamma }{2 \pi \nu }, \end{equation}

where $\nu$ is the kinematic viscosity, does not change. We start the simulation with core size $a_i = 0.8$ for both hub and satellite vortices, and stop it when the core size of the hub vortex reaches $a = 1$ , which occurs well after the establishment of a quasi-steady state. The distance between the hub and the satellite vortices, denoted by $R$ , is set to $R = 5$ . Unlike the growing core size, this distance remains nearly constant as the base flow reaches a quasi-steady state. Thus at the quasi-steady state, we have $a = 1$ and $R = 5$ . The initial circulation $\Gamma$ of the hub vortex is set to $2\pi$ ( $-\Gamma$ for the satellite vortices). We tracked the circulation, first-order and second-order moments of the vorticity field of the hub vortex throughout the simulation, which allowed us to compute its core radius and keep track of how it changes with time. The details of these calculations are provided in Appendix A.

In the simulations, we select ${Re} = 10^3$ to obtain the base flow, but choose a larger value ${Re} = 10^4$ for the perturbation analysis. This large value of the Reynolds number for the perturbation analysis will guarantee the existence of unstable modes. The choice of a smaller value of the Reynolds number for the base flow is just to reach the quasi-steady state of core size $a=1$ from $a_i=0.8$ more rapidly. This is justified because the quasi-steady state is not expected to depend on the Reynolds number (Le Dizès & Verga Reference Le Dizès and Verga2002). The most unstable mode for a given axial wavenumber $k$ , if it exists, is obtained by DNS from a random vorticity initial condition. The irrelevant modes decay, and the unstable modes grow exponentially, from which the growth rate can be obtained. A non-uniform grid is used where the number of grid points in $r$ is taken to be $695$ , and the number of Fourier modes in $\theta$ to be $512$ . Radial grid spacing $0.0475$ is used in the region of the hub vortex and satellite vortices.

4.1. Perturbation equations of linear disturbances

If one uses expressions (3.5) and (3.6) for the base flow, then the linearised Navier–Stokes equations for the velocity–pressure field $\boldsymbol{u} = (u,v,w,p)$ of the perturbations can be written in the form

(4.1) \begin{equation} \left ( \mathcal{L}\frac {\partial }{\partial t} + \mathcal{P}\frac {\partial }{\partial z} + \mathcal{M} - \frac {\mathcal{V}}{{Re}} \right )\boldsymbol{u} = \epsilon \left ( {\textrm{e}}^{3{\textrm{i}}\theta }\,\mathcal{N} + {\textrm{e}}^{-3{\textrm{i}}\theta }\,\overline {\mathcal{N}} \right )\boldsymbol{u}, \end{equation}

where the matrices $\mathcal{L}$ , $\mathcal{P}$ , $\mathcal{M}$ , $\mathcal{V}$ , $\mathcal{N}$ are given in Appendix B.

4.2. Triangular instability resonance mechanism

The mechanism for the growth of perturbations in a triangular-strained vortex is the same as that for the elliptic instability in the quadripolar-strained vortex (Moore & Saffman Reference Moore and Saffman1975; Tsai & Widnall Reference Tsai and Widnall1976). The only difference is that the triangular strain field corresponds to an azimuthal wavenumber $m=3$ perturbation of the vortex, while the elliptic strain field corresponds to an $m=2$ perturbation.

As for the elliptic instability, it is associated with a phenomenon of resonance of two quasi-neutral waves of the underlying vortex with the strain field. These waves, also called Kelvin modes, have an expression of the form

(4.2) \begin{equation} \boldsymbol{u} = \tilde {\boldsymbol{u}}(r)\,{\textrm{e}}^{{\textrm{i}}(m\theta + kz - \omega t)}, \end{equation}

where $m$ is the azimuthal wavenumber, $k$ is the axial wavenumber, $\omega$ is the frequency, and $\tilde {\boldsymbol{u}}(r)$ is the eigenfunction field of the Kelvin mode. They satisfy the perturbation equations for the unstrained vortex

(4.3) \begin{equation} \left ( \omega \mathcal{L} - k\mathcal{P} + {\textrm{i}}\,\mathcal{M}(m)-\frac {{\textrm{i}}}{{Re}}\,\mathcal{V}(m,k)\right )\tilde {\boldsymbol{u}} = 0, \end{equation}

where $\mathcal{M}(m)$ and $\mathcal{V}(m,k)$ are obtained by replacing $\partial /\partial \theta$ with ${i}m$ and $\partial /\partial z$ with ${i}k$ , in the matrix operators $\mathcal{M}$ and $\mathcal{V}$ , given in Appendix B.

To get the possible coupling of two Kelvin modes with the strain field, one should find two neutral Kelvin modes $(\omega _A, k_A, m_A)$ and $(\omega _B, k_B, m_B)$ satisfying

(4.4) \begin{equation} \omega _A = \omega _B,\;\;\; k_A = k_B,\;\;\;|m_A-m_B|=3. \end{equation}

Kelvin modes of the Lamb–Oseen vortex have been analysed by Fabre et al. (Reference Fabre, Sipp and Jacquin2006). In a viscous fluid, Kelvin modes are always damped ( $\textrm {Im}(\omega ) \lt 0$ ). Some of the Kelvin modes become neutral as the Reynolds number goes to infinity. But others continue to exhibit a large damping, which is associated with the presence of a critical layer (Sipp & Jacquin Reference Sipp and Jacquin2003). This feature makes the Lamb–Oseen vortex very different from the Rankine vortex, for which there is no critical layer damping.

4.2.1. Large- $k$ asymptotic prediction of the Kelvin modes

Le Dizès & Lacaze (Reference Le Dizès and Verga2005) have developed an asymptotic theory using WKB analysis to describe the various types of Kelvin modes that can be obtained in the infinite Reynolds number and large axial wavenumber limit. For the Lamb–Oseen vortex, they found two types of neutral modes: regular neutral core modes where no critical point is present, and singular neutral core modes for which a critical point is present on the real axis but the critical layer damping is asymptotically small. They showed that the frequency range of each mode can be obtained by analysing the three functions $\omega _\pm (r)$ and $\omega _c(r)$ defined by

(4.5) \begin{align}& \omega _{\pm }(r) = m\Omega _0(r) \pm \sqrt {2\Omega _0(r)\zeta _0(r)}, \end{align}

(4.6) \begin{align}&\qquad\qquad \omega _c(r) = m\Omega _0(r). \end{align}

The functions $\omega _\pm (r)$ are known as the epicyclic frequencies. They define the upper and lower bounds of the frequency range within which regular neutral core modes (blue regions in figure 7) can exist at a given radial location $r$ . Within this frequency range, the modes exhibit oscillatory behaviour for radial positions less than $r$ , and decay exponentially beyond that point. The function $\omega _c(r)$ , referred to as the critical frequency curve, gives the frequency at which a critical point occurs at $r$ . Regular modes cannot exist at these critical frequencies. However, the damping caused by the critical layer can become asymptotically small if the critical point for a given mode lies at a very large value of $r$ , far from the region where the mode behaves neutrally. Such modes are called singular neutral core modes (red regions in figure 7). For each azimuthal wavenumber $m$ , one can determine the corresponding frequency intervals in which these behaviours occur. For a more detailed explanation, readers are referred to Le Dizès & Lacaze (Reference Le Dizès and Lacaze2005). The condition of resonance can therefore be analysed by looking at the possible overlap of these frequency intervals for a couple of azimuthal wavenumbers $m$ and $m+3$ . This analysis leads to a unique possibility: only the couple $(m_A,m_B)=(1,- 2)$ (and $(m_A,m_B)=(- 1,2)$ , by symmetry) exhibits a frequency overlap. It corresponds to the frequency interval $ - 0.387 \lt \omega \lt 0$ , where $m_A=- 2$ singular core modes and $m_B=1$ regular core modes both exist, as shown in figure 7.

Figure 7. Plots of the epicyclic frequencies $\omega _+$ , $\omega _-$ (solid lines) and the critical frequency $\omega _c$ (dashed line) as functions of the radial coordinate $r$ , for (a) $m=1$ , (b) $m=- 2$ . In each plot, the blue regions (resp. red regions) indicate the frequency intervals of regular neutral core modes (resp. singular neutral core modes); the hatched region indicates the frequency interval where resonance between $m=1$ and $m=- 2$ modes is possible.

4.2.2. Numerical determination of the Kelvin modes

The Kelvin modes can also be obtained by numerically solving the eigenvalue problem (4.3), and we use these numerically computed modes for the remainder of the analysis. Our numerical solver employs a Chebyshev spectral collocation method, following the approach of Fabre & Jacquin (Reference Fabre and Jacquin2004). The eigenvalue problem defined in the domain $0 \lt r \lt \infty$ is extended to $- \infty \lt r \lt \infty$ and mapped onto a contour in the complex- $r$ plane. It is then solved in the Chebyshev domain $(- 1, 1)$ using $2(N+1)$ collocation points. A resolution $N = 200$ is found to be sufficient. For the inviscid limit, we adopt a complex mapping function, similar to the one used by Fabre & Jacquin (Reference Fabre and Jacquin2004), defined as

(4.7) \begin{equation} r = \frac {H\xi }{1-\xi ^2} + {\textrm{i}}\frac {A\xi }{\sqrt {1-\xi ^2}}, \end{equation}

where the parameter $H$ controls the radial spreading of the collocation points, and the parameter $A$ determines the inclination of the contour in the complex plane. We also take advantage of the parity properties of the eigenfunctions. For odd values of $m$ , we express $\tilde {w}$ and $\tilde {p}$ using odd polynomials, and $\tilde {u}$ and $\tilde {v}$ using even polynomials. Conversely, even values of $m$ , $\tilde {w}$ and $\tilde {p}$ are represented using even polynomials, while $\tilde {u}$ and $\tilde {v}$ are expressed using odd polynomials.

The above procedure has been done for the two azimuthal wavenumbers $m_A=1$ and $m_B=- 2$ , and a wavenumber $k$ varying between $0$ and $10$ at ${Re} = 10^4$ . The result gives the complex frequency curves that have been plotted in figure 8(a,b). One clearly sees in these plots the frequency interval where regular and singular core modes were expected from figure 7.

Figure 8. Dispersion curves of the Kelvin modes, obtained by solving the eigenvalue problem in (4.3) on the real axis for (a) $m_A = 1$ and (b) $m_B = -2$ at ${Re} = 10^4$ . The real part of the frequency, $\omega _r$ , is plotted against $k$ , while the damping is represented by the greyscale intensity of $-\omega _i$ , the imaginary part of the frequency. (c) Resonant Kelvin modes of the unstrained Lamb–Oseen vortex at ${Re} = 10^4$ , identified at the crossing points of the dispersion curves. Modes with positive growth rates at ${Re} = 10^4$ are circled in red, with their corresponding branch indices indicated in parentheses.

Figure 8(c) displays the crossing points of $m_A=1$ and $m_B=- 2$ branches. As the damping rate of the modes shown in this figure is small, each crossing point corresponds to a (quasi-)resonant configuration. The resonant configurations marked with red circles exhibited positive growth rates, while the other resonant modes showed no growth or negligible growth rates due to critical layer damping and volumic viscous damping. A detailed analysis is provided in the subsequent subsections. We label the branches of the dispersion curves as $(l_A,l_B)$ where $l_A$ and $l_B$ are the branch labels of $m_A=1$ and $m_B=- 2$ modes, respectively.

An interesting observation is that the branches of the dispersion curves for $m = - 2$ involved in the resonance correspond to the ‘L branches’ described in Fabre et al. (Reference Fabre, Sipp and Jacquin2006) (see their figures 19 and 20). In particular, the first branch in figure 8(c) for $m = - 2$ matches the ‘F branch’, also known as the flattening wave. Modes on this branch with axial wavenumber $k \lt 3.45$ are reported to experience significant critical layer damping, while modes with larger $k$ behave as regular, weakly damped waves. This trend is also evident here: the first mode, labelled $(1,1)$ , shows stronger critical layer damping compared to the modes labelled $(2,1)$ and $(3,1)$ , which are only weakly damped. The F branch also has an analogue in the Rankine vortex, where it is referred to as the ‘isolated branch’ by Eloy & Le Dizès (Reference Eloy and Le Dizès2001) and Fukumoto (Reference Fukumoto2003). The branches below the first one are referred to as ‘L2 branches’. In our results, the modes labelled $(3,2)$ and $(4,2)$ lie in the weakly damped regions of these L2 branches. Hence among the resonant modes considered, only $(1,1)$ exhibits strong critical layer damping (also see table 1).

Table 1. Values of the parameters of the growth rate equation (4.18) for the resonance of two Kelvin modes of azimuthal wavenumbers $m_A=1$ and $m_B=- 2$ at the different resonant points. Integration is done in the complex plane as explained in the text.

4.3. Theoretical expression for the triangular instability growth rate

The method for computing the growth rate associated with the resonant coupling of two Kelvin modes with the triangular straining field is the same as that used for the elliptic instability. It is based on an asymptotic analysis in the limit of small $\epsilon$ . For more details, we refer the reader to Moore & Saffman (Reference Moore and Saffman1975).

The idea is to consider the perturbation as a combination of two normal modes of azimuthal wavenumbers $m_A$ and $m_B$ (related by (4.4)) that corresponds at leading order in $\epsilon$ to a resonant configuration of Kelvin modes. We then write

(4.8) \begin{equation} \boldsymbol{u} \sim \left (A\tilde{\boldsymbol {u}}_A(r)\textrm{e}^{{\textrm{i}}(m_A\theta )} + B\tilde{\boldsymbol {u}}_B(r)\textrm{e}^{{\textrm{i}}(m_B\theta )}\right )\textrm{e}^{{\textrm{i}}(kz - \omega t)}, \end{equation}

where the (real) axial wavenumber $k$ and the (complex) frequency $\omega$ of the two modes are assumed to be close to a resonant point defined by an axial wavenumber $k_c$ and a real frequency $\omega _c$ (corresponding to one of the crossing points shown in figure 8 c).

Substituting (4.8) in (4.1), we get two equations for the components proportional to $e^{{i}m_A\theta }$ and $e^{{i}m_B\theta }$ , respectively:

(4.9) \begin{align} & A\left ( \omega \mathcal{L} - k\mathcal{P} + {\textrm{i}}\,\mathcal{M}(m_A)-\frac {{\textrm{i}}}{{Re}}\,\mathcal{V}(m_A,k)\right )\tilde{\boldsymbol {u}}_A = {\textrm{i}}B\epsilon\, \mathcal{N}(m_B)\,\tilde{\boldsymbol {u}}_B, \end{align}

(4.10) \begin{align}& B\left ( \omega \mathcal{L} - k\mathcal{P} + {\textrm{i}}\,\mathcal{M}(m_B)-\frac {{\textrm{i}}}{{Re}}\,\mathcal{V}(m_B,k)\right )\tilde{\boldsymbol {u}}_B = {i}A\epsilon\, \overline {\mathcal{N}}(m_A)\,\tilde{\boldsymbol {u}}_A, \end{align}

where $\mathcal{N}(m_B)$ and $\overline {\mathcal{N}}(m_A)$ are obtained by replacing $\partial /\partial \theta$ with ${i}m_B$ and ${i}m_A$ , respectively, in the $\mathcal{N}$ and $\overline {\mathcal{N}}$ operators given in Appendix B. These equations show how the two modes are coupled by the straining field that is responsible for the terms on the right-hand sides of these equations.

The equation giving the complex frequency $\omega$ is obtained from an orthogonality condition with the adjoint resonant Kelvin modes. The eigenfunctions $\tilde{\boldsymbol {u}}_A^{\dagger }$ and $\tilde{\boldsymbol {u}}_B^{\dagger }$ of these two adjoint modes are the solutions of the adjoint equation of (4.3) for $(m,k) =(m_A,k_c)$ and $(m,k)= (m_B,k_c)$ , respectively. These adjoint equations are obtained using the scalar product

(4.11) \begin{equation} \left \langle \boldsymbol{u_1},\boldsymbol{u_2} \right \rangle = \int _{0}^{\infty }\boldsymbol{u}^*_1(r)\,\boldsymbol{u}_2(r) r\,{\textrm d}r, \end{equation}

where $*$ denotes the complex conjugate. We have two options for applying the orthogonality condition: either we consider the inviscid expressions $(\tilde{\boldsymbol {u}}_A^{(\infty )},\tilde{\boldsymbol {u}}_B^{(\infty )})$ and $(\tilde{\boldsymbol {u}}_A^{\dagger (\infty )},\tilde{\boldsymbol {u}}_B^{\dagger (\infty )})$ of the Kelvin and adjoint modes, or we choose their expression for a given (large) Reynolds number.

In the first case, upon performing the scalar product of (4.9) with $\tilde{\boldsymbol {u}}_A^{\dagger (\infty )}$ and of (4.10) with $\tilde{\boldsymbol {u}}_B^{\dagger (\infty )}$ , we obtain

(4.12) \begin{align}& \left ( \omega -\omega _A^{(\infty )} -Q_A^{(\infty )}\left(k-k_c^{(\infty )}\right)-{\textrm{i}}\frac {V_A^{(\infty )}}{{Re}} \right )A = {\textrm{i}}\epsilon C_{AB}^{(\infty )}B, \end{align}

(4.13) \begin{align}& \left ( \omega -\omega _B^{(\infty )}-Q_B^{(\infty )}\left(k-k_c^{(\infty )}\right)-{\textrm{i}}\frac {V_B^{(\infty )}}{{Re}} \right )B = {\textrm{i}}\epsilon C_{BA}^{(\infty )}A, \end{align}

where the coefficients $Q_A^{(\infty )}$ , $Q_B^{(\infty )}$ , $V_A^{(\infty )}$ , $V_B^{(\infty )}$ , $C_{AB}^{(\infty )}$ , $C_{BA}^{(\infty )}$ are given by

(4.14) \begin{align}&\qquad Q_A^{(\infty )} = \frac {\left \langle \tilde{\boldsymbol {u}}_A^{\dagger (\infty )},\mathcal{P}\tilde{\boldsymbol {u}}_A^{(\infty )} \right \rangle }{\left \langle \tilde{\boldsymbol {u}}_A^{\dagger (\infty )},\mathcal{L}\tilde{\boldsymbol {u}}_A^{(\infty )} \right \rangle }, \qquad Q_B^{(\infty )} = \frac {\left \langle \tilde{\boldsymbol {u}}_B^{\dagger (\infty )},\mathcal{P}\tilde{\boldsymbol {u}}_B^{(\infty )} \right \rangle }{\left \langle \tilde{\boldsymbol {u}}_B^{\dagger (\infty )},\mathcal{L}\tilde{\boldsymbol {u}}_B^{(\infty )} \right \rangle }, \end{align}

(4.15) \begin{align}&\qquad V_A^{(\infty )}= \frac {\left \langle \tilde{\boldsymbol {u}}_A^{\dagger (\infty )},\mathcal{V}\tilde{\boldsymbol {u}}_A^{(\infty )} \right \rangle }{\left \langle \tilde{\boldsymbol {u}}_A^{\dagger (\infty )},\mathcal{L}\tilde{\boldsymbol {u}}_A^{(\infty )}\right \rangle }, \qquad V_B^{(\infty )} = \frac {\left \langle \tilde{\boldsymbol {u}}_B^{\dagger (\infty )},\mathcal{V}\tilde{\boldsymbol {u}}_B^{(\infty )} \right \rangle }{\left \langle \tilde{\boldsymbol {u}}_B^{\dagger (\infty )},\mathcal{L}\tilde{\boldsymbol {u}}_B^{(\infty )} \right \rangle }, \end{align}

(4.16) \begin{align}& C_{AB}^{(\infty )} = \frac {\left \langle \tilde{\boldsymbol {u}}_A^{\dagger },\mathcal{N}(m_B)\,\tilde{\boldsymbol {u}}_B^{(\infty )} \right \rangle }{\left \langle \tilde{\boldsymbol {u}}_A^{\dagger (\infty )},\mathcal{L}\tilde{\boldsymbol {u}}_A^{(\infty )} \right \rangle }, \qquad C_{BA}^{(\infty )} = \frac {\left \langle \tilde{\boldsymbol {u}}_B^{\dagger (\infty )},\overline {\mathcal{N}}(m_A)\,\tilde{\boldsymbol {u}}_A^{(\infty )} \right \rangle }{\left \langle \tilde{\boldsymbol {u}}_B^{\dagger (\infty )},\mathcal{L}\tilde{\boldsymbol {u}}_B^{(\infty )}\right \rangle } . \end{align}

The frequencies $\omega _A^{(\infty )}$ and $\omega _B^{(\infty )}$ are the inviscid estimates of the frequencies of the resonant Kelvin modes at $k=k_c^{(\infty )}$ . They can be expressed as

(4.17) \begin{equation} \omega _A^{(\infty )} = \omega _c^{(\infty )} + {i}\, \textrm {Im}\left(\omega _A^{(\infty )}\right) ,\quad\omega _B^{(\infty )} = \omega _c^{(\infty )} + {i}\, \textrm {Im}\left(\omega _B^{(\infty )}\right), \end{equation}

where $\textrm {Im}(\omega _A^{(\infty )} )$ and $\textrm {Im}(\omega _B^{(\infty )})$ are the critical layer damping rates of the modes.

For the azimuthal wavenumber couple $(m_A,m_B) = (1,-2)$ , only the Kelvin mode with azimuthal wavenumber $m_B=- 2$ is expected to experience critical layer damping at resonance. In particular, we have $\textrm {Im}(\omega _A^{(\infty )}) =0$ for $m_A=1$ at resonance.

Equations (4.12) and ( 4.13) with ( 4.17) give an equation for the complex frequency $\omega$ :

(4.18) \begin{align} & \left ( \omega - \omega _c^{(\infty )} - Q_A^{(\infty )}\left(k - k_c^{(\infty )}\right) - {\textrm{i}}\frac {V_A^{(\infty )}}{{Re}} \right ) \nonumber \\[2pt] &\left ( \omega - \omega _c^{(\infty )} - {\textrm{i}}\,\textrm {Im}\left(\omega _B^{(\infty )}\right) - Q_B^{(\infty )}\left(k - k_c^{(\infty )}\right) - {i}\frac {V_B^{(\infty )}}{{Re}} \right ) = -\epsilon ^2 \left(N^{(\infty )}\right)^2, \end{align}

where

(4.19) \begin{equation} N^{(\infty )} = \sqrt {C_{AB}^{(\infty )}C_{BA}^{(\infty )}}. \end{equation}

The growth rate $\sigma$ is defined as $\textrm {Im}(\omega )$ and is computed using the equation above. Both modes oscillate with a common resonant frequency given by ${{Re}}(\omega )$ . On the left-hand side of (4.18), we recognise the dispersion relation for the two resonant modes, $A$ and $B$ , close to the resonant point. More precisely, the left-hand side can be written as $(\omega - \omega _{A}(k,{Re}))(\omega -\omega _B(k,{Re}))$ , where $\omega _A(k,{Re})$ and $\omega _B(k,{Re})$ are the complex frequencies of the two Kelvin modes for a wavenumber $k$ close to $k_c$ , and a large Reynolds number. The right-hand side corresponds to the coupling term that is responsible for the instability.

The integrals in the inner products that appear in the coefficients of (4.14), (4.15) and (4.16) are evaluated along a complex path that avoids the critical point associated with the mode $m = - 2$ from below, since it moves in the lower part of the complex plane, as explained in Le Dizès (Reference Le Dizès2004). To achieve this, we use the complex mapping given in (4.7) with parameters $H = 4$ and $A = - 0.5$ . These values of $H$ and $A$ result in the contour passing slightly below the critical points for all the resonant points corresponding to $m=- 2$ . The integrals corresponding to $m = 1$ are also computed along this same contour. Since the $m = 1$ modes do not have a critical point, the values of the integrals remain unchanged. In table 1, we provide the values of the parameters appearing in (4.18) for all the resonant points that have been identified in figure 8 for the two azimuthal wavenumbers $m_A=1$ and $m_B=- 2$ .

Figure 9. Growth rates $\sigma$ are plotted against the axial wavenumber $k$ for the resonant modes, based on theoretical predictions. Solid black lines represent results from (4.18), while solid red lines correspond to (4.20), both computed for ${Re} = 10^4$ . Dashed black lines show results from (4.18) in the inviscid limit ( ${Re} = \infty$ ). The corresponding branch indices are also indicated.

Figure 10. Comparison of growth rates $\sigma$ against the axial wavenumber $k$ for the first four resonant modes computed using (4.20) (solid black lines) and DNS (circles). Values of the corresponding branch indices are reported.

An equation for the growth rate can also be obtained using the viscous expression of the resonant Kelvin modes. In that case, instead of (4.18), we get

(4.20) \begin{align}& \left ( \omega - \omega _c({Re}) - {\textrm{i}}\,\textrm {Im}(\omega _A({Re})) - Q_A({Re})(k - k_c({Re})) \right ) \nonumber\\[4pt] &\quad \times \left ( \omega - \omega _c({Re}) - {\textrm{i}}\,\textrm {Im}(\omega _B({Re})) - Q_B({Re})(k - k_c({Re})) \right ) = -\epsilon ^2 (N({Re}))^2, \end{align}

where the coefficients now depend on the Reynolds number. Note that the volumic viscous damping terms ${\textrm i}V_A^{(\infty )}/{Re}$ and ${\textrm i}V_B^{(\infty )}/{Re}$ in (4.18) are now included in the terms ${\textrm i}\,\textrm {Im}(\omega _A({Re}))$ and ${\textrm i}\,\textrm {Im}(\omega _B({Re}))$ . This second approach has also been used. For the Reynolds number ${Re}=10^4$ that we have considered, the results are similar, as can be seen in figure 9. This figure shows the five resonant modes for which a positive viscous growth rate was obtained. In this figure, we have also plotted the growth rate curves obtained from (4.18) for ${Re}=\infty$ . We clearly see that the viscous damping of the modes significantly affects the resonant configurations with the largest wavenumbers.

4.4. Numerical results

In this subsection, we first describe how growth rates and frequencies are obtained numerically, and then compare the results predicted by (4.18) or (4.20) with those obtained from DNS. The numerical results are obtained for ${Re}=10^4$ . Starting from a random initial energy distribution, modes associated with different azimuthal wavenumbers grow exponentially in time if they are unstable. The growth rate is determined as the slope (divided by 2) of the logarithm of the energy of the unstable modes. The frequency of the modes is computed by tracking their phase and measuring their rotation rate about the $z$ axis. The comparison of growth rates is shown in figure 10 for $\epsilon =0.008$ and $k$ ranging from 0.7 to 10. We observe that the first four unstable modes predicted by the theory within this wavenumber range are also captured in the DNS. The agreement in the growth rates is almost perfect. The values for the frequency and growth rate, obtained by both theory and DNS, for the four unstable modes at their resonant wavenumbers are given in table 2.

Table 2. Growth rate $\sigma$ and real part of the frequency $\omega$ of the most unstable mode for different axial wavenumbers $k$ , obtained by DNS and theory (4.20), for ${Re}=10^4$ and $\epsilon =0.008$ . The maximum growth rates predicted by (4.18) in the inviscid limit, denoted by $\sigma _{th}^{(\infty )}$ , are given in the last column. The chosen axial wavenumbers correspond to the resonant values for the Kelvin modes of azimuthal wavenumbers $(m_A=1,m_B=- 2)$ and branch labels $(l_A,l_B)$ .

In figure 11, we analyse the effect of $\epsilon$ on the growth rate curve of the most unstable mode corresponding to $(l_A,l_B)=(2,1)$ . As expected, the growth rate decreases when $\epsilon$ is halved, and the agreement between the growth rate curves is excellent, as illustrated in the figure.

Figure 11. Growth rate $\sigma$ versus $k$ for the resonant mode with $(l_A,l_B)=(2,1)$ at ${Re}= 10^4$ . Circles indicate DNS; solid black lines indicate theory; with black for $\epsilon = 0.008$ , and blue for $\epsilon = 0.004$ .

The structures and the energy ratio of the contributions from various azimuthal wavenumbers of the unstable modes are displayed in figures 12, 13, 14 and 15. In each figure, the leftmost panel shows the 2-D contours of the perturbation axial vorticity field in the hub vortex, obtained by DNS. The values are divided by $e^{\sigma t_s}$ , where $\sigma$ is the corresponding growth rate, and $t_s$ is the simulation time of the snapshot used. In the centre panels, we analyse the azimuthal composition of the unstable mode, by computing the energy percentage of each azimuthal component. This decomposition is compared to the theoretical prediction obtained for each resonant configuration. The theoretical energy ratio for mode $A$ , for instance, can be calculated by

(4.21) \begin{equation} \frac {E_A}{E_A + E_B} = \frac {\int |\tilde{\boldsymbol {u}}_A|^2\, r\,{\textrm d}r}{\int |\tilde{\boldsymbol {u}}_A|^2\, r\,{\textrm d}r + (B^2/A^2)\int |\tilde{\boldsymbol {u}}_B|^2\, r\,{\textrm d}r}, \end{equation}

Figure 12. Structure of the resonant combination with $k = 1.76$ for ${Re}= 10^4$ and $\epsilon =0.008$ , for mode $(1,1)$ . (a) Contours of the axial vorticity using DNS. (b) Energy ratio: red circles indicate DNS; blue crosses indicate theory. (c) Radial structure of the amplitude $|u(r)|$ of radial velocity: coloured lines indicate DNS; black lines indicate theory.

Figure 13. Structure of the resonant combination with $k = 5.18$ for ${Re}= 10^4$ and $\epsilon =0.008$ , for mode $(2,1)$ . (a) Contours of the axial vorticity using DNS. (b) Energy ratio: red circles indicate DNS; blue crosses indicate theory. (c) Radial structure of the amplitude $|u(r)|$ of radial velocity: coloured lines indicate DNS; black lines indicate theory.

Figure 14. Structure of the resonant combination with $k = 7.58$ for ${Re}= 10^4$ and $\epsilon =0.008$ , for mode $(3 ,1)$ . (a) Contours of the axial vorticity using DNS. (b) Energy ratio: red circles indicate DNS; blue crosses indicate theory. (c) Radial structure of the amplitude $|u(r)|$ of radial velocity: coloured lines indicate DNS; black lines indicate theory.

Figure 15. Structure of the resonant combination with $k = 9.64$ for ${Re}= 10^4$ and $\epsilon =0.008$ , for mode $(3,2)$ . (a) Contours of the axial vorticity using DNS. (b) Energy ratio: red circles indicate DNS; blue crosses indicate theory. (c) Radial structure of the amplitude $|u(r)|$ of radial velocity: coloured lines indicate DNS; black lines indicate theory.

where $|\tilde{\boldsymbol {u}}_A|^2 = \tilde {u}_A(r)^2 + \tilde {v}_A(r)^2 + \tilde {w}_A(r)^2$ , and $(B^2/A^2)$ can be calculated from (4.12) or (4.13). We can see that the unstable modes obtained from DNS are primarily composed of the azimuthal wavenumber pairs $m_A = 1$ , $m_B = - 2$ and $m_A = - 1$ , $m_B = 2$ , with negligible contributions from other azimuthal wavenumbers, as resonance is expected only between these wavenumbers. The energy ratio of the resonant azimuthal wavenumbers aligns well with theoretical predictions, except for the mode corresponding to the branch labels $(1,1)$ , which shows a slight discrepancy. It might be useful to note that the numerically obtained energy percentages of $m_A = 1$ and $m_B = - 2$ are respectively the same as those of $m_A = - 1$ and $m_B = 2$ . This is due to a certain symmetry that we enforce in the initial random vorticity distribution. Without this symmetry condition, these percentages might differ; however, the mutual energy ratio between the resonant pairs would still be the same. To further validate that the unstable resonant modes obtained from DNS share the same structure as those predicted by theory, the radial distributions of the absolute value of the radial disturbance velocity for both azimuthal wavenumbers are plotted in the right most panels of figures 12, 13, 14 and 15, comparing DNS results with theoretical predictions. The radial structures from DNS are extracted by performing a Fourier decomposition of the disturbance fields into components corresponding to different azimuthal wavenumbers. This comparison is conducted by normalising the amplitude of the $m_A=1$ azimuthal component to one at the origin in both DNS and theory. Good agreement is also observed here, except for a significant discrepancy in the amplitude of the mode $(1,1)$ corresponding to $m=- 2$ . We do not have a definitive explanation for these discrepancies in both the energy ratios and the eigenfunction radial structures for $m=- 2$ . However, we suspect that they arise from an imperfect representation of the critical layer, particularly strong for $(1,1)$ . Figure 12(c) shows additional oscillations in the theoretical predictions just beyond $r = 2$ . These oscillations result from the regularisation of the critical layer by viscous effects, which are absent in the numerical results. Therefore, we believe that the observed discrepancy is due to the inaccurate representation of the critical layer. From the leftmost and rightmost panels, it is also evident that modes associated with higher branches exhibit more zeros in the radial distributions of the disturbance fields.

4.5. Variation of maximum growth rate with ${Re}$ and $\epsilon$

Figure 16. Contours of the maximum growth rate in the $(\epsilon , {Re})$ parameter space, computed from theory using (4.18). Grey regions indicate where the mode $(1,1)$ is the most unstable, while blue regions indicate where the mode $(2,1)$ is the most unstable. Black contour lines represent growth rate levels from $0$ to $0.045$ , in steps of $0.005$ . The overall marginal stability curve, obtained from (4.22), is shown as a thick dashed black line and corresponds to the zero-growth-rate contour. The individual marginal stability curves for modes $(1,1)$ and $(2,1)$ are shown as black and blue dotted lines, respectively. Red stars mark the growth rates used in the main test cases: $\sigma = 0.0153$ for $\epsilon = 0.004$ , and $\sigma = 0.0345$ for $\epsilon = 0.008$ , at ${Re} = 10^4$ .

The theory allows us to easily analyse the effect of $\epsilon$ and ${Re}$ on the instability. Using (4.18), one can compute, for each resonant configuration, the maximum growth rate for any values of $\epsilon$ and ${Re}$ . Figure 16 shows the contours of the maximum growth rate in the $(\epsilon ,{Re})$ parameter space. For large ranges of $\epsilon$ ( $0\lt \epsilon \lt 0.01$ ) and ${Re}$ ( $100\lt {Re} \lt 10^5$ ), we find that the most unstable modes correspond to the branch label combinations $(l_A, l_B) = (1,1)$ and $(2,1)$ . Among these, the $(2,1)$ mode (shown in blue contours) is dominant over most of the parameter space, while the $(1,1)$ mode (shown in grey contours) becomes more unstable only at lower Reynolds numbers because of its small volumic viscous damping coefficients, as seen in the figure. The contour $\sigma = 0$ represents the marginal stability curve $\epsilon _s({Re})$ , which separates stable and unstable regions for these resonant modes and thus marks the onset of triangular instability. Using (4.18), this curve can be approximated as

(4.22) \begin{equation} \epsilon _s({Re}) \approx \left | \frac {\sqrt {V_A^{(\infty )} V_B^{(\infty )} + V_A^{(\infty )}\,\textrm {Im}(\omega _B^{(\infty )})\,{Re}}}{N^{(\infty )}\,{Re}} \right |, \end{equation}

with the relevant coefficients listed in table 1 for each resonant configuration. The figure also shows that for realistic values of $\epsilon$ (towards the right-hand side of the plot), triangular instability can arise at Reynolds numbers as low as ${Re} \sim 200$ . Likewise, even under weak straining conditions ( $\epsilon \lt 0.001$ ), instability may still occur at moderate Reynolds numbers ${Re} \sim 10{\,}000$ , albeit with small growth rates.

4.6. Comparison with triangular instability in the Rankine vortex

In this subsection, we briefly compare the characteristics of the triangular instability observed in the Lamb–Oseen vortex with those reported for the Rankine vortex. The latter has been studied theoretically by Eloy & Le Dizès (Reference Eloy and Le Dizès2001), and their results are used here for comparison.

The most notable difference lies in the absence of a critical layer in the Rankine vortex. This allows resonances for any pair of azimuthal wavenumbers $m$ and $m+3$ , even for asymptotically large $m$ . In contrast, for the Lamb–Oseen vortex, resonance is limited only between $m = 1$ and $m = - 2$ (or $m = - 1$ and $m = 2$ ), as critical layer damping suppresses other combinations.

For the Rankine vortex, the maximum growth rate is achieved in the large- $m$ , large- $k$ limit, with $\sigma _{max }^{(R)}= 49/32 \tilde {\epsilon } = 1.53 \tilde {\epsilon }$ , where $\tilde {\epsilon }$ is related to our external strain rate parameter $\epsilon$ by $\tilde {\epsilon } = 9/2 \epsilon$ . This gives $\sigma _{max }^{(R)} = 6.89 \epsilon$ . For the Lamb–Oseen vortex, the highest inviscid growth rate is $\sigma _{max }^{(L)} = 4.78 \epsilon$ for mode $(2,1)$ (see table 1). Alternatively, if one uses the strain rate parameter at the vortex centre, which is $\epsilon _0 = (3/2) \epsilon$ for the Rankine vortex, but $\epsilon _0 = s_0 \epsilon$ for the Lamb–Oseen vortex, we get $\sigma _{max }^{(R)} = 4.59 \epsilon _0$ and $\sigma _{max }^{(L)} = 2.70 \epsilon _0$ . In both cases, the values remain below the maximum growth rate reported for the Rankine vortex. However, focusing only on modes from the $m = 1$ and $m = - 2$ resonance, Eloy & Le Dizès (Reference Eloy and Le Dizès2001) reported $\sigma ^{(R)}= 0.95 \tilde {\epsilon } = 4.27 \epsilon = 2.85 \epsilon _0$ for the first principal mode (see their Figure 4 b). This value is comparable to the maximum inviscid growth rate found for the Lamb–Oseen vortex, which occurred for the non-principal mode $(2,1)$ . As shown in this study, the first principal mode is subject to significant critical layer damping and does not become the most unstable mode in the inviscid limit.