Numerical realization of helical vortices: application to vortex instability

The need to numerically represent a free vortex system arises frequently in fundamental and applied research. Many possible techniques for realizing this vortex system exist but most tend to prioritize accuracy either inside or outside of the vortex core, which therefore makes them unsuitable for a stability analysis considering the entire flow field. In this article, a simple method is presented that is shown to yield an accurate representation of the flow inside and outside of the vortex core. The method is readily implemented in any incompressible Navier–Stokes solver using primitive variables and Cartesian coordinates. It can potentially be used to model a wide range of vortices but is here applied to the case of two helices, which is of renewed interest due to its relevance for wind turbines and helicopters. Three-dimensional stability analysis is performed in both a rotating and a translating frame of reference, which yield eigenvalue spectra that feature both mutual inductance and elliptic instabilities. Comparison of these spectra with available theoretical predictions is used to validate the proposed baseflow model, and new insights into the elliptic instability of curved Batchelor vortices are presented. Furthermore, it is shown that the instabilities in the rotating and the translating reference frames have the same structure and growth rate, but different frequency. A relation between these frequencies is provided.


Introduction
The motion of an infinite helical vortex is one of the classical problems in applied mathematics and fluid mechanics that was perhaps first time addressed by Levy and Forsdyke [1]. Such a vortex moves through the fluid without changing its shape [2] and is an idealization of a common flow structure. Related to this problem is that of determining the corresponding velocity field that a helical vortex induces in a certain domain D. This is of prime interest during numerical simulations and is a prerequisite for stability calculations.
The velocity induced by a given vorticity distribution may be obtained from the Biot-Savart law, which by assuming that the vorticity is concentrated to a line vortex is expressed as where x is some point in D in which the velocity is to be evaluated, C is the line vortex parameterized by s, t is the tangent vector to the line, Γ is the vortex strength and x denotes a point on C [3] ( · refers to the Euclidean norm). The analysis of a flow field induced by Eq. (1) is obstructed by the presence of the singularity in the integrand when evaluated for a point on the vortex line. One way of circumventing this issue is to introduce a cut-off length δ and exclude a small region x − x < δ from the integration [4]. Another approach that was suggested by Rosenhead [5] and elaborated upon by Moore [6] is to desingularize the integral by modifying the denominator in (1) to ( x − x 2 + μ 2 ) 3/2 . The parameters δ or μ can be chosen such that the induced velocity at a point on the vortex matches that of a vortex ring for which an analytical expression is known [7]. Although these methods will give an accurate representation of the induced velocity far away from the vortex, they are ad hoc and will provide an ill-defined flow field in the vicinity of the vortex core, which might be the most important region from a stability point of view. To overcome this issue, two approximate solutions valid inside and outside the core, respectively, of a perturbed straight vortex were derived by [8] using matched asymptotic expansions. Later, Hardin [9] presented analytical expressions for the velocity field in the interior and exterior of a helix, which were furthered by Okulov [10]. These expressions however provide no details of the flow inside the vortex core. Elaborate expressions for this region have instead been presented by Callegari and Ting [11] and by Blanco-Rodríguez et al. [12]. Among the more recent contributions, we also mention the wake model proposed in [13], where the velocity field is evaluated through an asymptotic expansion of the Biot-Savart law after the vortex line has been determined iteratively.
The main motivation of this work is the desire to accurately represent a system consisting of M helical vortices numerically in order to perform a global stability analysis of this flow. A method of generating helical vortices that have gained much popularity in the wind turbine community is the actuator line (ACL) method [14]. Despite the success of the ACL method, it requires measured blade data to be available, which limits its application and makes it difficult to implement for an arbitrary helix. It also generates a root vortex, which is an inevitable feature in most real-world flows, but can be undesirable in a theoretical study. Another issue is the requirement of a Gaussian kernel to distribute the forcing from the actuator lines onto the grid points of the mesh, which may greatly influence the results and thus must be judiciously tuned [15].
A highly accurate baseflow is required when performing a numerical stability analysis. For the case of an open flow with a prescribed (helical) vortex, this implies that a detailed flow representation is required both inside the vortex core and in the far-field. Somewhat surprisingly, the above literature review has led us to conclude that a general procedure for generating such a field is missing. In this article, we therefore propose a step-by-step algorithm to generate a baseflow that is simple enough to potentially be implemented in a wide range of solvers, yet yields sufficiently accurate results to enable a meaningful stability analysis to be carried out. Our method is in certain aspects reminiscent of the technique recently discussed in [16], although more general. Selçuk et al. [16] considers a specific formulation of the Navier-Stokes equations that imposes helical symmetry and thereby reduces the three-dimensional flow problem to a two-dimensional one. Our method is, however, derived in the context of primitive variables in a Cartesian coordinate system. As a result, any vortex geometry and core structure is in principle possible in our framework. For the sake of brevity we will however limit the discussion to the case of two interlaced helical vortices. In spite of a substantial interest for this geometry within applications (e.g. rotor wakes etc.), there are still numerous open questions related to the transition of this flow, such as the competition between different types of instabilities and the appearance of these instabilities in different frames of reference, which motivate its further analysis. It is also a flow for which theoretical predictions of the eigenvalues and eigenvectors are available, which provides a direct means of validating the modeling approach. A recent experimental study by Quaranta et al. [17] has provided detailed measurements for the configuration at hand. This ensures that the flow case investigated is physically meaningful.

