Linear instability of a perturbed Lamb–Oseen vortex

The propagation of waves on vortices, as well as their stability, was first considered by Kelvin and by Rayleigh in the latter part of the 19th century. Of particular importance in the present paper is Rayleigh's 1916 stability criterion for axisymmetric disturbances, namely, that a sufficient condition for stability is that the square of the circulation increases in the radially outward direction for all r. Von Karman in 1934 showed how this result could be recovered using considerations of centripetal force and pressure gradient on a displaced fluid element. Given that we will employ the normal mode approach, a third derivation of Rayleigh's criterion by Synge in 1933 is, for this study, the most pertinent. Synge considered the linear disturbance equation for axisymmetric modes and then used Sturm–Liouville theory to extend Rayleigh's criterion by showing that it is a sufficient condition for instability. A detailed presentation of this background material can be found in section 15 of the monograph by Drazin and Reid (2004).

Although the stability of vortices is relevant to a number of phenomena in engineering and geophysical fluid dynamics, much of the recent interest has been motivated by the aircraft trailing vortex problem. The introduction of jumbo jets such as the A380 poses a major hazard to following aircraft given that the strength of the vortices shed from the wingtip of an airplane is related to its weight. With that application in mind, Batchelor (1964) found a similarity solution valid far downstream of an aircraft. Lessen et al (1974) were the first of many to investigate the linear stability of the Batchelor vortex, which they approximated as the locally parallel flow

$$\begin{equation} V (r) = \frac{q}{r} \left ( 1- e^{- r^2}\right) \hspace{0.15in} \mbox{and}\hspace{0.15in} W (r) = e^{- r^2}, \end{equation} \tag{ 1 }$$

where q is a swirl parameter. It was shown by these authors that the axial flow is stabilized for $q \geqslant 1.5$, approximately, and this result was seemingly confirmed by a number of subsequent investigations, the most detailed being that of Mayer and Powell (1992). Heaton (2007), however, showed that weakly amplified center modes (i.e. modes having a critical point singularity near r = 0) exist and these modes are unstable until $q \geqslant 2.31$.

For the purposes of this paper, it will be convenient to follow Le Dizès and Lacaze (2005) by rescaling the variables so that q = 1 and the axial velocity in (1) is $\bar{W} = W_{0}\, e^{- r^2}$. When $\bar{W} = 0$, the vortex is obviously stable and the profile $\bar{V} (r)$ in (1) reduces to the model known as the Lamb–Oseen vortex. Despite the absence of linear instability in the case of a single Lamb–Oseen vortex, normal mode perturbations are still of interest in the context of aircraft trailing vortices. That is because the inertial modes known as Kelvin waves propagating on a pair of counter-rotating vortices have been used in studies of the cooperative, short-wave instability mechanism ; a review of early analyses of this topic can be found in the monograph by Saffman (1992). This mechanism, now usually termed 'elliptic instability,' is a consequence of the strain field induced on one vortex by a second vortex, coupled with the resonance of two Kelvin modes with the same axial wavenumber k and frequency ω and azimuthal wavenumbers m and m + 2.

Assuming that the external strain field was weak, Tsai and Widnall (1976) employed the discontinuous Rankine vortex model in their investigation of elliptic instability and found that the largest growth rates occurred for a pair of stationary Kelvin modes with $m = \pm 1$. (The analysis of Kelvin modes propagating on a Rankine vortex is outlined in Saffman (1992)). A natural extension of that work would be to apply the theory to more realistic vortices and the Lamb–Oseen model is an obvious candidate. Such extensions were carried out, for example, by Le Dizès and Laporte (2002) and Sipp and Jacquin (2003). Critical point singularities arise in the inviscid problem when the vortex profile is continuous and it was shown in the latter paper that the presence of critical points changes some of the conclusions reached by Tsai and Widnall. The resonance involving the bending modes $m = \pm 1$, however, is still effective provided that ω = 0, because there is no critical point in that case.

In order to identify possible resonances that give rise to elliptic instability, it is necessary to have available dispersion curves for the Kelvin waves that are supported by Lamb–Oseen vortices. The goal of describing all the inertial modes that might propagate on a continuous vortex is, of course, of interest in its own right and was the object of two very detailed studies that appeared only one year apart. The first, Le Dizès and Lacaze (2005), was essentially an inviscid analysis assuming the axial wavenumber k to be large and employing the WKB method. Critical points are generally present for non-axisymmetric modes and an indented contour of integration, consistent with a viscous critical layer, was employed to integrate around such points. The second article, by Fabre et al (2006), dealt with the viscous equations and integrations were carried out at large Reynolds numbers. It was found that most of the viscous solutions approached the already known, generally singular, inviscid modes as $\textrm{Re} \rightarrow \infty$. In addition to these Kelvin modes, Fabre et al discuss slightly damped disturbances termed 'critical layer waves' that do not have an inviscid limit. The possible relevance of the foregoing results to turbulence and vortex breakdown was pointed out in Fabre et al (2006) and references are given therein to such applications.

The objective of this paper is to study the instability of a single unbounded vortex, a case that has received remarkably little attention compared with the instability of Taylor vortices, which has been studied exhaustively since 1923. In the introduction of their paper, Fabre et al (2006) remark that in most of the recent literature on vortex instabilities, the instability is generated either by an axial flow or by the presence of more than one vortex. A comprehensive review of the single vortex case has been given by Ash and Khorrami (1995), including many details about the stability analyses of the Batchelor vortex. Here, on the other hand, we will focus on the possible instability of a single vortex with no axial flow. Specifically, we will investigate the linear stability characteristics of the family of vortices

$$\begin{equation} \bar{V} (r) = \frac{1 - (1 - \varepsilon r^2 ) e^{- r^2}}{r}, \end{equation} \tag{ 2 }$$

