2.1 Viscous inner region

Consider a columnar vortex with velocity components $[0,V(r),W(r)]$ in cylindrical coordinates $(r,\unicode[STIX]{x1D703},z)$ in an incompressible and inviscid fluid. The assumption of axisymmetry of this base flow is reasonable since it describes the mean flow in the absence of strain and, thus, without non-axisymmetric effects.

Moore & Saffman (Reference Moore and Saffman1973) began by considering the inviscid roll-up of a wingtip vortex. Because the loading of the wing is an important unknown, a wing loading parameter $n$ was introduced, with $n=1/2$ corresponding to elliptical loading, $0<n<1$ . In the limit where $n\rightarrow 1$ , Batchelor’s solution is recovered. After making the light-loading approximation, the solutions of the inviscid governing equations were shown to be singular near the vortex core. This necessitates inclusion of viscosity to describe the interior of the vortex. Using boundary layer methods, Moore and Saffman then derived parabolic governing equations and obtained a solution using a similarity variable $\unicode[STIX]{x1D702}=-r^{2}/(4\unicode[STIX]{x1D708}t)$ , where $\unicode[STIX]{x1D708}$ denotes the kinematic viscosity and $t$ is time. The resulting family of similarity solutions for the azimuthal and axial mean velocities is given by

(2.1a,b ) $$\begin{eqnarray}v_{\unicode[STIX]{x1D703}}(r,t)=\frac{\unicode[STIX]{x1D6FD}}{(\unicode[STIX]{x1D708}t)^{n/2}}V_{n}(\unicode[STIX]{x1D702})\quad \text{and}\quad w(r,t)=\frac{\unicode[STIX]{x1D6FD}^{2}}{W_{\infty }(\unicode[STIX]{x1D708}t)^{n}}W_{n}(\unicode[STIX]{x1D702})\end{eqnarray}$$

respectively. The constant $\unicode[STIX]{x1D6FD}$ is related to the wing loading, angle of attack and aspect ratio, and $W_{\infty }$ is the speed of the wing at infinity. It can be seen from (2.1) that the velocity components are proportional to an azimuthal function $V_{n}(\unicode[STIX]{x1D702})$ and an axial function $W_{n}(\unicode[STIX]{x1D702})$ ; the exact expressions for both involve confluent hypergeometric functions. Moreover, the axial component $W_{n}(\unicode[STIX]{x1D702})$ is found to be the solution to a second-order ODE that must be solved numerically. Feys & Maslowe (Reference Feys and Maslowe2014) provide definitions of the functions $V_{n}(\unicode[STIX]{x1D702})$ and $W_{n}(\unicode[STIX]{x1D702})$ that were just introduced. Here, we limit ourselves to displaying the velocity components graphically.

The functions $V_{n}(r)$ and $W_{n}(r)$ are shown in figure 1, where $r$ has been non-dimensionalized by employing a length scale $r_{s}=\sqrt{4\unicode[STIX]{x1D708}t}$ ; the similarity variable $\unicode[STIX]{x1D702}$ is now related to $r$ by $\unicode[STIX]{x1D702}=-r^{2}$ . A particularly interesting aspect of the Moore–Saffman profiles is that the character of the axial component $W_{n}(r)$ changes near $n=0.44$ . For $n$ larger than this value, the vortex has a wake-like core with jet-like velocity only for large $r$ . This is consistent with Batchelor’s solution, which is recovered for $n\rightarrow 1$ (but corresponds to an unrealistic rectangular lift distribution) and which he described as being wake-like. As $n$ decreases, it can be seen that even for $n=0.7$ the axial velocity profile is not far from Gaussian. However, the magnitude of $W_{n}(0)$ decreases with $n$ , and for $n\lesssim 0.44$ , the flow relative to the wing becomes jet-like. Since experimental observations in close proximity to the wingtip have revealed both axial flow deficits and excesses, the wing loading parameter $n$ makes the Moore–Saffman family of profiles well suited to model the variety of axial flows present. Many factors influence the value of $n$ in practice because the wake- or jet-like characteristics of the vortex depend on the angle of attack and the tip configuration of the wing (see Chow et al. Reference Chow, Zilliac and Bradshaw1997; Anderson & Lawton Reference Anderson and Lawton2003). A high angle of attack or a round wingtip generally leads to an axial velocity excess at the vortex core (jet-like flow) which corresponds to a low value of  $n$ .

Figure 1.

Numerical solution giving (a) the tangential velocity profile $V_{n}(r)$ and (b) the axial velocity profile $W_{n}(r)$ obtained from the similarity solution of Moore & Saffman (Reference Moore and Saffman1973).

In the numerical investigation that follows, we will non-dimensionalize the velocity profiles given by the Moore–Saffman similarity solution in a way that facilitates comparison with results for the Batchelor vortex. To that end, we begin by defining the core circulation $\unicode[STIX]{x1D6E4}$ to be the circulation at the radius $r_{m}$ (or $\unicode[STIX]{x1D702}_{m}$ in the similarity variable) where the azimuthal velocity reaches its peak value. (Coincidentally, Spalart (Reference Spalart1998), in discussing the decay of trailing vortices, remarks that the core circulation ‘is as legitimate as any other measure’.) It follows from (2.1) that

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=\frac{\unicode[STIX]{x1D6FD}}{(\unicode[STIX]{x1D708}t)^{n/2}}r_{s}\unicode[STIX]{x1D6E4}_{n},\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{n}=2\unicode[STIX]{x03C0}\sqrt{-\unicode[STIX]{x1D702}_{m}}V_{n}(\unicode[STIX]{x1D702}_{m})$ is a dimensionless quantity wholly determined by the maximization of the function $V_{n}$ . We remark that $\unicode[STIX]{x1D6E4}_{n}$ depends on the value of $n$ , as can be inferred from figure 1(a), and must be obtained numerically since no closed-form expression is available. We will normalize all vortices by scaling them such that their core circulation is equal to the value for the Batchelor vortex, namely $\unicode[STIX]{x1D6E4}_{b}\simeq 4.4949$ .

We now define the length scale $r_{s}$ already introduced above and a characteristic velocity that will be used to scale both velocities, i.e.