Baseflow generation
The fluid motion may in general be described as a combination of an isotropic expansion, a volume preserving straining motion, and a rigid-body rotation [3]. In what follows, we will consider a flow field U that is assumed to be solenoidal (i.e. have zero rate of volume expansion) and consist of a set of M helical vortices that are embedded in a specified background flow. To determine U, the flow will be decomposed into an irrotational part U s and a rotational part U v such that U = U s + U v .

Vortex-induced flow component
We consider first the rotational part and assume that this satisfies where ω is a specified vorticity field. Depending on the frame of reference, it might be useful to express the vorticity field as the sum of a background field ω b and that due to a given vortex ω f , such that In this section, we focus on U v f and defer a discussion about U v b to the next. In order to determine U v f , the flow can be expressed as a vector potential, Assuming that ∇ · B v f = 0, a Poisson problem is obtained for the vector potential which has the formal solution in the domain D. This solution can be shown to be divergence-free if the vorticity field is either localized or tangential to the boundary [3].
In order to compute B v f , the global vorticity field ω f induced by the M-multiple of helical vortices must be determined. The center line of the jth individual vortex is a curve in R 3 parameterized by s and given by where ψ j = 2π( j − 1)/M is the angle separating the helices in the azimuthal direction ( j ∈ {1, . . . , M}), and R and L is the radius and the reduced pitch of the helices, respectively. The reduced pitch is related to the helical pitch, i.e. the streamwise separation h between two consecutive helix loops (of the same helix) as L = h/(2π), and the angle φ is related to s by where γ = √ R 2 + L 2 is introduced to simplify the notation. A right-handed orthogonal coordinate system that is aligned with the helix line C j for every s is provided by the Frenet-Serret formulae [18] is the torsion, curvature, and radius of curvature, respectively. Considering an axisymmetric vortex, the vorticity in a plane P(s, j) orthogonal to the local axial vector t j (s) and spanned by the local normal n j (s) and binormal vector b j (s), is expressed asω(r ), where r is a local radial coordinate (see Fig. 1). (The tilde-sign is used to denote variable quantities in the local coordinate systems.) Fig. 1 Sketch of the vortex (green), the Cartesian coordinate system and the Frenet-Serret coordinate system at s = φ = 0 and j = 1. Shown is also the P(s, j)-plane, which is spanned by the normal n j (s) and the binormal b j (s) vectors to the helix center line C j (s), and is perpendicular to the tangential vector t j (s) In order to evaluate the vorticity in a point x = (x, y, z) T ∈ D due to the jth vortex, the appropriate P-plane containing this point must be determined. This is done by finding the shortest distance from x to C j , which amounts to solving the nonlinear problem for φ(s). Equation (11) may have several roots, in which case the solution minimizing the distance is retained. Even so, there may exist several φ that yield the same minimum distance. Although this could be a potential issue, it is usually of minor concern provided that the vortices are concentrated such that the vorticity decays quickly with distance from the vortex center line C j . Given s and φ(s), the global Cartesian coordinates are mapped to the local Frenet-Serret coordinate system ξ (s, j) = (ξ(s, j), η(s, j), β(s, j)) T by, and expressed in cylindrical coordinates r = (r, θ, ξ) T . (Note that ξ = 0 since x ∈ P(s, j).) Once in the local cylindrical coordinate system, the vorticity is evaluated. In this work, the Batchelor vortex [19] is considered, which models a trailing edge vortex with an axial core flow. The axisymmetric vorticity profiles ω = (ω r ,ω θ ,ω ξ ) T read where Γ is the circulation, U c is the axial center line velocity and a is the radius of the vortex core. The corresponding velocity profilesŨ r = (Ũ r ,Ũ θ ,Ũ ξ ) T are given as where U 0 is a constant. Subsequently, the vorticity is mapped back from the cylindrical to the Frenet-Serret coordinate system, and then to the Cartesian coordinate system,