where the parameter ε will be varied and we know that the vortex is stable when ε = 0. It will be of particular interest to see if the computed growth rates are large enough to overcome the effect of viscosity which would likely be stabilizing. Figure 1 shows the velocity profile (2) for three values of ε. It should be remarked that these profiles retain the desired characteristic behavior for large and small values of r if the profile is to model trailing vortices. Specifically, $\bar{V}\sim r^{-1}$ in the far field, whereas $\bar{V}\sim (1 + \varepsilon) r$ near the vortex center, i.e. there is a slight increase in the angular velocity compared with the Lamb–Oseen vortex, but there is still rigid rotation locally. This property is apparent in figure 2(a) of Feys and Maslowe (2014) who determined numerically the velocity profiles yielded by the Moore-Saffman similarity solution for a trailing vortex and then investigated its stability for different values of the wing loading parameter.

Zoom In Zoom Out Reset image size

To show that the perturbed vortex (2) is at least locally unstable it is convenient to first note that the Rayleigh discriminant can be expressed in terms of the angular velocity $\bar{\Omega}$ and the axial vorticity ζ by writing $\Phi = 2\,\bar{\Omega}\,\zeta$. For the particular case of the velocity profile (2), we obtain

$$\begin{equation} \Phi = 2\,\bar{\Omega}\,\frac{(r\bar{V})^{^{\prime}}}{r} = 4\,\bar{\Omega}\, \left[ 1 + \varepsilon \left(1 - r^2\right)\right]\, e^{- r^2}. \end{equation} \tag{ 3 }$$

According to Sturm–Liouville theory, there will be instability, so long as Φ is negative for some value of r and we see from (3) that this condition can always be satisfied for the velocity profile (2). (For $\varepsilon \ll 1, \Phi$ will be negative only very far from the vortex center, so the instability will be extremely weak. Although instability is also possible for ε < 0, the corresponding growth rates were found to be small; we will therefore present results only for positive ε.)

Non-axisymmetric perturbations are also of interest in our study and for that case, a second theorem due to Rayleigh is pertinent. The latter theorem applies only if k = 0 and it provides a necessary condition for instability, namely, that $\zeta^{^{\prime}}$ change sign at some value of r. For the velocity profile (2), we obtain as the necessary condition

$$\begin{equation} \frac{d \zeta}{d r} = - 2 r \left[1 + \varepsilon \left(2 - r^2\right)\right] e^{- r^2} = 0, \end{equation} \tag{ 4 }$$

for some value of r. Equation (4), unlike (3), is not a sufficient condition for instability and we will see that although this condition is satisfied in our case, it does not lead to instability..

Michalke and Timme (1967) have reported an investigation of unstable vortices whose spirit is similar to that of the present paper. However, the application of primary interest to these researchers was to the stability of spanwise vortices that appear downstream of the initial instability in a two-dimensional mixing layer. They considered both broken-line profiles and the continuous vortex with velocity profile

$$\begin{equation} \bar{V} = \frac{[ 1 - \exp\,( - 2 r^6/3)]}{r}. \end{equation} \tag{ 5 }$$

The profile (5) was chosen to reproduce features observed in mixing layer experiments. In particular, the vorticity is concentrated in a band more or less centered at r = 1, which is where the maximum occurs. The Rayleigh discriminant $\Phi \geqslant 0$ for the above profile, so it is stable to axisymmetric perturbations. However, a sufficient condition for instability derived by Howard and Gupta (1962) suggests that instability is possible for helical perturbations with $k \ll m$. Michalke and Timme showed that, in fact, such instabilities do occur for their profile with the largest growth rates obtained for k = 0. Since $\zeta^{^{\prime}} = 16 r^3(1 - r^6) \exp( - 2 r^6/3)$, the necessary condition for instability (4) is satisfied for the velocity profile (5). It is, nonetheless, surprising that the largest growth rates were obtained for k = 0. We will see in section 3 that this is not the case for our velocity profile (2); further comparisons will be deferred until that section.

Although there are no general stability criteria for non-axisymmetric perturbations, an approximate generalization of Rayleigh's sufficient condition for instability was derived by Billant and Gallaire (2005). Their development assumes that k is large and employs the WKB method in a manner similar to Le Dizès and Lacaze (2005). We will see in section 3 that predictions using their analysis agree well with our numerical results for the onset of instability for the first two helical modes even though assumptions in their analysis about the eigenfunction are not satisfied.

Finally, there are important applications in geophysical fluid dynamics, where the stability of vortices is pertinent. Gent and McWilliams (1986) made the point that vortices with velocity profiles behaving as potential vortices in the far field decay too slowly to accurately describe, for example, tornadoes. For that reason, they investigated several models with velocity profiles such as $\bar{V} = r \exp(- r^{\alpha})$ that have zero total circulation in the far field. It will be seen in section 3 that although the far field velocity profile of a vortex does influence its stability, the most unstable modes of the perturbed Lamb–Oseen vortex do behave in a manner similar to those for the vortices considered by Gent and McWilliams.

Our objective, as stated above, is to investigate the instability of a vortex whose velocity components, in cylindrical coordinates, can be written $[0, \bar{V} (r), 0]$. To that end, a helical perturbation of $\textit{O}(\varepsilon)$ is superimposed on the mean flow and we write for the radial, azimuthal and axial velocity, respectively,

$$\begin{equation} [u_r,u_{\theta},u_z] = [0,\bar{V}(r),0]+\varepsilon [ u(r),v(r),w(r)] e^{i(kz+m\theta-\omega t)} \:. \end{equation} \tag{ 6 }$$

The pressure

$$\begin{equation} P = \bar{p}(r)+\epsilon p(r) e^{i(kz+m\theta -\omega t)} , \end{equation} \tag{ 7 }$$