(2.3a,b ) $$\begin{eqnarray}r_{s}=\sqrt{4\unicode[STIX]{x1D708}t}\quad \text{and}\quad V_{s}=\frac{\unicode[STIX]{x1D6E4}_{n}}{\unicode[STIX]{x1D6E4}_{b}}\frac{\unicode[STIX]{x1D6FD}}{(\unicode[STIX]{x1D708}t)^{n/2}}.\end{eqnarray}$$

The Reynolds number defined by these scales is $Re=r_{s}V_{s}/\unicode[STIX]{x1D708}$ . It depends on the wing loading parameter $n$ in two ways. First, the constants $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6E4}_{n}$ in the definition of $V_{s}$ vary slowly with $n$ . Second, the disparate nature of viscous diffusion on the length and velocity scales in (2.3), to wit $r_{s}\sim \sqrt{\unicode[STIX]{x1D708}t}$ and $V_{s}\sim (\unicode[STIX]{x1D708}t)^{-n/2}$ , causes the Reynolds number to be a different function of $\unicode[STIX]{x1D708}t$ for every value of $n$ . Regarding the latter comment, we will freeze the vortex flow in time (see § 3.1), which means fixing a time $t$ in (2.3).

In order to ensure that Batchelor’s model is approximated in the limit $n\rightarrow 1$ , we will insert the quantity $\bar{W}$ , which is the difference between the maximum and minimum of $W_{n}$ , into the expression for the axial velocity component only. The corresponding quantity for the Batchelor vortex is $\bar{W}_{b}=1$ . Taking $r_{s}$ and $V_{s}$ as in (2.3), the non-dimensionalized profiles that result are

(2.4a,b ) $$\begin{eqnarray}V(r)=\frac{\unicode[STIX]{x1D6E4}_{b}}{\unicode[STIX]{x1D6E4}_{n}}V_{n}(\unicode[STIX]{x1D702})\quad \text{and}\quad W(r)=W_{0}\frac{\bar{W}_{b}}{\bar{W}}W_{n}(\unicode[STIX]{x1D702}),\end{eqnarray}$$

with $\unicode[STIX]{x1D702}=-r^{2}$ and a dimensionless axial flow parameter $W_{0}$ given by

(2.5) $$\begin{eqnarray}W_{0}=\frac{\unicode[STIX]{x1D6E4}_{b}}{\unicode[STIX]{x1D6E4}_{n}}\frac{\unicode[STIX]{x1D6FD}}{W_{\infty }(\unicode[STIX]{x1D708}t)^{n/2}}\frac{\bar{W}}{\bar{W}_{b}}.\end{eqnarray}$$

The introduction of an axial flow parameter $W_{0}$ contrasts with the more common convention in the literature on vortex stability (Lessen, Singh & Paillet Reference Lessen, Singh and Paillet1974; Mayer & Powell Reference Mayer and Powell1992), where the magnitude of the axial velocity is normalized and a swirl parameter $q$ is introduced. Both are related in a straightforward way by $W_{0}=1/q$ . There are two reasons why the use of $W_{0}$ is preferable if one wishes to study the elliptical instability of the strained vortex. First, since the axial component $W(r)$ can be viewed as having a destabilizing effect on a stable vortex $V(r)$ in an unbounded domain, one expects the flow to be stable for $W_{0}<W_{0,c}$ , where $W_{0,c}$ is the stability threshold of the axial flow parameter. By using an axial flow parameter, calculations can be performed for a bounded interval of values of $W_{0}$ to get a complete picture of the elliptical instability in a linearly stable regime. We thus avoid having to deal with the unbounded interval $q>q_{c}$ . It should be noted that in the case of the Batchelor vortex, it was established by Lessen et al. (Reference Lessen, Singh and Paillet1974) that $q_{c}\simeq 1.5$ . This relates to a value of $W_{0,c}\simeq 0.66$ .

A second reason for using $W_{0}$ will become clear in § 3.1; the non-axisymmetric correction term induced by the strain field is, at first order, invariant to a rescaling of the azimuthal velocity. A change in $q$ therefore does not affect the strain field, making an axial flow parameter more appropriate for studying the stability of a strained vortex.

The Moore–Saffman family of velocity profiles describes vortices that develop from laminar flow over a wing. However, in free flight, where the Reynolds number may exceed $10^{6}$ , the airfoil boundary layer is likely to be turbulent. This motivated Phillips (Reference Phillips1981) to consider trailing vortices undergoing roll-up behind a lightly loaded wing on which the boundary layer is turbulent. Phillips determined the form of the Reynolds-stress distribution after identifying three regions in the vortex core. He then suggested that two modes for the decay of the peak azimuthal velocity in time are possible: ${\sim}t^{-n/2}$ , which is the rate for the profiles defined in this section, and ${\sim}t^{(n-2)/2}$ , a rate to which the decay may asymptote after some time (further behind the wingtip). Here, $t$ denotes time and $n$ is the wing loading parameter.

2.2 Inviscid outer region

The large- $r$ asymptotic behaviour of $V(r)$ and $W(r)$ in (2.4) is

(2.6a,b ) $$\begin{eqnarray}V(r)\sim r^{-n}\quad \text{and}\quad W(r)\sim r^{-2n},\end{eqnarray}$$

so that the circulation at infinity is unbounded for $n<1$ . This contrasts with other vortex models, such as the Batchelor vortex, which have a finite circulation. In order to make a comparison with those models more natural, the velocity profiles in (2.4) are matched to an inviscid outer solution by writing

(2.7a,b ) $$\begin{eqnarray}\tilde{V}(r)=\frac{V(r)}{\left(1+\left({\displaystyle \frac{r}{a}}\right)^{4}\right)^{(1-n)/4}}\quad \text{and}\quad \tilde{W}(r)=\frac{W(r)}{\left(1+\left({\displaystyle \frac{r}{a}}\right)^{4}\right)^{(1-n)/2}},\end{eqnarray}$$

where $a$ is a non-dimensional parameter that controls the viscous core radius. The velocity profiles in (2.7) involve a smooth transition from the Moore–Saffman profiles, in the vortex core near $r=0$ , to