Background flow component
The remaining flow components that need to be determined in order to complete the description of the baseflow are the irrotational component U s and the rotational background component U v b . The irrotational part of the velocity field has the properties This part can be written as the gradient of a scalar field (i.e. velocity potential) that satisfies the Laplace equation [3]. The properties of U v b are specified in (3c) and (3d). Both U s and U v b are dependent on the frame of reference in which the vortices are viewed, and must hence be specified accordingly. The motion of a helical filament can be described as a rotation in the φ-direction, a translation in the z-direction and a slipping along the filament [2]. Hence, two equilibrium frames exist in which the vortices will be steady (or quasi-steady), namely a rotating and a translating frame of reference. The inviscid motion of an equilibrium array of helical vortices has been considered by Okulov [10]. Such a system moves and rotates with the following translational and azimuthal velocities The angular velocity Ω is composed of the angular velocity induced by other vortices in the multiple, and the angular velocity induced by the vortex itself. In [10], the following expression for this quantity is derived where ζ(·) is the Riemann zeta function, = L/R is the dimensionless reduced pitch (introduced to simplify the notation) and a e is the vortex core size. Equation (18) assumes a constant distribution of vorticity throughout the core. In situations where more realistic vortices with a continuous vorticity profiles are considered, one can choose a e to be the equivalent core size of a Rankine vortex without axial flow [8,20]. With the aid of (17) and (18), the motion of the equilibrium frames can be determined. In the rotating frame of reference, which is common in numerical simulations [21,22], one has and (19b) are readily seen to satisfy (3c) and (3d). For the translating frame of reference on the other hand, one instead has The translating frame, although less common, is useful in many situations since it enables a straightforward comparison with experimental results.
It is noted that (17) with (18) only correspond to equilibrium velocities in the absence of viscosity. When viscosity is considered, the vortices will continuously diffuse, and hence one may expect the induced and self-induced velocities to change slightly over time. The timescale over which this attenuation takes place is however reasonably long, and therefore (17) and (18) may still yield a good estimate of the motion of the helices for a limited time horizon. This is indeed verified to be the case in the direct numerical simulations presented in Sect. 2.4.

Validation and implementation
The described approach has been implemented in the high-order spectral element solver Nek5000 [24], which does not enforce helical symmetry but solves the Navier-Stokes equations in a Cartesian coordinate system using primitive variables in a P N − P N −2 formulation [25]. Homogeneous Dirichlet boundary conditions are imposed on the outer rim and the flow is periodic in the z-direction (only one period of the helix is generated). A cross section of the computational mesh is shown in Fig. 2a. Equation (5) is here preconditioned with an additive Schwarz preconditioner [26] and solved with GMRES [27]. To solve (11), a standard bisection method with interval bracketing is used [28].
The numerical values chosen are listed in Table 1 and correspond to the experimental configuration recently investigated by Quaranta et al. [17]. The parameters are made dimensionless with the circulation measured in the experiment Γ * = 68.8 cm 2 /s and the measured terminal helix radius R * = 8.7 cm. The Reynolds number is defined as Re = Γ * /ν * where ν * is the kinematic viscosity. In addition, the experiments had a free-stream velocity of U * ∞ = 37 cm/s, and the tip-speed ratio for this case was λ = 5.4. The resulting two helices are visualized in Fig. 2b.
In Fig. 3 vorticity and velocity profiles across the vortex core are extracted from one of the generated vortices (φ = 0, ψ j = 0) in the translating frame of reference, and compared against the theoretical profiles given by (13) and (14). As seen, the method successfully generates a helix with the prescribed vorticity profiles. Inside the support region of the vorticity profiles, the extracted velocities are also seen to closely match those of the Batchelor vortex, but outside of this region they differ slightly. This is due to the fact that the vortices induce a 2Γ / h velocity difference between the outside and the inside of the helix [29], which is not accounted for by (14). For the purpose of comparing the theoretical and extracted axial velocity profiles in Fig. 3d, the coefficient U 0 in (14c) is chosen to be the average difference between the theoretical and the numerical data.  (13)- (14), and data extracted from the field obtained by solving (5) (filled circles). The data are extracted along the normal axis in the P(0, 1)-plane