where $\bar{p}(r)$ is the pressure of the unperturbed vortex. The axial wavenumber is denoted k, m is the azimuthal wavenumber and $\omega = \omega_r+i \omega_i$ is the complex frequency. Substituting these expressions into the inviscid momentum equations, the continuity equation, and then linearizing, one can eliminate all variables but u to obtain the stability equation derived by Howard and Gupta (1962), namely,

$$\begin{equation} \gamma^2 D\{S D_{*}u\} - \displaystyle{\left\{\gamma^2 + \frac{m\gamma}{r^2}\left(D[SD(r\bar{V})] - 3\,S\,\zeta(r)\right) - 2\,k^2\,S\,\bar{\Omega}\,\zeta(r)\right\}} u = 0 , \end{equation} \tag{ 8 }$$

where

$$\begin{align*} \displaystyle{D} & = \displaystyle {\frac{d}{dr},\hspace{0.11in} D_{*} = \frac{d}{dr} + \frac{1}{r},\hspace{0.11in}\bar{\Omega} = \frac{\bar{V}}{r}, \hspace{0.1in} \hspace{0.11in} S = \frac{r^2}{m^2 + k^2r^2} \hspace{0.12in}} \nonumber\\[6pt] & \mbox{and} \hspace{0.12in} \nonumber\\[6pt] \displaystyle{\zeta(r)} & = \displaystyle{\frac{D(r\bar{V})}{r} = D_{*}\bar{V}.} \end{align*}$$

$\bar{\Omega}$ is the angular velocity of the mean flow, $\zeta(r)$ is its axial vorticity and

$$\begin{equation*} \gamma (r) = m\,\bar{\Omega} - \omega . \end{equation*}$$

At this point, the numerical procedure to solve the eigenvalue problem will only be outlined, because it will be seen that the axisymmetric case requires a different method from that used to deal with the helical case $m \geqslant 1$. A more detailed presentation will be given in section 3. In all cases, though, the range of integration is from r = 0 to $r \rightarrow \infty$. There is a regular singular point at r = 0, whereas the singularity at infinity is irregular. Series solutions will be employed to begin the integration at either end and the solutions matched at an interior point by requiring the Wronskian to vanish.

Near the origin, the solution can be represented by a Frobenius expansion having the form

$$\begin{equation} u = u_{0}\, r^{|m|-1}\ \left[ 1 + a_2 r^2 + O\left(r^4\right)\right], \end{equation} \tag{ 9 }$$

so that the radial perturbation velocity vanishes at the center of the vortex except in the case $|m| = 1$, the so called bending mode, and then it is continuous. The coefficient a₂ of the r² term is given by

$$\begin{align*} a_{2} & = \left\{(\delta + 3)(\delta + 1)\left(\frac{k \gamma_{0}}{m}\right)^2 - 4\left[k (1 + \varepsilon)\right]^2 \right. \nonumber\\ & \qquad \left. - 4 m \gamma_{0} \varepsilon_{2}\left(1 + \frac{k^2}{m^2}\right) - \varepsilon\frac{k^2}{m^2}\right\}/ 4 (2 + \delta) \gamma_{0}^2, \end{align*}$$

where $\delta = |m| - 1, \gamma_{0} = m (1 + \varepsilon) - \omega$ and $\varepsilon_{2} = 1 + 2\,\varepsilon$.

Far from the center of the vortex, the velocity profile can be be approximated as a potential vortex so that $\bar{V} \sim r^{-1}$. Employing that approximation, it can be shown that w, the axial velocity perturbation, satisfies a modified Bessel equation and K_m is the solution that decays as $r \rightarrow \infty$. From the linearized equations of motion, $u \sim w^{^{\prime}}$ and using the asymptotic expansion for $K_{m}^{^{\prime}}$, we find that for large r

$$\begin{equation} u \sim u_{\infty}\,\frac{e^{-k r}}{\sqrt{k r}}\ \Big ( 1 + \frac{4 m^2 + 3}{8 k r} +\cdots \Big ). \end{equation} \tag{ 10 }$$

A derivation of the far field solution is given in Ash and Khorrami (1995), but their asymptotic solutions (8.2.130b) are incorrect because they are missing the square roots.

2.1. Theory for axisymmetric modes

Setting m = 0 greatly simplifies the Howard–Gupta equation (8) and if, in addition, we let $g = u r$, the following Sturm–Liouville form results:

$$\begin{equation} \frac{d}{dr}\left(\frac{1}{r}\, \frac{dg}{dr}\right) + \left( \lambda\,\frac{\Phi(r)}{r} - \frac{k^2}{r}\right)\,g = 0, \end{equation} \tag{ 11 }$$

where $\lambda = k^2/\omega^2$ and $\Phi(r)$ is the Rayleigh discriminant. The boundary conditions are homogeneous with $g \sim r^2$ near r = 0 and $g \rightarrow 0$ exponentially as $r \rightarrow \infty$. Before discussing the stability conditions with respect to the Rayleigh discriminant $\Phi (r)$, we will prove the following necessary condition for instability involving the complex constant ω, whose imaginary part is the amplification factor.

Multiplying (11) by $r\,\omega^2\,g^{*}$ and integrating between the boundaries, we obtain for the imaginary part of the resulting equation

$$\begin{equation} - 2\, \omega_{i}\,\omega_{r}\, \int_{0}^{\infty}\,\left[\,\left( \frac{r + 1}{r}\right) \left|\frac{d g}{d r}\right|^2\, + k^2\, \left| g \right|^2 \right]d r = 0, \end{equation} \tag{ 12 }$$

where the homogeneous boundary conditions have been employed. Assuming that there is instability, i.e. $\omega_\mathrm{i} \gt 0$, (12) can be satisfied only if $\omega_\mathrm{r} = 0$. We will see that this is true only for axisymmetric modes and is generally not the case for the helical modes.