(2.8a,b ) $$\begin{eqnarray}\tilde{V}(r)\sim \frac{1}{r}\quad \text{and}\quad \tilde{W}(r)\sim \frac{1}{r^{2}},\quad \text{for}\;r\rightarrow \infty .\end{eqnarray}$$

The form of the first formula in (2.7) closely resembles that of the second vortex model considered by Fabre & Jacquin (Reference Fabre and Jacquin2004a ). In all that follows, we will omit writing the tilde in (2.7) and (2.8).

The foregoing matching procedure impacts the far field and has little bearing on the viscous solution in the vortex core. For instance, for $n=0.5$ , the denominator in (2.7) affects the azimuthal velocity of the Moore–Saffman solution at $r=a/2$ by less than 1 % and its axial velocity at the same radius by approximately 1.5 %. The difference between the Moore–Saffman profiles from § 2.1 and their matched counterparts in (2.7) vanishes for $a\rightarrow \infty$ . Moreover, the concentration of the velocity profiles near the vortex centre on display in figure 1, combined with the comments at the start of this paragraph, suggests that the dependence of numerical results on $a$ may already be weak for $a\lesssim 6$ . Even so, the parameter $a$ must be chosen sufficiently large to prevent the outer region from encroaching on the viscous core. For more details about the effect of $a$ on the results of our calculations, we refer the reader to § 3.3.

2.3 Linear stability of a single vortex

Varying the wing loading parameter $n$ in the Moore–Saffman vortex model significantly impacts the linear stability characteristics (Feys & Maslowe Reference Feys and Maslowe2014). The key stability properties are reviewed here. They differ from those in the aforementioned investigation because the present paper includes the core circulation in the velocity scale (see § 2.1). Even so, the perturbation form of the governing equations considered in this section is the same as in Feys & Maslowe (Reference Feys and Maslowe2014). The mean profiles as $n\rightarrow 1$ correspond to the Batchelor vortex mean flow profiles up to an accuracy of $10^{-3}$ . The growth rates and regions of instability for the Batchelor vortex have been documented extensively (Mayer & Powell Reference Mayer and Powell1992; Ash & Khorrami Reference Ash, Khorrami and Green1995), and our findings are in very good agreement with those results. It should be noted that we have exchanged the swirl parameter $q$ for the axial flow parameter $W_{0}=1/q$ , as was stated in § 2.1, and that this causes the growth rates to scale by a factor of $W_{0}$ . The effect of viscosity on the inviscid modes of instability is illustrated. Many linear modes of instability of the Moore–Saffman vortex possess a critical layer singularity. The existence and role of the critical layer are briefly discussed at the end of this section.

Figure 2.

For azimuthal wavenumber $m=-1$ , both (a) the growth rate $\unicode[STIX]{x1D714}_{i}$ and (b) the corresponding axial wavenumber $k$ for the most unstable inviscid modes for all values in the $(W_{0},n)$ -plane are displayed. The contour spacing is 0.05 in (a) and 0.03 in (b), with the outermost contour representing the neutral curve. The horizontal dotted lines mark $n=0.44$ and $n=0.15$ .

With regard to the inviscid instability of the Moore–Saffman family of velocity profiles, figure 2 displays the growth rates and axial wavenumbers of the most amplified $m=-1$ disturbances at every point in the $(W_{0},n)$ -plane. First, we discuss figure 2(a). The contours on this figure are curves of constant growth rate, and the outermost contour represents the neutral curve (corresponding to $W_{0,c}$ ). The fastest growing mode was calculated independently for every value of the axial flow parameter and the wing loading parameter; each point on the figure represents the growth rate $\unicode[STIX]{x1D714}_{i}$ of the most unstable perturbation for fixed values of $W_{0}$ and $n$ . For example, the most unstable disturbance for $W_{0}=1$ and $n=0.50$ has a growth rate of $\unicode[STIX]{x1D714}_{i}\simeq 0.2595$ . Several conclusions are true for all values of the wing loading parameter $n$ . The axial component $W(r)$ has a destabilizing effect. For $W_{0}=0$ , the purely swirling flow is stable for all values of $n$ in the sense that all disturbances decay in time. The flow destabilizes when the axial flow parameter $W_{0}$ is increased past a critical value $W_{0,c}$ . As mentioned in § 2.1, $W_{0,c}\simeq 0.66$ for $n=0.95$ , a value of the wing loading parameter for which the Moore–Saffman profiles closely correspond to those of the Batchelor vortex. For $n=0.50$ , the threshold has decreased to approximately $W_{0,c}\simeq 0.24$ . Increasing $W_{0}$ past the critical value results in a further destabilization of the linear inviscid modes and an increase of the growth rate of the most amplified disturbance. Two important features of the Moore–Saffman vortex model are visible. First, a decrease in $n$ is associated with a decrease in the critical value $W_{0,c}$ . This is true until $n=0.44$ , a constant that is indicated by the upper horizontal dotted line. The special significance of this value of the wing loading parameter can be understood from figure 1(b): the change in sign of the core axial velocity at $n=0.44$ corresponds to a transition from wake-like to jet-like axial flow.

Returning to figure 2(a), the vortex is found to become more stable when $n$ is decreased below $n=0.44$ . No unstable disturbances are found for values of the Moore–Saffman parameter lower than $n\simeq 0.15$ , depicted on the figure by the bottom dotted line. The neutral curve can be seen to approach the line $n=0.15$ asymptotically. The fact that $W_{0,c}\rightarrow \infty$ means that no amount of axial flow will cause the damped modes to become unstable. Second, the growth rates are found to be greater for $0.44\leqslant n<0.95$ than for the Batchelor vortex, a fact that becomes more pronounced for large $W_{0}$ and small $n$ .