Vortex adaptation
Although the vortex visualized in 2b fulfills all the desired properties, it is not a solution to the incompressible Navier-Stokes equations We insert the derived field as an initial condition for (21) and advect the solution for a short period of time. To ensure an accurate time evolution, a third-order semi-implicit extrapolation/backwards differentiation (EXT3/BDF3) time-integration is used. The advection term is evaluated using a convective formulation that is stabilized using over-integration [30] and 1% of filtering [31] applied on the peripheral elements. The term F on the right-hand side of (21a) is a volume force that in the rotating frame of reference corresponds to the sum of the centrifugal and the Coriolis force, and in the translating frame of reference vanishes. (Note that the centrifugal force can be written as the gradient of a scalar field, i.e. −ω b × (ω b × (x, y, 0) T ) = ∇c where c = 2 fr (x 2 + y 2 )/2, and hence be included into the pressure-term on the right-hand side of (21a).) The evolution of a counter-and a co-rotating vortex-pair under the influence of viscosity has been studied numerically by Sipp et al. [32] and by Le Dizès and Verga [33], respectively. Recently, this process has also been investigated for the case of a Gaussian helical vortex with different core sizes and reduced pitches [16]. As shown in these studies, the evolution of the vortex until merging may be divided into two different stages, namely an initial non-viscous adaptation followed by a viscous diffusion of the vortex. During the initial stage, every vortex undergoes a rapid adaptation to the strain field produced by the curvature and the neighboring vortex filaments. This is a non-viscous process characterized by strong oscillations. Once the oscillations have decayed, the vortex settles into a slowly varying mean state where it evolves on a viscous timescale. During this stage the vortex cores expand slowly until the merging stage is reached and neighboring vortex loops rapidly approach each other [34].
To study the development of the present helical vortex, the aspect ratio of the vortex core in the P(0, 1)plane will be monitored. In this plane, the vortex centroid (η,β) T and the vortex dispersion a about the centroid are obtained fromη whereω ξ denotes the t 1 -component of the vorticity vector, the circulation is Γ = Sω ξ dηdβ, and the integrals are taken over S , which is a part of the P(0, 1)-plane that contains the vortex core [3]. The vortex dispersion (23c) expresses the spread of the vorticity and may be interpreted as a measure of the vortex core radius. Similarly, one can interpret as the vortex dispersion in the different coordinate directions and define the aspect ratio of the core as e = a β /a η [32]. The temporal evolution of the aspect ratio in the translating frame of reference is shown in Fig. 4a, in which both of the above-mentioned stages can be identified. (The evolution is similar in the rotating frame of reference.) Since the initial profile is relatively smooth, the oscillations damp out quickly such that the adaptation stage only extends up till t ≈ 0.4. From this point onwards, the solution will be close to a steady solution to the Euler equations [32,33]. Figure 4 also shows the time evolution of the core velocity U c and the vortex core radius a. The former is estimated by interpolating the velocities in S to a polar grid centered at (η,β) and averaging in the azimuthal direction, whereas the latter is estimated from (23c). As seen in Fig. 4c, the expansion of the vortex core closely follows the viscous diffusion law of two-dimensional vortices given by where a 0 is the initial core radius given in Table 1. In Fig. 5, the vorticity in the cross section at the beginning and at the end of the adaptation stage is compared. During the adaptation, the vortex core develops from an initially axisymmetric shape into an ellipse. This has been shown to be a stable shape for an inviscid Rankine vortex in an irrotational strain field, given that the strain to vorticity ratio is small [35].

Stability analysis
In this section, the stability of the generated baseflow will be studied. The stability of helical filaments has previously (among others) been considered analytically by Widnall [36] for a single helix, and by Gupta and Loewy [37] and Okulov [10] for a helix multiple.
We consider a perturbation with the ansatz function (u(x, t), p(x, t)) = (û(x),p(x))e −iωt + c.c., where ω ∈ C is an eigenvalue andû ∈ C 3 ,p ∈ C are velocity and pressure eigenvectors, respectively. This perturbation ansatz substituted into the linearized Navier-Stokes equations, leads to an eigenvalue problem governing the temporal stability of the flow. Additional information regarding the receptivity of different eigenvectors can be obtained by considering eigensolutions of the adjoint Navier-Stokes equations [38] − which are of the form where (û † ,p † ) are adjoint velocity and pressure eigenvectors. (Superscript asterisk denote complex conjugation.) These eigenvalue problems are solved using a matrix-free approach [39] together with P_ARPACK [40,41].
In similarity to (21a), the force term f on the right-hand side of (26a) corresponds to a volume force that must be included into the governing equation when considering the rotating frame of reference. However, in contrast to the baseflow calculation, there is no centrifugal force present in the perturbation formulation and in this case f = f cor = 2 fr (u y , −u x , 0) T (and correspondingly f † = 2 fr (−u †y , u †x , 0) T ).
We emphasize that the baseflow U does not correspond to a stationary solution of (21) as is commonly assumed. However, as was noted in the previous section, it is indeed an approximate solution to the steady Fig. 6 Eigenvalue spectrum. Translating frame of reference: long-wavelength instabilities (red plusses), short-wavelength instabilities (red crosses). Rotating frame of reference: long-wavelength instabilities (stars), short-wavelength instabilities (bullets). Growth rates in the rotating frame of reference together with the frequencies in the translating frame of reference estimated using (28) are plotted with open circles. The numbers refer to the azimuthal wavenumbers k in the φ-direction (invariant with respect to reference frame). A close-up on the eigenvalues corresponding to the short-wavelength instabilities in the translating frame of reference are shown in the inset plot Euler equations and can be considered to evolve on a larger timescale than the instabilities. We therefore extract the solution after the relaxation process in Sect. 2.4 is completed (i.e. at t = 0.4) and treat it as being frozen in time. Such a procedure has also been employed in previous studies e.g. [42][43][44][45]. It should be noted that since no stabilization or restriction (aside from streamwise periodicity) is imposed on the flow, the baseflow may contain a low amplitude component of possibly unstable eigenmodes at the end of the adaptation stage. However, given that the relaxation stage is short and the level of the numerical noise is low, the deformation of the baseflow due to any such modes will be negligible.
The 43 most unstable eigenvalues are sought, and in Fig. 6 the eigenvalue spectra in the rotating and the translating frame of reference are plotted. (Since the eigenvalues are complex conjugate only the positive frequencies are shown.) Comparing the eigenvalues in the two reference frames, the growth rates are seen to be the same, whereas the frequencies are related through the expression where k is the azimuthal wavenumber in the φ-direction annotated in Fig. 6. The subscripts r and t are here used to denote the rotating and the translating reference frames, respectively. See "Appendix B" for a derivation of (28). The flow is unstable to both long-and short-wavelength instabilities, which correspond to deformations of the vortex filaments and of the vortex cores, respectively. As seen, the short-wavelength instabilities are weaker than the long-wavelength instabilities, which have up to twice as large growth rate. The differences in frequencies between the two types of instabilities are several orders of magnitude. They are further discussed in the following sections.
To ensure that all the relevant flow structures are resolved, the eigenvalue computations are performed using three different polynomial orders, N = {11, 14, 17}. By comparing the different spectra, N = 14 is found to be sufficient (see the second quadrant of the mesh in Fig. 2a for a visualization of the mesh density). In addition to the resolution study, the eigenvalue computations are repeated using different temporal separation between the Krylov vectors [41] in order to certify that no aliasing is present in the spectra [39].