Having established now that $\omega = i\,\omega_\mathrm{i}$ for axisymmetric perturbations, we can see that if k is $\textit{O} (1)$, $\lambda = - k^2/\omega_\mathrm{i}^2$ will be a large negative number that becomes infinite as $\omega_\mathrm{i} \rightarrow 0$. It is therefore of interest to consider the product $\lambda\,\Phi(r)$ appearing in (11) and figure 2 below shows $\Phi(r)$. Loosely speaking, $\lambda\,\Phi(r)$ must be positive for some range of r in order that the eigenfunction will be oscillatory there, thus permitting the homogeneous boundary conditions to be satisfied. We see in figure 2 that in the region where Φ is negative that its magnitude is small. It will not be surprising then to see in the results below that λ must be very large for instability of any magnitude to occur. The numerical consequences of this will now be discussed.

Zoom In Zoom Out Reset image size

2.2. Computational issues for axisymmetric modes

It will be seen below that the eigenfunctions grow very rapidly as we integrate outward from r = 0. That being the case, it is easier to solve the equation for u than to solve (11) for g, because near the center of the vortex $u \sim r$, whereas $g \sim r^2$. The equation satisfied by u in the axisymmetric case m = 0 is the following:

$$\begin{equation} r^2\,\frac{d^2 u}{d \, r^2} + r\, \frac{d u}{d r} - \left[\left(k^2 - \lambda \Phi\right)r^2 + 1 \right] u = 0. \end{equation} \tag{ 13 }$$

Near the origin, it is obvious from (3) that the Taylor series expansion of Φ is of the form

$$\begin{equation} \Phi \approx \Phi_{0} + r^2 \Phi_{2} + \cdots. \end{equation} \tag{ 14 }$$

It follows that to lowest order in r, u satisfies the modified Bessel equation. Only one solution is bounded at r = 0, namely,

$$\begin{equation} u = I_{1} \left(\sqrt{k^2 - \Phi_{0} \lambda}\,r\right) = I_{1} \left(\frac{k}{\omega_\mathrm{i}} \sqrt{\Phi_{0} + \omega_\mathrm{i}^2}\,r\right). \end{equation} \tag{ 15 }$$

For large values of its argument, I₁ grows exponentially and that will be the case here because k is typically $\textit{O} (1)$ and $\omega_\mathrm{i} \ll 1$. Only when r is very small does u grow linearly in r and the Frobenius series below can be used to initiate the integration.

Unlike (8), the governing equation for helical modes, there are no singular points in (13) apart from the boundaries, so a straightforward Runge–Kutta procedure can be used to carry out the integration. Taking account of the rapid growth of u discussed above, we can utilize a power series to begin the integration outward towards a value of r where a matching condition to the outer solution is imposed. For the example used to illustrate the situation, the integration was initiated at r = 0.01 and a small step size, again equal to 0.01 was employed. Near the vortex center, we find that the solution can be written

$$\begin{equation} u = u_{0}\,r\, \left[ 1 + a_2 r^2 + a_4 r^4 + O\left(r^6\right)\right], \end{equation} \tag{ 16 }$$

where

$$\begin{align*} a_2 & = - \frac{1}{8} \left[4 \lambda\left(1 + \varepsilon\right)^2 - k^2 \right] \hspace{0.15in} \\[6pt] & \mbox{and} \hspace{0.15in}\\[6pt] a_4 & = \frac{\lambda}{2} (1 + \varepsilon)\left(\frac{1}{2} + \varepsilon\right) + \frac{1}{192}\left[4 \lambda (1 + \varepsilon)^2 - \left.k^2\right)\right]^2. \end{align*}$$

The asymptotic far field boundary condition for the outer solution can be obtained from equation (10) by simply setting m = 0. Given that $k r \gg 1$, it is sufficient to approximate u by the leading order term in (10). The second condition required to initiate the inward integration of (13) is obtained by differentiating the asymptotic form for u, which leads to

$$\begin{equation} u^{^{\prime}}(r) \sim \sqrt{k}\left( 1 - 2 \sqrt{k r} \right )\, \frac{u}{2}. \end{equation} \tag{ 17 }$$

To solve the eigenvalue problem, we first integrate outward from r = 0.01 and then integrate inward from r = 6 to an interior matching point that can be chosen by trial and error. It is most convenient to fix ε and k while iterating on λ until the Wronskian of the two solutions vanishes at the chosen interior point. Either arbitrary constant can then be multiplied by an appropriate factor to arrive at an eigenfunction that is continuous.

The shooting method described above was usually successful for values of $|\lambda| \lt 100$. However, it was difficult to obtain a converged solution for larger values of $|\lambda|$. The eigenfunction in figure 3 below illustrates the difficulty, i.e. the rapid growth of u as we integrate away from r = 0. For this case, $\varepsilon = 0.90, k = 0.70$ and the solvability condition of vanishing Wronskian was imposed at $r_{0} = 0.86$. A converged solution was obtained for $\lambda = - 103.71$, which corresponds to an amplification factor $\omega_\mathrm{i} = 0.069$.

Zoom In Zoom Out Reset image size

As predicted by (15), the solution for u begins to grow exponentially not far from the vortex center. Moreover, the numerical solution shows that this extremely rapid growth continues even beyond r = 1. Clearly, our shooting method needs to be modified in order to find solutions when $|\lambda|$ is very large. Such a method will be described in the following section. We note, in passing, that similar problems were encountered by Smyth and McWilliams (1998) in their investigation of the linear instability of a vortex in the context of geophysical fluid dynamics. Specifically, they were unable to treat unstable modes, known to be present, when the density stratification was large compared to the planetary rotation.