The effect of a change in $n$ on the axial wavenumber $k$ is clarified in figure 2(b). This figure plots the axial wavenumber of the most amplified mode at every point in the $(W_{0},n)$ -plane for which the growth rate was shown in figure 2(a). The contours are curves of constant $k$ . Only unstable modes have been computed; the calculation was halted at the neutral curve, which is visible as the leftmost curve. The most unstable disturbance for $W_{0}=1$ and $n=0.50$ , with the aforementioned growth rate of $\unicode[STIX]{x1D714}_{i}\simeq 0.2595$ , has $k\simeq 0.6852$ . Interesting variations in $k$ can be observed. Besides an increase in $k$ when $n$ is decreased, a diagonal ridge is visible which marks a region of rapid change in the axial wavelength of the most unstable disturbance. These topographical structures demonstrate that the introduction of an extra parameter $n$ by Moore and Saffman generalizes the known results for the Batchelor vortex in a non-trivial way. The normal modes presented in figure 2 are all primary modes of instability; other growing disturbances, with a lower growth rate, can exist for the same set of parameter values $m$ , $n$ and $k$ . This is true for all unstable azimuthal wavenumbers, and this is worth keeping in mind because triadic resonances between these different modes and the strain field have been observed in numerical simulations (Lacaze et al. Reference Lacaze, Ryan and Le Dizès2007).

Unstable perturbations with $m=1$ exist for all values of $n$ , but their growth rate is small compared with disturbances with $m<0$ . The Moore–Saffman vortex is stable for $m>1$ and $m=0$ . Even so, several of these stable azimuthal wavenumbers may participate in the elliptical instability process. It is, for example, possible for cooperative instability to cause growth in disturbances with azimuthal wavenumber $m=0$ through the resonance $(m_{1},m_{2})=(-2,0)$ . In fact, we will show in § 4 that this particular resonant pair manifests itself for the vortex model that is the subject of this paper.

As a further matter, viscosity will be present in what follows. The effect of a finite Reynolds number on the inviscid modes of instability, into which both Mayer & Powell (Reference Mayer and Powell1992) and Ash & Khorrami (Reference Ash, Khorrami and Green1995) offer insight, is of a strictly damping nature. Although the amplification rate of unstable modes is modified markedly at Reynolds numbers of the order of 100 and below, the stability properties of the flow for the range of Reynolds numbers considered in this paper ( $Re$ greater than $10^{3}$ ) bear a strong resemblance, both qualitatively and quantitatively, to those illustrated in figure 2 for the inviscid case. One minor difference concerns $W_{0,c}$ , the value of $W_{0}$ associated with the onset of instability, which equals approximately 0.203 for $Re=3180$ , a very slight increase from the inviscid case. For $W_{0}=0.482$ , instability occurs for $0.34\lesssim n\lesssim 0.81$ , with a maximum growth rate of $\unicode[STIX]{x1D714}_{i}\simeq 0.1131$ .

To conclude, we highlight the role of critical layer singularities in the linear modes of instability. The importance for the cooperative instability mechanism lies in the fact that the critical layer singularity damps neutral modes, thereby reducing the number of possible resonant pairs exhibiting growth (Sipp & Jacquin Reference Sipp and Jacquin2003). Nevertheless, resonant coupling between a neutral mode and a damped critical layer mode is possible if the coupling is sufficiently strong (Lacaze et al. Reference Lacaze, Ryan and Le Dizès2007). The location of the critical layer is determined by a positive real root $r_{c}$ of the expression

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}(r)=m\frac{V(r)}{r}+kW(r)-\unicode[STIX]{x1D714},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D714}_{i}$ is the complex frequency. In the absence of axial flow ( $W(r)\equiv 0$ ), it is clear that a singularity can arise only if the angular frequency $\unicode[STIX]{x1D714}/m$ of the normal mode is in the range of the angular momentum of the vortex. Moreover, changing the sign of $m$ only causes the sign of the real frequency $\unicode[STIX]{x1D714}_{r}$ to invert. The situation is complicated by a non-zero axial flow term which breaks this symmetry. The eigenvalue problem for $\unicode[STIX]{x1D714}$ has to be solved in that case and care has to be taken near the critical radius $r_{c}$ . To solve for the eigenmode, the contour of integration must be deformed away from the real positive axis over or under the singularity depending on whether the sign of $\unicode[STIX]{x1D6FE}^{\prime }(r_{c})$ is negative or positive. The viscous critical layer analysis of Le Dizès (Reference Le Dizès2004) shows how to obtain the correct eigenvalue condition for a neutral mode.

3.1 Governing equations

The governing equations are obtained in two steps. To begin with, a first-order expression for the elliptical deformation of the base mean flow is derived. Then, the perturbation equations are obtained by linearizing about the elliptical mean flow.

The non-axisymmetric correction to a base mean flow $[0,V_{0}(r),W_{0}(r)]$ , for the Moore–Saffman vortex in this study given by (2.7), subject to a stationary dipolar strain field of magnitude $\unicode[STIX]{x1D700}$ is detailed in Moore & Saffman (Reference Moore and Saffman1975) in the absence of axial flow ( $W_{0}(r)\equiv 0$ ). As remarked by Lacaze et al. (Reference Lacaze, Ryan and Le Dizès2007), the presence of an axial velocity component does not complicate matters since it is decoupled from the other velocity components at order $\unicode[STIX]{x1D700}$ . The strain field causes an elliptical deformation of the base mean flow, resulting in a velocity field $\bar{\boldsymbol{u}}$ equal to

(3.1) $$\begin{eqnarray}[\bar{u}(r,\unicode[STIX]{x1D703}),\bar{v}(r,\unicode[STIX]{x1D703}),\bar{w}(r,\unicode[STIX]{x1D703})]=[0,V_{0}(r),W_{0}(r)]+\unicode[STIX]{x1D700}[U_{1}(r,\unicode[STIX]{x1D703}),V_{1}(r,\unicode[STIX]{x1D703}),W_{1}(r,\unicode[STIX]{x1D703})].\end{eqnarray}$$