Mutual inductance instabilities
The obtained long-wavelength instabilities correspond to the mutual inductance instability [36], which arise from the self-and mutually induced velocities of the perturbed vortex filaments. By developing a single helical vortex in the (Rφ)z-plane and locally approximating the filaments as a single row of co-rotating infinitely long slanted vortices, an analogy between the vortex pairing due to the mutual inductance instability and an array of point vortices was established by Quaranta et al. [20]. This analogy was furthered in [17], where the growth rates for a two-helix configuration computed using the Biot-Savart law (1) following [37] were successfully compared against the growth rates of three-dimensional perturbations in a single row of infinitely long vortices derived by Robinson and Saffman [46].
Following [17], we here compare the numerical growth rates of the mutual inductance instability with the predictions by Robinson and Saffman [46]. (Details on adapting the theoretical results to the helix geometry are omitted here and the reader is referred to [17].) To enable the comparison, the numerical growth rates in Fig. 6 have been scaled with 2h 2 / , whereby the maximum growth rate becomes close to π/2 (i.e. the maximum growth rate for point vortices). As seen in Fig. 7, there is a close agreement between the numerical and the theoretical data sets. Note, however, that the present numerical setup only enables integer azimuthal wavenumbers in the φ-direction to be recovered due to the periodic boundary conditions in the streamwise direction. If the length of the domain would be extended such that more periods of the helix are taken into account, the eigenvalue spectrum would become richer.
As discussed in [20], the largest growth rate in the case of an array of point vortices is obtained for perturbations that are out-of-phase on neighboring vortices. Considering three-dimensional disturbances on two helical vortices, Fig. 7 shows that successive peaks in the growth rate curve correspond to perturbations that are either in-or out-of-phase (i.e. have phase difference equal to 0 or π radians, respectively). A sample of eigenvectors corresponding to the three most unstable modes is visualized in Fig. 8. The most unstable mode is shown in Fig. 8a. It corresponds to the eigenvalue ω = 4.742i with wavenumber k = 0. Since it has a zero frequency, the eigenvalue is not visible in Fig. 6 (due to the logarithmic scaling of the horizontal axis). Moreover, it is invariant with respect to the frame of reference. This mode corresponds to a relative displacement of the two helices relative to each other, whereby one helix contracts and the other expands. This induces a uniform pairing of the helices. The effect of the core size, the Reynolds number and the reduced pitch on the growth rate of this mode has recently been investigated in [44]. In Fig. 8b, c the eigenvectors corresponding to k = 1 and k = 2 are visualized. These modes induce a local pairing of the neighboring vortex loops at 2k locations in the azimuthal direction, and correspond to perturbations that are in-and out-of-phase on the two helices, respectively.
By visualizing the axial vorticity of the eigenvectors in Fig. 8 in a coordinate system aligned with the vortex filament (as done in Fig. 9a for the mode with k = 0), they are seen to consist of two lobes with opposite sign vorticity that is concentrated within the vortex core (r < a) and decrease in strength outside of it. This is in contrast to the shape of the corresponding adjoint eigenvector (see Fig. 9b), which vorticity field features opposite signed slender elongated structures. These structures are located far away from the cores and are seen to bridge neighboring vortex filaments. Since all the eigenvectors are normalized to have unit L 2 -norm, the resulting amplitude of the direct eigenvector can not be inferred from the values of the adjoint eigenvector [38]. However, by comparing the relative strength of the different velocity components, it is noted that inplane momentum forcing in the n j -and b j -directions tends to most effectively excite the direct eigenvectors. Fig. 7 Comparison of the growth rates (2h 2 / )Im(ω) of the mutual inductance instabilities (stars) with the theoretical growth rates of three-dimensional perturbations on infinitely long straight vortices adapted to the helix geometry for in-phase (solid line) and out-of-phase (dashed line) perturbations [17,46]. As reference, the value of π/2 is plotted (dotted line) As reference, a circle with radius a| t=0.4 (see Fig. 4c) that represents the vortex core is plotted These results are consistent with the prevailing understanding of the mutual inductance instability, since such perturbations would seek to displace the filaments orthogonally to the tangent of the vortex center lines, and thus promote pairing of the filaments.