As suggested above, the numerical method employed in the axisymmetric case must be different from that used to deal with helical modes. The two principal differences are consequences of the result $\omega_\mathrm{r} = 0$ when m = 0. At first sight, the helical case would appear to be more difficult, because an indented contour is required to deal with the critical point singularity and the eigenvalue iteration involves two constants (which we will choose to be $\omega_\mathrm{r}$ and $\omega_\mathrm{i}$), rather than just $\omega_\mathrm{i}$. The greater challenge, however, is with the axisymmetric problem, owing to the magnitude of the eigenfunction u(r), as illustrated in figure 3. Even with the arbitrary constant $u_{0} = 2 \times 10^{-4}$, the eigenfunction u reaches a maximum value slightly larger than $91 \times 10^{6}$. A numerical method will now be presented that is able to treat this situation.

3.1. Numerical method and results for axisymmetric modes

It was pointed out below (15) that there is a region not far from the origin where u grows exponentially and an example was presented to illustrate the consequences of this growth. Given that $u^{^{\prime}}$ also grows exponentially at the same rate, we can define a new dependent variable proportional to $u^{^{\prime}}/u$ that should be better behaved. This method, used previously to treat the stability of two-dimensional free shear layers (Michalke 1965, Maslowe and Kelly 1971), leads to a nonlinear first order equation and is known as a Riccati transformation. An apparent complication in the axisymmetric case, however, is that $u^{^{\prime}}/u \sim 1/r$; this can be seen from the Frobenius expansion (16). To remove the singularity at r = 0, we modify the usual transformation by writing

$$\begin{equation} u = r\,e^{\textstyle \int_{0}^{r} Q(r)\,d r}. \end{equation} \tag{ 18 }$$

Substituting now the change of dependent variable defined by (18) into the governing equation for u, i.e. (13), leads to the Riccati equation

$$\begin{equation} \frac{d Q}{d r} + \frac{3 Q}{r} = - Q^2 - (k^2 - \lambda \Phi ). \end{equation} \tag{ 19 }$$

The procedure employed to solve the eigenvalue problem associated with (19) involves an inward and an outward integration, as before. Instead of requiring the Wronskian of the two solutions to vanish at the matching point, however, we now require that Q be continuous.

There are several possibilities to begin the outward integration, details of which can be deduced by first differentiating (18) and then introducing the Frobenius expansion (16) to approximate $u^{^{\prime}}/u$. We obtain

$$\begin{equation} Q = \frac{u^{^{\prime}}}{u} - \frac{1}{r} \approx 2\, a_{2}\, r + 2\, ( 2\, a_{4} - a_{2}^2 )\, r^3 + \cdots , \end{equation} \tag{ 20 }$$

where a₂ and a₄ are given below (16). In most cases, the outward integration was initiated at r = 0.02 with the initial value of Q provided by (20). Noting that $Q(0) = 0$, (19) could also be linearized around r = 0 and the equation integrated. If Φ is expanded as in (14), then the first derivative at r = 0 is given by

$$\begin{equation} \frac{d Q}{d r} (0) = (k^2 - \lambda \Phi_{0})\,\frac{1}{4} \end{equation} \tag{ 21 }$$

and the integration could begin at r = 0. This, of course, is exactly the result obtained by differentiating (20) and setting r to zero.

Figure 4 below presents the growth rates as a function of axial wavenumber for three different values of ε. We see that the instability is substantial and it increases monotonically with k in all cases. The qualitative behavior of ω_i as a function of k is similar to that reported by Billant and Gallaire (2005) for a family of zero circulation vortices relevant in geophysical applications with velocity profile $\bar{V} = r \exp(- r^{\alpha})$. Such behavior, termed 'ultraviolet catastrophy' in the vortex study of Smyth and McWilliams (1998), makes it clear that viscosity should be incorporated to obtain the correct behavior for short waves.

Zoom In Zoom Out Reset image size

To obtain the results in figure 4, we had to deal with values of $|\lambda|$ as large as 2000. Before employing the Riccati method to investigate the stability of the perturbed Lamb–Oseen vortex, the program was tested on the case $\bar{V} = r \exp(- r^{2})$ and the results seemed to agree well with those in figure 1(a) of Billant and Gallaire (2005). A Chebyshev collocation method was employed in the latter study. We do not know if the authors chose $\omega_\mathrm{r} = 0$ from Sturm–Liouville theory or if they determined it by iteration, because nothing is said in their paper about $\omega_\mathrm{r}$. It will be seen below that $\omega_\mathrm{r} \neq 0$ in the case of nonaxisymmetric modes, so the question of its behavior should not be ignored.

The Riccati method was also used by Michalke and Timme (1967) to investigate the stability of the velocity profile (5). Their computations, however, did not face the same challenges as those in the present paper, because the vortex described by (5) is stable to axisymmetric perturbations. As noted in the introduction, the focus of their paper was on unstable nonaxisymmetric disturbances, so comparisons will be deferred to the following sub-section.

3.2. Numerical method and results for helical modes

We present in this section results obtained by solving the Howard–Gupta equation (8) for non-axisymmetric perturbations with m = 1 or m = 2. It will be seen that the growth rates decrease rapidly with increasing m, so there is little point in considering values of m > 2. The numerical procedure used will essentially be the shooting method described above in section 2.2. However, u and ω are now complex, so that the integration contour must be indented to avoid the critical point singularity, where $\gamma (r_\mathrm c) = 0$ (recall that $\gamma (r) = m\,\bar{\Omega} (r) - \omega$). The viscous critical layer analysis of Le Dizès (2004) shows that the correct eigenvalue condition for a neutral mode is obtained by considering the limit $\omega_\mathrm i \downarrow 0$. For our vortex, this means that the integration contour must pass above the singularity, a consequence of $\gamma_\mathrm{c}^{^{\prime}}$ being negative.