An important consequence of the external strain field is the appearance of a radial mean velocity component $U_{1}(r,\unicode[STIX]{x1D703})$ at order $\unicode[STIX]{x1D700}$ . The correction terms in (3.1) are given by

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle U_{1}(r,\unicode[STIX]{x1D703})=-\frac{F(r)}{r}\sin (2\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle V_{1}(r,\unicode[STIX]{x1D703})=-\frac{F^{\prime }(r)}{2}\cos (2\unicode[STIX]{x1D703}) & \displaystyle\end{eqnarray}$$

and

(3.4) $$\begin{eqnarray}W_{1}(r,\unicode[STIX]{x1D703})=-\frac{F(r)W_{0}^{\prime }(r)}{2V_{0}(r)}\cos (2\unicode[STIX]{x1D703}),\end{eqnarray}$$

where $F(r)$ satisfies the second-order ODE

(3.5) $$\begin{eqnarray}F^{\prime \prime }(r)+\frac{1}{r}F^{\prime }(r)-\left(\frac{4}{r^{2}}+\frac{\unicode[STIX]{x1D714}_{0}^{\prime }(r)}{V_{0}(r)}\right)F(r)=0.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6FA}_{0}=V_{0}(r)/r$ is the angular velocity of the base mean flow and $\unicode[STIX]{x1D714}_{0}=(rV_{0}(r))^{\prime }/r$ is the axial component of its vorticity. (It should be noted that (3.5) corrects the sign of the first term in the brackets in (2.15) of Moore & Saffman (Reference Moore and Saffman1975). It is not clear in their paper how combining (2.9) with (2.12) implies (2.13) having all positive signs. As far as we know, all papers since Moore & Saffman (Reference Moore and Saffman1975) use the correct expression given in Saffman (Reference Saffman1992).) The boundary conditions for $F(r)$ are prescribed by the behaviour of $V_{0}(r)$ for $r\rightarrow 0$ and $r\rightarrow \infty$ . With regards to $r\rightarrow 0$ , the rigid rotation exhibited by the vortex at its core leads to $F(r)=O(r^{2})$ or, equivalently, $F(r)\sim S_{0}r^{2}$ for some constant $S_{0}$ . In the far field, the behaviour of the azimuthal velocity profile in this study is identical to that of the Batchelor vortex due to the inclusion of an inviscid outer flow (see § 2.2). The ensuing requirement is that $F(r)\rightarrow r^{2}$ as $r\rightarrow \infty$ . This corresponds to the assumption of a unitary external strain rate and hence determines the normalization of $F(r)$ . The two aforementioned conditions allow us to determine $F(r)$ in a unique way by solving (3.5) numerically.

Figure 3.

For a unitary external strain field, the strain rate $S_{0}$ at the vortex centre is shown as a function of the wing loading parameter $n$ .

In practice, the fact that $F(r)$ becomes unbounded for $r\rightarrow \infty$ poses a computational obstacle. It is resolved by defining $G(r)=F(r)/r^{2}$ and solving the ODE for $G(r)$ instead. One finds $G(r)\rightarrow 1$ and $G^{\prime }(r)\rightarrow 0$ as $r\rightarrow \infty$ , while a boundary condition at $r=0$ can be obtained easily from that for $F(r)$ . The constant $S_{0}$ is then determined using a fourth-order Runge–Kutta shooting method by integrating the equation governing $G(r)$ inward from a large radius and outward from $r=0$ , and matching both at an intermediate radius. As mentioned in the previous paragraph, the solution for $F(r)$ is unique. It follows that the constant $S_{0}$ , which measures the strain rate at the vortex centre, is unique as well. For the Moore–Saffman family of profiles introduced in § 2, the magnitude of $S_{0}$ depends on the value of the wing loading parameter $n$ . Figure 3 shows the precise nature of this dependence. An increase of the strain rate at the vortex centre can be observed for decreasing $n$ , from a value of $S_{0}\simeq 2.54$ for $n=0.99$ to $S_{0}\simeq 3.15$ for $n=0.20$ . Physically, this increase signifies that the unitary external strain field is felt more strongly at the vortex centre for $n<1$ . We will refer back to this observation at the point of interpreting our numerical results in § 4.

In order to obtain equations for a velocity–pressure perturbation $(\tilde{\boldsymbol{u}},\tilde{p})$ of order $\unicode[STIX]{x1D705}$ , the three-dimensional incompressible Navier–Stokes equations

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

with $Re$ the Reynolds number as defined in § 2.1, are linearized about the non-axisymmetric mean flow given by (3.1). The linearized governing equations are then expanded in cylindrical coordinates in a manner similar to what was done in Ash & Khorrami (Reference Ash, Khorrami and Green1995). The fluctuations to the steady state are assumed to be of the form

(3.8) $$\begin{eqnarray}\displaystyle & \boldsymbol{u}(r,\unicode[STIX]{x1D703},z,t)=\bar{\boldsymbol{u}}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D705}\tilde{\boldsymbol{u}}(r,\unicode[STIX]{x1D703},z,t), & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & p(r,\unicode[STIX]{x1D703},z,t)=\bar{p}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D705}\tilde{p}(r,\unicode[STIX]{x1D703},z,t). & \displaystyle\end{eqnarray}$$