Elliptic instabilities
The short-wavelength instabilities observed in the stability analysis correspond to the elliptic instability, which has been a topic of intense research over the past decades (see [34] for a recent review, and [47] for details regarding the curvature instability, which is another short-wavelength instability in curved Batchelor vortices).
The elliptic instability originates from a resonance between two Kelvin modes on an originally axisymmetric vortex core that has been distorted by an external strain field [48]. In a local cylindrical coordinate system aligned with the vortex filament, these modes may be described as where α, m ∈ R are the axial and the azimuthal wavenumbers, ω ∈ C is the frequency and the growth rate, anď u ∈ C 3 represents the eigenvector. The elliptic deformation of the vortex core (see Fig. 5b) due to the strain induced by the neighboring vortices can be viewed as a stationary (Re(ω) = 0) and axially homogeneous (α = 0) wave with azimuthal wavenumber m = ±2. Therefore, the condition for resonance reads where the subscripts 1 and 2 refer to the two interacting Kelvin modes [34]. Possible resonance combinations of frequencies and wavenumbers are determined by the dispersion relation, where the most unstable configurations typically correspond to crossing points between branches of the dispersion relation that have the same label l, so called principal modes [49]. The label l represents the number of zero-crossing in the radial direction of u [50]. For straight Lamb-Oseen vortices, the most unstable mode corresponds to a combination of Kelvin modes with azimuthal wavenumbers m = 1 and m = −1. However, in the case of Batchelor vortices, this configuration has been shown to stabilize with increased axial flow, and be replaced by other unstable resonance configurations [42,43]. Recently, the elliptic instability of curved Batchelor vortices has been investigated by Blanco-Rodríguez and Le Dizès [51]. These authors computed the growth rates for certain combinations of Kelvin modes using the following equation which can be applied to different vortex configurations after appropriately specifying the case specific strain rate parameter S. Blanco-Rodríguez et al. [12] gave an expression for this parameter in the case of M helical vortices as S = − 3 16 ln 2 + 1 3/2 Mε − 4 2 + 1 3 M 2 + 2 20 4 + 12 2 + 9 96 2 where and I j (·) and K j (·) are modified Bessel functions of the first and second kind. One also has that ε = Ra/(R 2 + L 2 ), and as before, = L/R is the dimensionless reduced pitch. Approximate relations for the other coefficients appearing in (31) are given in Appendix C of [51] (note also the corrigendum to some of these relations appended to [47]). The Reynolds number used in (31) is defined as Re = Re/(2π).
Since the elliptically deformed vortex cores are incorporated into the baseflow, the unstable resonance conditions discussed above will appear as unstable eigenmodes in the present stability calculations. A close-up Fig. 10 Elliptic instability modes. a Single eigenvector with azimuthal wavenumber k = 45, and b a superposition of two eigenvectors with k = 46 that correspond to the same eigenvalue. Streamwise componentû z is visualized with white and black contours indicating positive and negative velocities, respectively. The baseflow is shown in green using the λ 2 -criterion [23] (a) (b) Fig. 11 Cross section of the mode with azimuthal wavenumber k = 45 in the P(0, 1)-plane showing contours of aω ξ and b ω †ξ . As reference, a circle with radius a| t=0.4 that represents the vortex core is plotted of the numerically determined eigenvalues associated with the elliptic instability is shown in the inset of Fig. 6. In Fig. 10a, the eigenvector corresponding to the most unstable elliptic instability with the global wavenumber k = 45 is visualized. This mode assumes the shape of a pair of secondary helices that winds around the vortex filament.
In order to be able to compare the numerically obtained growth rates with the theoretical results from (31), the axial velocity and core size at the end of the relaxation process must be used. From Fig. 4b, c U c | t=0.4 = −1.5249 and a| t=0.4 = 0.0347 are obtained, which give the scaled vortex core velocity 2π (aU c )| t=0.4 = −0.3328. For a vortex core velocity of this magnitude, the resonance configuration that is expected to be most unstable is (m 1 , m 2 , l) = (−2, 0, 1) [51]. The axial vorticity inside the vortex coreω ξ for the eigenvector visualized in Fig. 10a is shown in Fig. 11a, which indeed closely resembles the structure of the mode (−2, 0, 1) (c.f. figure 3c of [51] as well as figure 15a and 16 of [42]). All other modes that appear in the spectrum of Fig. 6 and correspond to the elliptic instability are of the same type. In Fig. 12, the growth rates of these eigenmodes are plotted together with the largest root of (31) for different axial wavenumbers a| t=0.4 α (where α = k/γ ). This figure shows that the flow is unstable around the wavenumber for ideal resonance [52] α c ≈ 1.49 in the interval 1.286 ≤ a| t=0.4 α ≤ 1.788. The agreement between the numerics and the theory is overall good, although the theory slightly underestimates the growth rates of the instabilities. This small discrepancy might be associated with the approximate expressions used to evaluate the coefficients in (31), possible changes in the helix radius R (which is assumed constant during the relaxation time), the estimation of the vortex core size a, or the estimation of the axial center line velocity U c . The theoretical results are noted to be sensitive with respect to these parameters. For instance, a minor increase in the estimated core radius about 2% yields a nearly perfect match.
In Fig. 11b, the adjoint of the (−2, 0, 1)-mode is visualized. In a P(s, j)-plane, this mode assumes the shape of a spiral structure located outside the vortex core. To gain further insight into the receptivity mechanisms Fig. 12 Comparison of the numerically determined growth rates for the elliptic instability scaled with 2π a| t=0.4 (stars) and the theoretical growth rates obtained from (31) (solid line) of this mode, the components of the strain rate tensor in a Frenet-Serret frame may be examined. In such a local coordinate system, the strain rate tensor is expressed as where (∇U) is the velocity gradient tensor of the baseflow in the Cartesian coordinate system, and T(s, ψ j ) is the transformation operator between the Cartesian and the Frenet-Serret coordinate systems (see "Appendix A"). The components along the diagonal of the strain rate tensor represent the strain in the t j -, n j -and b j -directions, whereas the off-diagonal elements represent shear. These components are visualized in Fig. 13. (Note thatS(s, ψ j ) is symmetric and hence only has six unique components.) From this figure, it is readily seen that the componentsS ηη ,S ηβ andS ββ spatially coincide with the adjoint (−2, 0, 1)-mode, whereas the other components either attain their largest amplitude inside the core or are relatively week. This suggests that velocity perturbations outside the vortex core that alters the strain and shear of the vortex in a plane orthogonal to the helix line, most effectively excite the elliptic instability of the present type. Another interesting remark is that every complex eigenvalue −iω in Fig. 6 that represent an elliptic instability has a geometric multiplicity of two, meaning that there are two linearly independent eigenvectors corresponding to every such eigenvalue. The two eigenvectors either differ by (i), a simple reversal in sign of the velocity on one of the helices, or (ii), a reversal in sign of the velocity on one of the helices for the complex conjugate eigenvector phase shifted by π/2 radians. If we consider a superposition of two eigenmodes and denote the phase difference between them by χ . Then, as visualized in Fig. 10b, eigenvectors with k = {41, 43, 46, 47} that satisfy property (i) will cancel out on one helix and enhance each other on the other one for χ = 0, and vice versa when χ = π. Correspondingly, eigenvectors with k = {42, 44, 45, 48} that satisfy property (ii) will enhance each other on both helices when χ = ±π/2. Such constructive and destructive interference between the modes suggests that the breakdown to turbulence due to the elliptic instability is likely to appear different on the two vortices.