The stability calculations for the m = 1 mode are of particular interest for several reasons. This mode is distinctive because it is the only mode with nonzero radial velocity at the vortex center. That characteristic makes its presence clearly visible in the elliptical instability experiments of Leweke and Williamson (1998). Although most related studies indicate that other values of $|m|$ lead to larger growth rates in the inviscid case, Mayer and Powell (1992) present some results for the Batchelor vortex suggesting that this may not be the case when viscosity is included.

In figure 5 below, results are presented for both the real and imaginary parts of ω as a function of ε with the wavenumber fixed at the value k = 1.50. We see from figure 5(a) that unlike the axisymmetric mode, the frequency ω_r is nonzero and it increases quite rapidly with ε. With regard to the amplification rates, at the chosen value of k, the bending mode is damped very slightly until ε reaches about 0.22. Although it increases rapidly beyond this value of $\varepsilon, \omega_\mathrm{i}$ is still smaller than it is for the m = 0 case. For example, at $\varepsilon = 0.80, \omega_\mathrm{i} = 0.0642$ for the bending mode compared with 0.0960 for the axisymmetric disturbance. Given that the wavenumber corresponding to the largest amplification rate increases very slowly with ε, figure 5(b) depicts quite closely the maximum value of $\omega_\mathrm{i}$ as a function of ε.

Zoom In Zoom Out Reset image size

Continuing now with the behavior of the bending mode, we examine the eigenfunctions and show that they do not resemble those of the axisymmetric case. They are, of course, complex and are nonzero at the vortex center. In addition, they oscillate slightly, particularly for $\textit{O}(1)$ values of r. Although the author has examined eigenfunctions for many different axial wavenumbers and even some for higher modes, the one illustrated in figure 6 is sufficient to exhibit these differences compared with the axisymmetric mode shown in figure 3. The magnitude of $|u|$ is directly related to the term $k^2/(m\,\bar{\Omega} - \omega )^2$, which multiplies $\bar{\Omega} \zeta$ in (8). This term is much smaller in the case of the bending mode, as a direct result of the fact that m and $\omega_\mathrm{r}$ are nonzero when $m \geqslant 1$. In the axisymmetric case, $k^2/(m\,\bar{\Omega} - \omega )^2$ becomes equal to $-\lambda$, which is very large when k is large and/or $\omega_\mathrm{i}$ is small.

Zoom In Zoom Out Reset image size

To conclude our presentation of results on the stability of the perturbed Lamb–Oseen vortex to helical perturbations, we compare growth rates for the m = 2 mode with those of the bending mode m = 1. Figure 7 presents a comparison of amplification rates for these two modes at a value of ε = 0.6. It is, perhaps, surprising that the amplification rate of the m = 2 mode is so much smaller. However, this is consistent with results in figure 1(a) of Billant and Gallaire (2005) for the vortex $\bar{V} = r \exp(- r^{2})$. A qualitative difference that should be pointed out, however, is that both stability curves in figure 7 contain a value of k where the growth rate is maximum. This is, of course, normally the case, but it is not true for the vortex $\bar{V} = r \exp(- r^{2}),$ where ω_i continues its monotonic growth with k. A related point of some interest is the behavior of these stability curves at what we can term the critical values of k. As k increases near 1.04, the m = 1 mode goes from being damped to becoming amplified. Near the other end, however, the stability curve never intersects the k-axis and this mode becomes a Kelvin wave. The same behavior is observed for the m = 2 helical mode.

Zoom In Zoom Out Reset image size

Finally, we compare the results obtained in this section with those of Michalke and Timme (1967), who investigated the velocity profile (5). As noted in the Introduction, this velocity profile is stable to axisymmetric perturbations. However, the possibility of instability to non-axisymmetric perturbations is suggested by (21) of Howard and Gupta (1962). In the absence of axial flow, the latter result can be written

$$\begin{equation} k^2\,\Phi - \frac{1}{4}\,m^2\,\bar{\Omega}^{^{\prime} 2} \geqslant 0 \end{equation} \tag{ 22 }$$

is a sufficient condition for stability. Clearly, for helical modes there will always be a value of k small enough that this condition is violated. That does not guarantee instability, however, and computations must be carried out for specific cases to see if and when instability occurs. That said, let us now briefly summarize the results presented by Michalke and Timme for the continuous profile (5).

First of all, the behavior of $\omega_\mathrm{i}$ with k is the exact opposite of that found for either the perturbed Lamb–Oseen vortex or the profile $\bar{V} = r \exp(- r^{2})$. The largest growth rates were obtained at k = 0 and ω_i was found to decrease rapidly with increasing k. Secondly, the azimuthal wavenumber leading to the largest growth rates was m = 4; for m = 3 and m = 5, the maximum value of $\omega_\mathrm{i}$ was approximately 25% smaller and for m = 1 only a neutral mode was found. Possibly, these results are related to the behavior near r = 0 of the model employed in Michalke and Timme (1967), i.e. $\bar{V}\sim \frac{2}{3}\,r^5$, as opposed to the rigid rotation characteristic of most vortices. (We note, in passing, that the ordinate of figure 4 in their paper is labeled incorrectly. This figure shows the growth rate for a discontinuous model, so it is not directly related to this paper.)

3.3. Extension of Rayleigh's criterion to helical modes

As cited in the Introduction, an approximate generalization of Rayleigh's criterion for axisymmetric disturbances was derived by Billant and Gallaire (2005). Their development assumes that k is large and the WKB method was utilized. The eigenfunction in their analysis was assumed to have two turning points, implying an oscillatory region along the real axis between the Stokes lines emanating from these turning points. Outside of the oscillatory region, the eigenfunction was represented by the usual exponentially growing and decaying WKB functions. In the large k limit, it was ascertained by Billant and Gallaire (2005) that the onset of instability is determined by first finding the value $r = r_{0}$, where $\partial \omega/\partial r = 0.$ The desired value of r₀ is obtained by solving the implicit equation