They lead to the following equations governing the perturbation:

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\tilde{\boldsymbol{u}}}{\unicode[STIX]{x2202}t}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\bar{\boldsymbol{u}}+(\bar{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\tilde{\boldsymbol{u}}=-\unicode[STIX]{x1D735}\tilde{p}+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\tilde{\boldsymbol{u}}, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0. & \displaystyle\end{eqnarray}$$

The numerical method used to solve the linearized viscous equations given by (3.10) and (3.11) is described in the following section. The domain is unbounded and the only boundary condition that is applied is the decay of the perturbation quantities at infinity. As is usually presumed (Fabre, Sipp & Jacquin Reference Fabre, Sipp and Jacquin2006; Lacaze et al. Reference Lacaze, Ryan and Le Dizès2007), the slow viscous diffusion of the mean flow $\bar{\boldsymbol{u}}$ has been neglected in (3.8) and (3.9). The background flow is effectively frozen in time, so that any temporal growth can be directly attributed to the perturbation  $\tilde{\boldsymbol{u}}$ .

3.2 Numerical method

The three-dimensional linearized Navier–Stokes equations, (3.10) and (3.11), are to be solved in cylindrical coordinates using a spectral method in space and a third-order splitting method in time. The problem of spatially discretizing the perturbation quantities $(\tilde{\boldsymbol{u}},\tilde{p})$ requires some thought because the physical domain is unbounded and includes the pole at $r=0$ . Clustering of collocation points near the coordinate singularity at the origin can restrict the size of the explicit time step. Care is taken to avoid this by choosing a suitable orthogonal basis of polynomials in the radial direction. First, the details of the spectral discretization in space are given and compared with other standard numerical approaches. Then, the time-stepping mechanism is described. Finally, the parameters of the numerical simulation are given.

Spectral methods are a popular choice due to their high accuracy and relative ease of use. The periodicity of the azimuthal variable $\unicode[STIX]{x1D703}$ allows a Fourier expansion. The axial variable $z$ does as well, provided that the axial dimension, which is inherently unbounded, is assumed to be periodic. In the radial direction, however, there are numerous families of polynomials that form suitable orthogonal bases. Our numerical routine employs mapped associated Legendre polynomials in the radial variable $r$ . The specifics of this spectral expansion were first formulated by Matsushima & Marcus (Reference Matsushima and Marcus1997). The orthogonal polynomials satisfy the pole condition exactly at the coordinate singularity $r=0$ , and their behaviour as $r\rightarrow \infty$ is suitable for expanding smooth functions that decay algebraically or exponentially. The associated Legendre polynomials, denoted by $P_{n}^{m}$ , generalize the more familiar Legendre polynomials and satisfy a similar recursion. Formulae detailing the recursion relations necessary for manipulating the polynomials $P_{n}^{m}$ in spectral function space can be found in Abramowitz & Stegun (Reference Abramowitz and Stegun1972).

Following Matsushima & Marcus (Reference Matsushima and Marcus1997), the radial variable $0\leqslant r<\infty$ is mapped to a computational variable $-1\leqslant \unicode[STIX]{x1D707}<1$ via the algebraic map

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{r^{2}-L^{2}}{r^{2}+L^{2}},\end{eqnarray}$$

where $L>0$ is a map parameter that can be chosen to optimize convergence. After a considerable amount of algebra, (3.12) eliminates $r$ from (3.10) and (3.11) by substituting in the variable $\unicode[STIX]{x1D707}$ instead. All perturbation quantities are then expanded into Fourier/associated Legendre series as follows:

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{u}}(\unicode[STIX]{x1D707},\unicode[STIX]{x1D703},z,t)=\mathop{\sum }_{k=-K}^{K}\mathop{\sum }_{m=-N}^{N}\mathop{\sum }_{n=|m|}^{N}\tilde{\boldsymbol{u}}_{k,m,n}(t)P_{n}^{m}(\unicode[STIX]{x1D707})\text{e}^{\text{i}m\unicode[STIX]{x1D703}}\text{e}^{\text{i}kz}, & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{p}(\unicode[STIX]{x1D707},\unicode[STIX]{x1D703},z,t)=\mathop{\sum }_{k=-K}^{K}\mathop{\sum }_{m=-N}^{N}\mathop{\sum }_{n=|m|}^{N}\tilde{p}_{k,m,n}(t)P_{n}^{m}(\unicode[STIX]{x1D707})\text{e}^{\text{i}m\unicode[STIX]{x1D703}}\text{e}^{\text{i}kz}. & \displaystyle\end{eqnarray}$$

Here, the integers $N$ and $K$ determine the number of basis functions in the expansion. The complex coefficients $\tilde{\boldsymbol{u}}_{k,m,n}(t)$ and $\tilde{p}_{k,m,n}(t)$ are the expansion coefficients for the velocity and pressure respectively. The truncation of the radial basis in (3.13) and (3.14) is triangular in the sense that higher azimuthal wavenumbers $m$ are associated with fewer radial basis functions. This happens naturally as the associated Legendre polynomials satisfy $P_{n}^{m}(\unicode[STIX]{x1D707})=0$ for $n<|m|$ . Boyd (Reference Boyd2000) mentions that for spherical harmonics, a triangular truncation gives the most uniform resolution (and allows the longest time step). Although this has not been proven for the associated Legendre polynomial basis, it seems plausible that it is preferable here as well. Quadrature in the variable $\unicode[STIX]{x1D707}$ is performed using Gauss–Legendre collocation points. When mapped to the variable $r$ , half of the collocation points are guaranteed to lie between 0 and $L$ . In this way $L$ controls the resolution near the core of the vortex.

The application of associated Legendre polynomials to pseudospectral methods is relatively novel, but was shown to be successful for simulating cooperative instability in Matsushima & Marcus (Reference Matsushima and Marcus1997). More recently, Bristol et al. (Reference Bristol, Ortega, Marcus and Savaş2004) have investigated the instabilities of a pair of parallel unequal-strength vortices using the same numerical method. The spectral expansion used by these authors differs from the one used in this paper in an important way. Both Matsushima & Marcus (Reference Matsushima and Marcus1997) and Bristol et al. (Reference Bristol, Ortega, Marcus and Savaş2004) performed a toroidal and poloidal decomposition of the velocity field which eliminates the incompressibility condition (3.11) and reduces the number of unknowns. In lieu of this, we have taken adjoint complex combinations of the radial and azimuthal velocity profiles, thus simplifying the vector Laplacian. Similar considerations apply to the spectral tau-method proposed by Sakai & Redekopp (Reference Sakai and Redekopp2009). That method employs one-sided Jacobi polynomials as radial expansion functions for solving vector field equations defined in polar coordinates.

The spectral method described above has significant advantages over other numerical techniques. When compared with an approach in Cartesian coordinates using a Fourier expansion in all directions, such as the one taken by Faddy & Pullin (Reference Faddy and Pullin2005) in conjunction with a special mapping for the axial coordinate developed by Cain, Ferziger & Reynolds (Reference Cain, Ferziger and Reynolds1984), our expansion in cylindrical coordinates requires fewer basis functions to achieve a desired accuracy and allows a direct inspection of the azimuthal wavenumbers. The disadvantage is that there is no fast transform to the space of associated Legendre polynomials, thus increasing computational expense. In cylindrical coordinates, there are other polynomial bases that can be used to expand in the radial direction. Several of these are listed in Boyd (Reference Boyd2000), including Bessel functions and Chebyshev polynomials. Bessel functions have been shown to be ill-suited because of slow convergence of the expansion. Gottlieb & Orszag (Reference Gottlieb and Orszag1977) suggest using Chebyshev polynomials for a cylindrical geometry because a fast transform is available. Clustering of the grid points occurs near $r=0$ , however, which severely limits the time step.

Stepping in time in our numerical procedure is performed using a third-order stiffly stable splitting method devised by Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991). This scheme is similar to Adams-type explicit schemes for time integration, but it has a significantly broader stability region due to the fact that the time derivative is discretized using backwards differentiation. Moreover, the convective terms in (3.10) are handled explicitly while the diffusive terms are implemented implicitly for stability reasons. The lack of a predictive equation for the evolution in time of the pressure stems from the fact that, in an optimization sense, the pressure $\tilde{p}$ is a Lagrange multiplier for the momentum equations associated with the divergence-free constraint. As a consequence, at every time step (of size $\unicode[STIX]{x0394}t$ ), a Poisson equation of the form

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\tilde{p}=\frac{1}{\unicode[STIX]{x0394}t}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}\end{eqnarray}$$