Discussion
An accurate baseflow is the starting point for every stability analysis. In most cases, the basic state will be completely determined by the geometry of the problem and the boundary conditions, and can in many situations be obtained by merely integrating the governing equations in time. Simple as it might seem, this task turns out not to be entirely trivial for a flow that features an infinitely long vortex multiple. In this situation, the obvious boundary conditions would be to impose a constant streamwise velocity or a solid body rotation in the far-field, and have periodic boundary conditions in the streamwise direction. This implies that the vortex necessarily must be positioned in the flow through the initial condition. The question thus becomes how to determine an initial field that will enable an accurate representation of the flow inside the vortex core, and still capture the external strain field due to the surrounding vortex filaments. In this respect, all the techniques that are known to us appear to fall short in one way or the other, and therefore disqualify for a study of the instability mechanisms in the flow.
In this article, a straightforward step-by-step algorithm to generate a suitable baseflow is proposed. The method has been implemented in a general purpose Navier-Stokes solver and is discussed in the context of a pair of helical vortices whose characteristics recently have been determined in an experiment [17]. The resulting 13 Strain rate tensorS evaluated for the baseflow and visualized in the P(0, 1)-plane. The superscripts ξ , η and β that appear in the subcaptions refer to the directions given by the vectors t 1 , n 1 and b 1 , respectively. As reference, a circle with radius a| t=0.4 that represents the vortex core is plotted flow field is seen to satisfy all the specified properties, and a stability analysis is observed to capture both longwavelength mutual inductance instabilities and short-wavelength elliptic instabilities. Both the translating and the rotating frame of reference have been discussed, and the stability properties determined in these different frames been compared. Although both the long-and the short-wavelength instabilities have been analyzed before, they are for the first time obtained through a single computation, which is remarkable given that the difference in frequency is several orders of magnitude. To the best of the authors' knowledge, the elliptic instability is for the first time captured and analyzed numerically in a helical vortex multiple. Subsequently, the first validation of the theoretical growth rate expressions derived by Blanco-Rodríguez and Le Dizès [51] is presented. Good agreement between the numerical and the theoretical results is reported. The geometric multiplicity of the corresponding eigenvalues suggests that this instability will develop differently in the different vortices of the helix multiple. Regarding the mutual inductance instabilities, good agreement is also observed between the numerical eigenvalues and the theoretical predictions by Robinson and Saffman [46] adapted to the helix geometry [17]. The fact that both of these instabilities are captured simultaneously, enables their relative importance to be studied. It is shown that the mutual inductance instability for the present configuration has twice as large growth rate as the elliptic instability, and hence can be expected to dominate the transition process. This finding agrees well with the experimental observations by Quaranta et al. [17], where the mutual inductance instability always was observed for the parameters listed in Table 1. In addition, the receptivity of the two instability types have been studied and adjoint eigenvectors been reported. It is shown that the mutual inductance instability is most sensitive to disturbances that enhance vortex pairing by altering the distance between neighboring filaments. By comparing the spatial support of the adjoint eigenvectors for the elliptic instability with the distribution of the individual components of the strain rate tensor for the baseflow, it is found that the (− 2, 0, 1)-mode obtained for the present flow configuration is most sensitive to perturbations outside of the vortex core that affects the straining and shearing of the vortex in a plane orthogonal to the vortex center line.
Taken together, the presented way of imposing a vortex inside the numerical domain appears rather appealing and capable of increasing our understanding about instabilities in vortex systems. Since only periodicity in the streamwise direction is assumed, the present approach can readily be used to investigate many other configurations of higher complexity. Examples involve curved helices that are inscribed in a torus, which may be used to investigate the instabilities and breakdown of horizontal axis wind turbine wakes during yawed inflow conditions [53]. Another interesting configuration that might be considered is that of twin helices arising in rotor applications where vane tips are implemented [54]. A direct means of studying the elliptic instability in a geometry with curvature and torsion is also provided. For studies where an inflow/outflow arrangement is more suitable than a periodic configuration, the derived flow field may be used to extract an accurate inflow condition. This would in turn enable the optimal forcing of the flow to be computed, as well as transient analysis involving impulse response and optimal perturbation to be carried out.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
The operators mapping vector quantities between these systems read Appendix B: Disturbances in the rotating and the translating frame of reference In this section, a relation between the frequencies of a disturbance in the rotating and the translating frame of reference is derived. We consider a rotating reference frame that rotates with the angular speed fr around a fixed axis z f , and a translating reference frame that translates along z f with the speed U fr . These coordinate systems are thus related as ⎛ ⎝ x r y r z r ⎞ ⎠ = ⎡ ⎣ cos( fr t) − sin( fr t) 0 sin( fr t) cos( fr t) 0 0 0 1 where the subscripts r and t refer to the rotating and the translating frame, respectively. By considering points along a helix curve C j (s) (that is stationary in the both frames), and inserting (7) into (38), one obtains upon simplification that φ(s r ) = φ(s t ) + fr t, φ(s r ) = φ(s t ) + U fr L t. (39a,b) Equating (39a) and (39b) yields the relation fr = U fr /L between the motions of the two reference frames, which were used in Sect. 2.2. In the local helix-oriented cylindrical coordinate system, the perturbations are assumed to be homogeneous in the θ -and s-direction and inhomogeneous in the radial direction, as specified by (29),ũ (r, θ, s, t) =ǔ(r )e i(αs+mθ −ωt) + c.c.
The flow considered in this work is helically symmetric, which means that the structures are invariant under a rotation δ/L around the z f -axis, and a translation δ along z f . The mapping between the reference frames defined by (38) corresponds to such a pair of operations. This can be realized by letting δ = U fr T denote the displacement during the time T , for which (39) yields fr T = δ/L. As a result, it follows that θ r = θ t and r r = r t . Substitution of these relations, (39a) and (8) into (29), gives the following expressions for the perturbation in the rotating and the translating frame of referencẽ u r (r r (r t ), θ r (θ t ), s r (s t ), t) =ǔ r (r t )e i(α r s t +m r θ t −[ω r r −α r γ fr] t)+ω i r t + c.c., (40a) where the superscripts r and i denote the real and imaginary parts, respectively. In order for the helical symmetry to hold, we require the spatial dependency of the perturbation to be invariant with respect to reference frame. Hence, α r = α t ≡ α, m r = m t ≡ m andǔ t =ǔ r ≡ǔ. Furthermore, we require the temporal dependency of the perturbation to be the same in the two frames, i.e. ω i t = ω i r and ω r t = ω r r − αγ fr .
Due to the assumption of periodicity in the s-direction, it follows from (8) that the perturbation also will be periodic in the global azimuthal φ-direction with the wavenumber k = αγ . It is noted from (41) that positive frequencies in the rotating frame of reference will be mapped to negative frequencies in the translating frame of reference if k is sufficiently large. This is for instance the case for the elliptic instabilities discussed in Sect. 3.2. As a result, the frequency of the mode in the rotating frame of reference increases monotonously with the wavenumber k, whereas in the translating frame of reference it does not (see Fig. 6). Since the Navier-Stokes operator governing the stability problem is real, the sets of eigenvalues {−iω r } and {−iω t } will appear in complex conjugate pairs. Therefore, we write (41) as ω r t = ω r r − sgn(ω r r )k fr , which is equivalent to (28) and is the relation between the frequencies in the rotating and the translating frame of reference plotted in Fig. 6.