$$\begin{equation} \Phi^{^{\prime}}(r_{0}) = - 2\, i\,m \,\bar{\Omega}^{^{\prime}}(r_{0})\,\sqrt{-\Phi (r_{0})}. \end{equation} \tag{ 23 }$$

A generalized Newton iteration was used to solve (23) for r₀ and it was found to converge rapidly even when the initial guess was inaccurate. Once r₀ has been obtained, ω, the complex frequency, is given by the equation

$$\begin{equation} \omega (r_{0}) = m\,\bar{\Omega}(r_{0}) + i\,\sqrt{-\Phi (r_{0})}. \end{equation} \tag{ 24 }$$

Some results for azimuthal wavenumbers m = 1 and m = 2 are presented in table 1. We see that the onset of instability is predicted by (24) to occur at ε ≈ 0.22 for the m = 1 bending mode and ε ≈ 0.51 for the m = 2 helical mode. Using the method described in section 3.2, the actual computed results for the critical values of ε were found to be 0.21 for m = 1 and ε = 0.52 for m = 2. In both cases, the predictions based on large-k asymptotics are seen to be remarkably accurate. This is less surprising in the m = 2 case, because the corresponding value of k is 3.28 and the expansion is in powers of k⁻². However, for m = 1 the wavenumber at the stability boundary was k = 1.5, so the method seems to work well even for waves that are not extremely short. The eigenfunctions (see e.g. figure 6) computed by the author do not resemble the configuration assumed by Billant and Gallaire, so we must conclude that their method does not require a faithful representation for u(r) in order to yield accurate predictions for the onset of instability.

Table 1. Frequency and amplification factors for helical modes m = 1 and m = 2 near the critical values of ε.

m	ε	r0	ωr	ωi
1	0.20	$2.4081 + 0.00460 i$	0.1633	−0.000035
1	0.22	$2.2833 + 0.00004 i$	0.1745	$-4.6\times10^{-9}$
1	0.24	$2.1815 + 0.08190 i$	0.1836	0.00101
2	0.50	$1.5154 + 0.0009 i$	0.6338	$- 3.09\times10^{-6}$
2	0.51	$1.5350 + 0.0060 i$	0.6402	$9.87 \times10^{-6}$
2	0.52	$1.5315 + 0.0503 i$	0.6468	0.00036
2	0.54	$1.5170 + 0.0815 i$	0.6598	0.0016

In the introduction of this paper, it was pointed out that despite an extensive literature dealing with the stability of vortices, very few articles investigate the instability of a single vortex with no axial flow. Most studies of wave-like perturbations on vortices have instead concerned the propagation of Kelvin modes. An example is the discontinuous Rankine vortex comprised of an inner region of rigid rotation and an outer region, where the behavior is that of a potential vortex $\bar{V}\sim r^{-1}$. In the case of the more realistic Lamb–Oseen vortex, some new modes exist owing to the presence of a critical layer. However, these modes turn out to be damped, which is disappointing from the standpoint of hydrodynamic stability theory.

It appears that the paper by Michalke and Timme (1967) was the first to actually investigate numerically the instability of a continuous, unbounded vortex. As discussed in section 3.2, the approach of Michalke and Timme was somewhat different from our own in that their goal was to demonstrate that a vortex, stable according to Rayleigh's circulation theorem, could be destabilized by a non-axisymmetric disturbance. As an example of a continuous model, they showed that the vortex with velocity profile (5) was unstable to helical perturbations having azimuthal wavenumbers $2 \leqslant m \leqslant 6$. In addition to being of interest in its own right, this result shows that the analysis in section 4 of Billant and Gallaire (2005) is sometimes incorrect. The latter purported to demonstrate that axisymmetric modes always have higher growth rates than helical perturbations, which is clearly not the case for the vortex profile (5).

In this paper, the author has shown that a slight perturbation of the stable Lamb–Oseen vortex can give rise to instabilities of significant growth rate. The results presented in sections 3.1 and 3.2 show that the growth rate $\omega_\mathrm{i}$ of these disturbances becomes large when $\varepsilon,$ the perturbation parameter, is increased beyond infinitesimal values. This is particularly true for axisymmetric modes as k, the axial wavenumber, becomes large. The monotonic increase in growth rate with k shown in figure 4 was previously observed in the case of the velocity profile $\bar{V} = r \exp(- r^{2})$, as reported by Billant and Gallaire (2005). For the helical modes of azimuthal wavenumber m = 1 and m = 2, on the other hand, ω_i behaves in a more conventional way in the case of the perturbed Lamb–Oseen vortex. As shown in figure 7, there is a finite range of unstable wavenumbers, whereas the growth of $\omega_\mathrm{i}$ continues to be monotonic with increasing k for the vortex $\bar{V} = r \exp(- r^{2})$.

Even though the hydrodynamic stability characteristics of a single vortex are known for only a few cases, we can still make some general observations on what has been learned up to this point. First of all, the largest growth rates are usually associated with axisymmetric modes, where there is a monotonic increase of $\omega_\mathrm{i}$ with k in all cases of a vortex whose angular velocity is nonzero at its center. This suggests that centrifugal instability is the dominant mechanism, as opposed to shear instability. That being the case, one can assume that the velocity of the vortex at the value of r where Φ changes sign is related to the growth rate. For the perturbed Lamb–Oseen vortex, equation (3) shows that Φ, the Rayleigh discriminant, becomes negative for $r_\mathrm{c} \gt \sqrt{1 + 1/\varepsilon}$. As ε increases from zero, $r_\mathrm{c}$ moves in from infinity, $\bar{V} (r_\mathrm{c})$ increases and so does the amplification factor of the instability, as shown in figure 4.