must be solved to determine the pressure necessary to make the velocity field $\tilde{\boldsymbol{u}}$ satisfy (3.11). After being initialized to random noise, the disturbance velocity field $\tilde{\boldsymbol{u}}$ is made divergence-free by applying a pressure gradient that is calculated in this manner. The disturbance velocity and pressure fields are then advanced in time and quickly become a combination of exponentially growing and decaying modes. In order to compute the azimuthal wavenumbers $m$ associated with maximum amplification for a given value of $k$ , we track the evolution in time of the kinetic energy $E(m)$ of each azimuthal Fourier mode in (3.13). The gain $\unicode[STIX]{x1D714}_{i}(m)$ in its energy is obtained using the following simple relation (Faddy Reference Faddy2005):

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{i}(m)=\lim _{t\rightarrow \infty }\frac{1}{2}\frac{\text{d}\ln E(m)}{\text{d}t}.\end{eqnarray}$$

Convergence of this quantity is found to be rapid in most cases. A comparison of the magnitude of $\unicode[STIX]{x1D714}_{i}(m)$ for all modes in the truncation admits identification of the two azimuthal wavenumbers $m_{1}$ and $m_{2}$ with the largest growth rate $\unicode[STIX]{x1D714}_{i}(m_{1})\simeq \unicode[STIX]{x1D714}_{i}(m_{2})$ .

The equations necessary for time-stepping the expansion coefficients $\tilde{\boldsymbol{u}}_{k,m,n}(t)$ and $\tilde{p}_{k,m,n}(t)$ given in (3.13) and (3.14), which are obtained from (3.10) and (3.15) via the splitting method formulation, have a convenient form in the sense that coefficients with different azimuthal wavenumbers $m$ are mostly decoupled. More specifically, both the pressure gradient and the viscous diffusion term are applied to every azimuthal wavenumber independently. There is, however, coupling between azimuthal wavenumbers (which are exactly 2 apart, i.e. $m-2$ , $m$ and $m+2$ ) through the effect of the non-axisymmetric correction to the mean flow in the convection term of the governing equation. Nonetheless, the equations allow for an extensive parallelization of the numerical code which in turn reduces simulation runtime drastically. The method is implemented in Fortran 95 using the OpenMP interface for parallel processing.

We now discuss the range of values of the computational parameters that were used in the numerical simulations. The resolution parameter $L$ was found to have a significant impact on the rate of convergence and it is thus necessary to choose this parameter carefully. The optimal value of $L$ is found to lie between 0.5 and 6 and depends primarily on the axial wavenumber $k$ . The radial decay of the modes increases for decreasing $k$ and the collocation points have to be spread out accordingly to ensure accurate results. The time step $\unicode[STIX]{x0394}t$ was chosen to be of order 0.001. This size of time step was found to be sufficient to integrate the governing equation (3.10) as far as $t=1000$ . The asymptotic behaviour of the perturbation velocity field was thereby recovered and the dominant modes could be obtained. Axial wavenumbers were usually chosen to be between 1 and 7 with steps of 0.01, unless otherwise stated. The integer $N$ in (3.13) and (3.14) was chosen to be of the order of 60, although accurate results could already be obtained for $N$ as low as $25$ . Numerically, $N$ determines the number of azimuthal (and radial) wavenumbers. For example, $N=25$ corresponds to 51 Fourier modes in the azimuthal direction. When $N$ is increased, which corresponds to increasing the number of azimuthal and radial modes, the time step $\unicode[STIX]{x0394}t$ is decreased. Specifically, $\unicode[STIX]{x0394}t$ is taken to be inversely proportional to  $N$ .

3.3 Numerical convergence and code verification

Several checks were performed in order to assure that the numerical routine functions as desired. In this section, we compare some of our results with previously published data and obtain excellent agreement. In addition, we present convergence analyses of the numerical method in order to show that in most cases a very low number of spectral modes is needed to get an accurate estimate of the instability growth rate. Results are first presented for instabilities in the absence of strain. After this, a strain field is introduced and elliptical instability growth rates are verified. Finally, the effect of the physical parameter $a$ , introduced in (2.7), on the numerical results is illustrated.

To begin with, we reproduce several published results for modes of instability of the Batchelor vortex in the absence of strain ( $\unicode[STIX]{x1D700}=0$ ). In table 1, comparisons are presented for two unstable modes: an inviscid mode in the presence of viscosity (case I) and a viscous centre mode (case II). For case I, the governing equations were integrated to $t=300$ . For case II, due to the smaller growth rate of the instability, we integrated until $t=800$ . As is well known, the growth rate of centre modes is at least an order of magnitude smaller than that of inviscid modes. In both cases, convergence is observed with an increase in the number of basis functions. Looking at case I, the growth rate of the instability is already determined quite accurately for $N=10$ . Little gain in accuracy is perceived in runs with progressively higher values of $N$ . We find the growth rate to be $\unicode[STIX]{x1D714}_{i}=0.32446063$ , all digits being significant. Excellent agreement with the numerical results of Fabre & Jacquin (Reference Fabre and Jacquin2004b ) and Mao (Reference Mao2010) is obtained. Case II is a viscous centre mode, which is a mode of instability with a perturbation field concentrated in the vortex centre and which occurs only in the presence of viscosity. Convergence is found to be somewhat slower than for case I, which is expected due to the nature of its perturbation velocity field and small growth rate. Viscous centre modes are not expected to participate in the elliptical instability mechanism; case II was examined only to verify our code. The numerical routine yields a growth rate of $\unicode[STIX]{x1D714}_{i}=0.019734$ , all digits being significant. Many recent computations have been made for this type of mode, for example by Fabre & Jacquin (Reference Fabre and Jacquin2004b ), Heaton (Reference Heaton2007b ) and Mao (Reference Mao2010). Table 1 demonstrates that our numerical result agrees well with all three of these investigations.

Table 1.

Code verification against published results for the Batchelor vortex for two cases of vortex instability. Case I is an unstable inviscid mode in the presence of viscosity for $Re=1000$ , $q=0.761$ , $m=-3$ and $k=1.659$ . Case II corresponds to an unstable viscous centre mode for $Re=14000$ , $q=2$ , $m=-1$ and $k=0.268$ . The change column presents the relative error (in per cent) between the growth rate in every row and the following one.

Next, we look at the Batchelor vortex subject to a strain field of magnitude $\unicode[STIX]{x1D700}=0.0625$ . Figure 4 shows numerical convergence results representative of those for all simulations of elliptical instability performed in this study. The results displayed in this figure were obtained for a specific set of flow parameter values and $n=0.95$ , a value of the wing loading parameter for which the mean flow profiles closely resemble those of the Batchelor vortex. We refer to this set-up as case III. For the flow parameters chosen, a resonance occurs between the strain field and the pair of azimuthal wavenumbers $(m_{1},m_{2})=(-2,0)$ . For each sample point on the curve, a random initial perturbation velocity field was advanced in time. As was mentioned in § 3.2, the perturbation growth rate is extracted by running the simulation for a sufficiently long time. To ensure accurate results, the equations for case III presented here were solved up to time $t=400$ . The most unstable perturbation component quickly develops exponential growth, which allows us to extract a numerical growth rate. The growth rate $\unicode[STIX]{x1D714}_{i}$ of the instability was found to be close to 0.0851.

Figure 4.

(a) The convergence of the growth rate $\unicode[STIX]{x1D714}_{i}$ with an increasing number of modes $N$ for case III. The following parameters are fixed: $n=0.95$ (approximately Batchelor), $\unicode[STIX]{x1D700}=0.0625$ , $W_{0}=0.482$ , $Re=3180$ and $k=1.88$ . The resulting elliptical instability is characterized by a resonance between the strain field and the azimuthal wavenumbers $(m_{1},m_{2})=(-2,0)$ . In (b), the region near the asymptotic value of $\unicode[STIX]{x1D714}_{i}$ has been magnified. In both plots, the dots represent sample points.

The present results compare favourably with the single previous investigation of this case of which we are aware. The numerical result given by Lacaze et al. (Reference Lacaze, Ryan and Le Dizès2007) is $\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D700}\simeq 1.23$ , which corresponds to $\unicode[STIX]{x1D714}_{i}=0.0769$ . The same authors find a theoretical prediction equal to $\unicode[STIX]{x1D714}_{i}/\unicode[STIX]{x1D700}\simeq 1.32$ , and this corresponds to $\unicode[STIX]{x1D714}_{i}\simeq 0.0825$ . With regard to convergence, figure 4(a) shows that our computation rapidly approaches an asymptotic value between $N=16$ and $N=22$ . For higher values of $N$ , a slight oscillation in $\unicode[STIX]{x1D714}_{i}$ is visible, but its amplitude is small and decreasing. A magnified version of the same curve is presented in figure 4(b), and it shows the convergence behaviour more clearly. The oscillation amplitude can be seen to diminish rapidly to less than $10^{-4}$ for $N\approx 50$ . This provides clear evidence of the correct functioning of our numerical scheme. A look at figure 5 confirms this conclusion. This figure displays the perturbation vorticity field associated with the mode of resonance in case III. The contours shown in this figure can be compared with those in figure 15(a) from Lacaze et al. (Reference Lacaze, Ryan and Le Dizès2007); the modal structures are seen to resemble one another closely.

Figure 5.

The perturbation vorticity field for the unstable mode named case III. The elliptical instability structure is visible as a combination of modes with azimuthal wavenumbers $m_{1}=-2$ and $m_{2}=0$ . The contours, which are isovorticity lines, are spaced 0.2 apart. The zero-vorticity contour is not displayed. Two poles of negative vorticity (with a minimum value of $-1$ ) and two poles of positive vorticity (with a maximum value of approximately 0.6) exist.

To conclude this section, we briefly revisit the discussion in § 2.2 about the dependence of our numerical results on $a$ , which is the non-dimensional parameter that dictates the point of transition from the Moore–Saffman profiles, in the inner region of the vortex, to the inviscid outer solution described by (2.8). Throughout our calculations we have observed a negligible dependence on $a$ when its value is sufficiently large. To illustrate this, figure 6 shows the growth rates $\unicode[STIX]{x1D714}_{i}$ as a function of $a$ for two trial cases (both modes of elliptical instability). The physical parameters for the two plots in this figure are quite different, yet we observe the same behaviour in both: the growth rates change somewhat for $a\lesssim 3$ and converge when $a$ is greater than that value. The fact that the results remain mostly invariant under a change of $a$ is a subtle but direct consequence of our choice in § 2.1 to normalize the core circulation instead of the far-field circulation. We used $a=5$ to generate the numerical results in § 4, unless stated otherwise.

Figure 6.

The dependence of the growth rate $\unicode[STIX]{x1D714}_{i}$ on the core radius parameter $a$ is displayed for elliptical modes of instability for (a) $\unicode[STIX]{x1D700}=0.0625$ , $W_{0}=0.482$ , $Re=3180$ , $k=1.88$ and $n=0.50$ and (b) $\unicode[STIX]{x1D700}=0.0625$ , $W_{0}=0.15$ , $Re=20000$ , $k=4.5$ and $n=0.70$ .