The occurrence of a maximum growth rate as $k\rightarrow\infty$ is associated with centrifugal instability and contrasts totally with the behavior of the non-axisymmetric modes for the vortices investigated by Michalke and Timme (1967). The latter authors studied three vortex models, two of which had discontinuous vorticity profiles, and in all 3 cases the largest growth rates were obtained for k = 0. This behavior suggests a shear layer mode of instability. The one instance where similar behavior was observed in Billant and Gallaire (2005) was for the case of a vortex with velocity profile $\bar{V} = r \exp(- r^{4})$, where the growth rate was nonzero at k = 0 when m = 2 and m = 3. The monotonic growth of ω_i with k was found, however, to still occur as $k \rightarrow \infty$.

From the results obtained to date, it is obvious that numerical solution of the Howard–Gupta equation (8) is required to determine the stability properties of a single vortex having a continuous velocity profile. Even in cases where the Rayleigh criterion guaranties instability one does not know the value of the axial wavenumber k corresponding to the onset of the instability. As for the value of k yielding the maximum growth rate, the examples reviewed above reveal that any value is possible, including k = 0. The results presented in figures 4–7 of this paper suggest a number of potentially interesting extensions. The most obvious would be to see if the addition of a slight viscosity would damp the short-wave instability evident in the axisymmetric mode shown in figure 4. Other effects whose incorporation would be of interest could include density stratification and/or nonlinearity.

Finally, there are important applications in geophysical fluid dynamics, such as tornadoes, that were mentioned briefly in this paper. Gent and McWilliams (1986) argued that vortices with velocity profiles behaving as potential vortices in the far field decay too slowly to accurately simulate these phenomena. Instead, they investigated several models with velocity profiles such as $\bar{V} = r \exp(- r^{\alpha})$ that have zero total circulation. Given that computational issues limited the extent of the results reported in both Gent and McWilliams (1986) and Smyth and McWilliams (1998), it would be interesting to expand those results using a more robust scheme; for example, the Riccati method described in section 3.1.

The author is pleased to acknowledge the assistance of several colleagues in the preparation of this article. These include Peter Bartello, David Jeffrey, Bruce Yin and Gantumur Tsogtgerel.

There are many factors, in addition to integration step size, that influence the accuracy of computations of the sort presented above. At the suggestion of a referee, this appendix has been added to the paper in order to summarize the tests made to verify our results. We begin by a comparison with a result for a damped $m = k = 1$ mode on a Lamb–Oseen vortex reported by Fabre et al (2006). Table 1 of their paper is intended to illustrate how, as the Reynolds number Re is increased, the solution tends towards the inviscid value given in the table caption as $\omega = 0.15324 - 0.02067 i$. Using the indented contour method described in section 3.2, the author obtains $\omega = 0.15324 - 0.02267 i$. Given the perfect agreement for the frequency $\omega_\mathrm{r}$, it was surprising to note the significant disagreement for $\omega_\mathrm{i}$. Fortunately, the values provided in table 1 of Fabre et al for Re$= 10^3, 10^4$ and 10⁵ make it clear that there is a typographical error in their inviscid result for $\omega_\mathrm{i}$. In other words, once that error is corrected, there is actually perfect agreement to five significant figures.

A similar comparison involving a non-axisymmetric m = 2 mode was not entirely possible because only results for Re$= 10^3$ were provided by Fabre et al (2006). It was nonetheless interesting to compare an inviscid computation with their result in figure 21(d) for k = 5, where it was found that $\omega = 0.4077 - 0.0730 i$. This compares with our inviscid result $\omega = 0.4177 - 0.0450 i$. The same tendencies are present as in the case of the bending mode with a somewhat greater damping effect of viscosity. This is due to the axial wavenumber k = 5 being much greater than k = 1, as viscosity has a greater effect on short waves.

Our final comparison will be with results reported by Smyth and McWilliams (1998) and Billant and Gallaire (2005) for the vortex with velocity profile $\bar{V} = r \exp(- r^{. 2})$. Given that the aforementioned papers contain figures, but no tabular results, the comparison must be qualitative. As shown in figure A1, the monotonic growth of ω_i with k for large values of k is similar to that illustrated in figure 1(a) of Billant and Gallaire (2005) and figure 4 of Smyth and McWilliams (1998). And, we recall that figure 4 of this paper for the case of the perturbed Lamb–Oseen vortex displays the same behavior.

Zoom In Zoom Out Reset image size

The computations shown in figure A1 were carried out using the Riccati transformation presented in section 3.1. For a number of values of k, the governing eigenvalue problem was also solved using the conventional shooting method described in section 2.2. Comparing the values for $\omega_\mathrm{i}$ to 4 significant figures, it was found that the discrepancy was never greater than 0.0001. Taking k = 1.5, for example, the Riccati method yielded $\omega_\mathrm{i} = 0.13849$ compared with $\omega_\mathrm{i} = 0.13843$ when solving equation (13), the second order ODE.

The factors varied to test the accuracy of our computations included integration step size, starting point of the integration in either direction, number of terms in the Frobenius solution and upper bound of the Wronskian at the value of r where the solutions were matched. The integrations were carried out in double precision and the results turned out to be surprisingly insensitive to all of these factors. As an example, the step size was varied with k = 3 when integrating equation (13) outward, again with the velocity profile $\bar{V} = r \exp(- r^{2})$. The integration was normally initiated at r = 0.01 which was also the step size. As a test, the integration was initiated at r = 0.02 with a step size of $\delta r = 0.04$ and the solutions matched at r = 0.86. The Wronskian at the matching point was required to be less than 10⁻³ and the two solutions for $\omega_\mathrm{i}$ turned out to be identical to five significant figures!