1. Introduction
Vortices are present in nature as well as in many engineering applications. They are systematically created by lifting surfaces, such as wings. They take the form of trailing vortices for airplanes but exhibit a helical shape when they are generated by a rotor as for helicopters (Leishman Reference Leishman2006) or wind turbines (Hansen Reference Hansen2015). As they could significantly contribute to the velocity felt by their generating surface, it is important to know their structure and dynamics. In the present work, we consider the rotor wake model that was first introduced by Joukowski (Reference Joukowski1929). In this so-called Joukowski model, the wake is modelled by a set of vortices created by each blade: one bound vortex fixed on the blade and two free vortices emitted from the tip and root of the blade (see figure 1). Our objective is to provide the structure of the wake and its stability for all wind turbine regimes and all vertical flight regimes of a helicopter.
Joukowski vortex model (from Joukowski Reference Joukowski1929).
In the wind turbine community, most researchers have been using the blade element momentum theory for the design of wind turbines. This theory does not require a precise description of the flow. It is based on a balance between, on the one hand, thrust and torque applied to the rotor and, on the other hand, the change of axial and angular momentum experienced by the flow. It has been progressively improved to consider the presence of multiple blades, three-dimensional effects and a prescribed vortical wake (see Sørensen Reference Sørensen2016). Several wake models have been considered, either in the form of vortex sheet (Betz Reference Betz1920; Goldstein Reference Goldstein1929; Chattot Reference Chattot2003; Branlard & Gaunaa Reference Branlard and Gaunaa2015) or vortex filaments (Okulov & Sørensen Reference Okulov and Sørensen2010; Wood, Okulov & Bhattacharjee Reference Wood, Okulov and Bhattacharjee2017) but always with a prescribed cylindrical or helical symmetry. For this reason, these models cannot capture the contraction of helicopter wakes and the expansion of wind turbine wakes. In our model, as no helical symmetry is prescribed, the spatial evolution of the wake is fully described. We shall see that strongly deformed wake structures, evolving on either side of the rotor can in particular be obtained.
Our framework is the same as in Durán Venegas & Le Dizès (Reference Durán Venegas and Le Dizès2019). The vortices are concentrated in thin filaments of small but finite core size. This allows us to use a cutoff approach to compute the self-induced velocity of the vortices from Biot–Savart law (see Saffman Reference Saffman1992, for details). The structure and dynamics of the vortices are then computed using a free-vortex method, each vortex being discretized in straight-line segments advected by the flow (Leishman Reference Leishman2006). To obtain the wake structure, we first focus on steady solutions. More precisely, we look for solutions in a uniform normal incident flow that are stationary in the rotor frame. The solutions are also assumed to keep a 
$2{\rm \pi} /N$ azimuthal symmetry for an 
$N$-blade rotor, and become in the far-field uniform interlaced helices.
Perfect helical structures have been the subject of an enormous amount of work that goes back to Kelvin (Reference Kelvin1880) and Da Rios (Reference Da Rios1916). They have the particularity to rotate and translate without changing shape. Recent theoretical progress has mainly been made thanks to the expressions obtained by Hardin (Reference Hardin1982). However, there are still some debates about helix self-induced motion (Velasco Fuentes Reference Velasco Fuentes2018; Okulov & Sørensen Reference Okulov and Sørensen2020). Our objective is here to construct vortex solutions that take into account the presence of the rotor but match uniform helices in the far wake.
The stability of the solutions is also addressed. So far, most works have been concerned with the stability of uniform helices. Widnall (Reference Widnall1972) and Gupta & Loewy (Reference Gupta and Loewy1974) considered the full induction deduced from the Biot–Savart law and obtained the long-wavelength stability characteristics of a single helix and multiple helices, respectively. The effect of an axial flow in the vortex cores was analysed by Fukumoto & Miyazaki (Reference Fukumoto and Miyazaki1991). More recently, global temporal stability characteristics were also derived by Okulov (Reference Okulov2004) and Okulov & Sørensen (Reference Okulov and Sørensen2007) using Hardin's expressions. In all these works, the core size is assumed small and short-wavelength instabilities (Blanco-Rodríguez & Le Dizès Reference Blanco-Rodríguez and Le Dizès2016, Reference Blanco-Rodríguez and Le Dizès2017; Hattori, Blanco-Rodríguez & Le Dizès Reference Hattori, Blanco-Rodríguez and Le Dizès2019) that may develop in vortex cores are neglected.
A number of authors have also considered the stability problem using numerical simulations. Ivanell et al. (Reference Ivanell, Mikkelsen, Sørensen and Henningson2010) and Sarmast et al. (Reference Sarmast, Dadfar, Mikkelsen, Schlatter, Ivanell, Sørensen and Henningson2014) studied the instability of the wind turbine wake using the actuator line method (Sørensen & Shen Reference Sørensen and Shen2002) to model the rotor. They were able to calculate the spatial growth rate of individual modes through a spectral analysis of the nonlinear flow. A temporal linear stability analysis was recently performed by Selçuk, Delbende & Rossi (Reference Selçuk, Delbende and Rossi2018) and Brynjell-Rahkola & Henningson (Reference Brynjell-Rahkola and Henningson2020), obtaining the growth rate of different unstable modes using time-stepping methods on the linearized Navier–Stokes equations.
Experiments have also been performed in the wind turbine regime (Leweke et al. Reference Leweke, Quaranta, Bolnot, Blanco-Rodríguez and Le Dizès2014). By perturbing the flow of a single blade rotor, Quaranta, Bolnot & Leweke (Reference Quaranta, Bolnot and Leweke2015) were able to measure the growth rates of the main unstable modes and demonstrated a good agreement of their measurements with the theoretical predictions of Widnall (Reference Widnall1972). They also provided a convincing explanation of the instability mechanism in term of local pairing. In Quaranta et al. (Reference Quaranta, Brynjell-Rahkola, Leweke and Henningson2018) a similar analysis was performed with a two-blade rotor, and the unstable pairing modes predicted by Gupta & Loewy (Reference Gupta and Loewy1974) was captured.
In our work, we shall also consider the spatial development of the perturbations, and analyse how the stability changes from convective to absolute (Huerre & Monkewitz Reference Huerre and Monkewitz1990). A similar study has been performed for an array of vortex rings in Bolnot, Le Dizès & Leweke (Reference Bolnot, Le Dizès and Leweke2014).
The paper is organized as follows. In § 2, we introduce the framework of the vortex method applied to Joukowski's model. We provide the numerical procedures to obtain the base flow and analyse its stability. The parameters describing our problem are also provided. In § 3, the base flow is characterized as function of the parameters for a fixed blade number and a fixed vortex core size. We show that the geometry of the vortex structure strongly varies close to the rotor as the direction of the incident velocity changes. The characteristics of the induced velocity field are analysed in the rotor plane and in the far wake. Thrust and power coefficients are also computed and compared to momentum theory predictions. In § 4, the linear stability properties of the solutions are considered. The spatio-temporal development of a localized perturbation is analysed for four characteristic regimes. We first show that the theoretical temporal growth rate curves of Gupta & Loewy (Reference Gupta and Loewy1974) are recovered in the far wake for all cases. We then analyse the change of nature from convective to absolute of the instability. We show that in the far field, the spatial growth is reasonably well described by a two-dimensional (2-D) point vortex model. In the near wake, a more complex dynamics is, however, observed, especially when the vortex structure is present on both sides of the rotor plane.
2. Framework
The Joukowski rotor model is a particular model for the wake generated by a rotor. It assumes a bound vortex on the blade and two free vortices of opposite circulation emitted from the root and tip of each blade, as illustrated in the sketch by Joukowski (Reference Joukowski1929) shown in figure 1.
This model is adequate in the near wake of the rotor if the circulation profile on the blade is uniform along the span. We are rarely in these conditions, so Joukowski model should be considered instead as an approximate model for the intermediate wake where the roll-up processes of the vorticity sheet shed from the blade trailing edge have already formed two main vortices close to the tip and root for each blade. The hub vortex of the Joukowski model then results from the merging of the root vortices of each blade. By symmetry, this vortex can be positioned on the rotor axis as soon as there are at least two blades. Its circulation is equal to 
$-N\varGamma$ if there is a tip vortex of circulation 
$\varGamma$ associated with each of the 
$N$ blades. The procedure used to construct a Joukowski wake model for any circulation profile has been described in Durán Venegas, Le Dizès & Eloy (Reference Durán Venegas, Le Dizès and Eloy2019). Note, however, that if the rotor has a single blade, the root vortex cannot be positioned on the rotor axis and is expected to become helically deformed. More complex solutions are expected in this case. The far wake associated with these solutions has been analysed in Durán Venegas & Le Dizès (Reference Durán Venegas and Le Dizès2019).
To describe the vortex system, we use a free-vortex method. Vortices are assumed to be thin filaments of small core size that move in the fluid as material lines advected by the velocity field
\begin{equation} \frac{\mathrm{d} \boldsymbol{\xi}}{\mathrm{d} t} = \boldsymbol{U}( \boldsymbol{\xi}) = \boldsymbol{U}_{\infty}+ \boldsymbol{U}^{ind}( \boldsymbol{\xi}) , \end{equation}
where 
$\boldsymbol {\xi }$ is the position vector of the vortex filament, 
$\boldsymbol {U}$ the velocity field,composed of the external velocity field 
$\boldsymbol {V}_{\infty } = V_{\infty } \boldsymbol {e_z}$ and the induced velocity field 
$\boldsymbol {U}^{ind}( \boldsymbol {\xi })$ generated by the vortex filaments.
Because vortices are continuously created, it is necessary to introduce a second variable that characterizes the age of the vortices. As explained by Leishman et al. (Reference Leishman, Bhagwat and Bagai2002), a convenient way to describe the vortex system of a rotor is actually to parametrize the vortex position by two angular coordinates 
$\psi$ and 
$\zeta$ measuring the angular positions of the blade at time 
$t$, and at the time when the vortex element was created, respectively (see figure 2). As we shall see, this second angle can be viewed as a spatial coordinate along the vortex structure. This transforms (2.1) into
\begin{equation} \varOmega_R \left[\frac{\partial \boldsymbol{\xi} (\psi,\zeta)}{\partial \psi} + \frac{\partial \boldsymbol{\xi} (\psi,\zeta)}{\partial \zeta} \right] = \boldsymbol{U}_{\infty}+ \boldsymbol{U}^{ind}( \boldsymbol{\xi}), \end{equation}
where 
$\varOmega _R >0$ is the rotation speed of the rotor.
Local reference frame used by Leishman, Bhagwat & Bagai (Reference Leishman, Bhagwat and Bagai2002).
The induced velocity 
$\boldsymbol {U}^{ind}(\boldsymbol {\xi })$ can be decomposed into four different contributions
\begin{equation} \boldsymbol{U}^{ind}(\boldsymbol{\xi})=\boldsymbol{U}^{tip}(\boldsymbol{\xi})+\boldsymbol{U}^{hub}(\boldsymbol{\xi}) +\boldsymbol{U}^{blade}(\boldsymbol{\xi}) +\boldsymbol{U}^{FW}(\boldsymbol{\xi}), \end{equation}
where 
$\boldsymbol {U}^{tip}$ , 
$\boldsymbol {U}^{hub}$, 
$\boldsymbol {U}^{blade}$ and 
$\boldsymbol {U}^{FW}$ are the velocity fields induced by the tip vortices, the hub vortex, the bound vortex on the blade and the far wake, respectively. All these contributions are computed using the Biot–Savart law. The divergence of the Biot–Savart integral on the vortices is solved by introducing a small vortex core size 
$a$. This allows us to obtain the self-induced velocity by the so-called cutoff method (Saffman Reference Saffman1992). For simplicity, the vortex core size is assumed constant and identical for all the vortices.
Both hub and bound vortices have prescribed positions on the rotor axis and blade respectively. The difficulty is then to compute to position of tip vortices which are attached to the blade on one side and connected to a perfect helical structure in the far wake.
The first objective of the analysis is to obtain steady wake solutions in the frame rotating with the rotor. This amounts to find solutions of the following system of equations:
\begin{equation} \frac{\mathrm{d} r}{\mathrm{d} \zeta} = \displaystyle \frac{1}{\varOmega_R} U_r^{ind}, \quad \frac{\mathrm{d} \phi}{\mathrm{d} \zeta} = \displaystyle \frac{1}{\varOmega_R} (\varOmega^{ind} - \varOmega_R), \quad \frac{\mathrm{d} z}{\mathrm{d} \zeta} = \frac{1}{\varOmega_R} (U_z ^{ind}+ V_{\infty}). \end{equation}
where 
$r$, 
$\phi$ and 
$z$ are the radial, angular and axial positions of the vortex at 
$\zeta$ and 
$U_r^{ind}$, 
$\varOmega ^{ind}$ and 
$U_z^{ind}$ the radial, angular and axial components of the induced velocity. This system can be directly obtained from (2.2) by requiring solutions to be independent of 
$\psi$ in the rotating frame. It could have also been obtained by enforcing the vector tangent to the vortex trajectory described in cylindrical coordinates by 
$(r(\zeta ), \phi (\zeta ), z(\zeta ))$ to be always parallel to the velocity field 
$\boldsymbol {U}^{\infty }+ \boldsymbol {U}^{ind} - r \varOmega _R \boldsymbol {e_{\phi }}$ in the rotating frame. The vortex structure and the velocity field are thus steady but the vortex elements are not: they are moving along the steady vortex structure.
Each tip vortex structure is attached to the blade at a distance 
$R_b$ from the rotation axis. Assuming a 
$2{\rm \pi} /N$ azimuthal symmetry for an 
$N$ blade rotor, we can consider a single tip vortex structure, whose starting position is chosen to be
\begin{equation} r(0)= R_b,\quad \phi(0) = 0,\quad z(0)=0. \end{equation} As in Durán Venegas & Le Dizès (Reference Durán Venegas and Le Dizès2019) for generalized helical pairs, (2.4a–c) are solved by an iterative method using a finite difference scheme. The tip vortex structure is discretized in small elements with at least 25 points by turn. This allows us to use explicit expressions for the induced velocity with distant vortex segments contributions and a local self-induced contribution associated with local curvature (see Durán Venegas & Le Dizès Reference Durán Venegas and Le Dizès2019). A sufficiently large number of turns are considered such that the structure is almost helical at the end of the computation domain. This number depends on the operational conditions and has varied between 15 and 40 in our study. At the end of the computational domain, estimates for the helix pitch 
$h_{FW}$ and radius 
$R_{FW}$ are obtained and used to construct the far-wake structure. At least 15 turns of this far-wake structure are used to compute the induced velocity of the far wake (see figure 3).The dependence of the results with respect to the numerical parameters (number of discretization points by turn, number of turns in the computation domain, number of turns in the far wake) is analysed in appendix A. The numerical parameters have been chosen such that the results are fully converged. The iterative method is initiated by a perfect helix solution or by a solution obtained for close parameters.
Schematic of the complete vortical structure generated at the wake (to simplify the figure, only the tip vortex from a single blade is represented). The solid line shows the calculation domain and the dashed line shows the prescribed far-wake structure. A smaller number of turns are shown compared to that used for the computation.
The vortices are defined by their circulation 
$\varGamma$ and their core size 
$a$. Both are assumed constant and uniform. Tip vortices are emitted at the radial coordinates 
$R_b$. If the rotor has 
$N$ blades, rotates at the angular velocity 
$\varOmega _R$ in an incident flow of axial velocity 
$V_{\infty }$, we can then define four non-dimensional parameters
\begin{equation} \lambda=\frac{R_b \varOmega_R}{V_{\infty}}, \quad \eta=\frac{\varGamma}{R_b^{2} \varOmega_R}, \quad \varepsilon=\displaystyle\frac{a}{R_b}, \quad N. \end{equation}
The parameter 
$\lambda$ is known as the tip-speed ratio. It can be either positive or negative. By convention, we consider that 
$\lambda$ is negative for a helicopter in ascent flight (
$V_{\infty }<0$), and positive in descent flight (
$V_{\infty }>0$). The wind turbine regime is then obtained for positive values of 
$\lambda$. The parameter 
$\eta$ measures the relative strength of the vortices: it is related to the so-called blade loading. This parameter can be chosen positive if we define the circulation of the bound vortex as 
${ \boldsymbol {\varGamma }_{bound}} = \varGamma \boldsymbol {e_r}$. The parameter 
$\eta$ is in general smaller than 0.1. The parameter 
$\varepsilon$ characterizes the vortex core size with respect to the rotor size. This parameter will always be considered as small, and typically equal to 0.05 or 0.01. The number, 
$N$, of blades will be kept fixed to 2.
3. Description of the steady solutions
3.1. Geometry
For fixed 
$\varepsilon$ and 
$N$, when the tip-speed ratio 
$\lambda$ and vortex strength parameter 
$\eta$ are varied, different vortex geometries are obtained. These geometries are illustrated for a fixed 
$\eta$ and different 
$\lambda$ in figures 4 and 5.
Vortex geometry for a two-bladed rotor in wind turbine regimes; (a) 
$\lambda =3.3$, (b) 
$\lambda =4.3$, (c) 
$\lambda =5.1$, (d) 
$\lambda =6.1$. For all cases 
$\eta = 0.05$ and 
$\varepsilon = 0.01$.
Vortex geometry for a two-bladed rotor in helicopter regimes; (a) 
$\lambda =-10$, (b) 
$\lambda =\infty$, (c) 
$\lambda =19$, (d) 
$\lambda =11$. For all cases 
$\eta = 0.05$ and 
$\varepsilon = 0.01$.
Two main types of solutions are obtained according to the behaviour of the solution away from the rotor: solutions in the wind turbine regime for which the vortex structure expands and goes upwards (figure 4); solutions in the helicopter regime for which the vortex structure contracts and goes downwards (figure 5). This corresponds to the two domains of parameters, which are indicated in figure 6 for 
$\varepsilon = 0.01$ and 
$N=2$ on either side of the grey region where no solution has been obtained.
Map of the rotor regimes as a function of 
$\lambda$ and 
$\eta$. The grey region represents the region where no solution has been obtained. On the left, helicopter regimes: the wake contracts and is advected downwards. On the right, wind turbine regimes: the wake expands and is advected upwards. In both regimes, the dashed line corresponds to the boundary of the domain where the tip vortices are on either side of the rotor plane. The other boundary of that domain is the grey region.
Before describing the geometry of the vortices close to the rotor, let us first look at the characteristics of the far wake. By construction, the tip vortices form perfect helices in the far wake. The radius 
$R_{FW}$ and pitch 
$h_{FW}$ of these helices are shown in figure 7 for the parameters of figure 6. By convention, we have chosen the pitch to be positive when the vortices go upwards, and negative when they go downwards. We can see that both the wake contraction in the helicopter regime and the wake expansion in the wind turbine are reduced as the tip-speed ratio 
$|\lambda |$ decreases. An opposite tendency is by contrast observed for the far-wake helix pitch that increases with 
$|\lambda |$.
Contour maps of the far-wake helix characteristics for 
$\varepsilon = 0.01$; 
$N=2$, non-dimensionalized by 
$R_b$. (a) Radius 
$R_{FW}/R_b$; (b) pitch 
$h_{FW}/R_b$. By convention, 
$h_{FW}$ is negative (respectively positive) when the vortices are advected downwards (respectively upwards). The dashed lines represent the theoretical formulae (3.3a) (yellow) and (3.3b) (red).
An estimate for 
$h_{FW}$ can be deduced from simple arguments. As the helical structure is stationary (in the frame rotating at the angular speed 
$\varOmega _R$), its rotation must be balanced by its axial displacement, that is 
$V_{FW}/ \Omega _{FW} = - h_{FW}/(2{\rm \pi} )$ (assuming 
$\varOmega _{FW}<0$). An approximation to the axial displacement speed 
$V_{FW}$ can be obtained from the displacement speed of a double array of opposite point vortices separated by a distance 
$h_{FW}/N$ in a background axial flow 
$V_{\infty }$, that is 
$V_{FW}= -N \varGamma /(2 h_{FW}) + V_{\infty }$ in the wind turbine regimes and 
$V_z = N \varGamma /(2 h_{FW}) + V_{\infty }$ in the helicopter regimes. If we neglect the self-induced rotation of the structure, that is assume that 
$\Omega _{FW}\approx - \Omega _R$, the condition of stationarity then gives
\begin{gather} h_{FW} = \frac{2{\rm \pi}}{\varOmega_R} \left(\frac{N\varGamma}{2 h_{FW} } + V_{\infty}\right), \quad \text{for helicopter regimes}\ (h_{FW} <0), \end{gather}
\begin{gather}h_{FW} = \frac{2{\rm \pi}}{\varOmega_R} \left(-\frac{N\varGamma}{2 h_{FW} } + V_{\infty}\right) , \quad \text{for wind turbine regimes}\ (h_{FW} >0). \end{gather}
This gives a second-order equation for 
$X=h_{FW}/R_b$ of the form
\begin{equation} X^{2} - \frac{2{\rm \pi}}{\lambda} X \pm N {\rm \pi}\eta =0, \end{equation}
where the sign 
$-$ is for helicopter regimes, and the sign 
$+$ for wind turbine regimes. The adequate solution for each case is
\begin{gather} h_{FW}/R_b = \frac{\rm \pi}{\lambda} - \sqrt{\frac{{\rm \pi}^{2}}{\lambda^{2}} + N {\rm \pi}\eta} ,\quad \text{for helicopter regimes}\ (h_{FW} <0), \end{gather}
\begin{gather}h_{FW}/R_b = \frac{\rm \pi}{\lambda} + \sqrt{\frac{{\rm \pi}^{2}}{\lambda^{2}} - N {\rm \pi}\eta} ,\quad \text{for wind turbine regimes}\ (h_{FW} >0). \end{gather}
These contours are plotted with dashed lines in figure 7(b). We can see that these formulae provide a reasonably good approximation of 
$h_{FW}/R_b$, especially for small 
$\eta$.
The wind turbine regime is associated with large values of 
$1/\lambda$. It also corresponds to the so-called windmill brake regime of helicopters in rapid descent flight. In this regime, the vortices move in the direction of the external wind, and the vortex structure expands. As seen in figure 8(b), the radius of the structure increases monotonously up to its far-wake value. By contrast, the vortex pitch does not seem to change much along the structure at least as long as 
$\lambda$ is not too large (figure 9b). The most interesting feature is the existence, for each 
$\eta$, of a critical value of 
$\lambda$ above which the tip vortices are first emitted downwards before being advected upwards.This transition is materialized by a dashed line on the right of the grey region in figure 6. As 
$\lambda$ gets larger, the expansion continues to increase and finally gets too large to fit in the computation domain. This provides the limit of existence of the wind turbine solutions.
Radial position of the wake in the axial direction for different tip-speed ratio values. With 
$\eta =0.05$, 
$\varepsilon =0.01$ and 
$N=2$. (a) Helicopter regimes. (b) Wind turbine regimes.
Evolution of the local pitch along the wake for different tip-speed ratio values. With 
$\eta =0.05$ and 
$\varepsilon =0.01$. (a) Helicopter regimes. (b) Wind turbine regimes. By convention, 
$h_{FW}$ is negative (respectively positive) when the vortices are advected downwards (respectively upwards).
The change of direction of the vortex emission as 
$\lambda$ increases has been observed experimentally. In figure 10(a,c) we visualize using dye the time-averaged trajectories of the vortices in a rotor side view for values of 
$\lambda$ before and after the transition. We clearly see on the bottom view obtained for 
$\lambda =11.2$ that the trajectory of the tip vortex goes below the rotor plane. Unfortunately, vortex circulation and core size have not been measured in the experiment, so we cannot perform a quantitative comparison. Nevertheless, figure 10 demonstrates that we can obtain a good agreement for the geometry with the numerical model for similar values of 
$\lambda$.
Time-averaged dye visualization of the wake behind half of the rotor plane in a wind turbine configuration at 
$\lambda = 7.23$ (a,b) and 
$\lambda =11.2$ (c,d) (from Quaranta Reference Quaranta2017). (a,c) Qualitative comparison of the wake geometry with the numerical model for wind turbine configuration at high tip-speed ratio (up: 
$\lambda =8$, down: 
$\lambda =10$; 
$\eta = 0.02$, 
$\varepsilon =0.01$). The rotor disc is represented in black dashed line. In black solid line, the projection of the tip vortices and, in red, their radial evolution.
In the helicopter domain of parameters (left side of the grey region in figure 6), we cover the normal flight conditions of a helicopter: climbing flight for negative 
$1/\lambda$, hovering for 
$1/\lambda =0$ and weak descent flight for positive (and small) values of 
$1/\lambda$. In all these cases, the vortex structure matches away from the rotor a helical structure that moves downwards, although the external wind is downwards for negative 
$\lambda$ and upwards for positive 
$\lambda$.
The main modifications of the geometry as 
$\lambda$ varies are observed close to the rotor plane. In figures 8(a) and 9(a) we display the variation in the wake of the radius and pitch of the vortex structure, respectively. As already mentioned above, we can observe in figure 8(a), that the contraction of the vortex structure becomes more important as 
$1/\lambda$ increases. As for the wind turbine regime, we observe that there exists a critical value of 
$1/\lambda$ above which tip vortices are no longer emitted towards its far wake, but in the opposite direction. This transition is indicated by a dashed curve on the left of the hatched region in figure 6. For values of 
$1/\lambda$ larger than this critical value, the vortices go upwards before being advected downwards. Accordingly, the local pitch changes sign as we move along the structure: it is positive when the vortices move upwards close to rotor, and becomes negative afterwards when vortices start to be advected downwards (see figure 9a).
It is interesting to note that the transition curve is close to the line 
$1/\lambda =0$ but not identical to this line. For small 
$\eta$ (low loading), this transition occurs in a climbing flight configuration, whereas this transition is observed in the descent flight regime for large 
$\eta$ (strong loading). As 
$1/\lambda$ is increased above the transition, the vortex excursion above the rotor plane becomes more important. It becomes more and more difficult to follow the solution as the number of turns above the rotor plane grows as 
$1/\lambda$ increases.
For all the parameters tested, we have not found any overlap between the wind turbine domain and the helicopter domain. However, we cannot exclude this possibility if a larger computational domain is considered for wind turbine solutions or if helicopter solutions are continued further.
3.2. Induced velocities
For applications, especially the design of rotors, it is useful to know the velocity field induced by the vortex structure. In this section, we provide the induced velocity field for four different configurations typical of each identified flow regime. We fix 
$\eta$, 
$\varepsilon$ and 
$N$ (
$\eta =0.05$, 
$\varepsilon =0.01$, 
$N=2$) and consider different values of 
$\lambda$ (
$\lambda =-20$, 
$\lambda =11$, 
$\lambda =6.2$, 
$\lambda =4.3$). For convenience, we normalize the induced axial velocity by the outer velocity 
$V_{\infty }$, and the induced azimuthal velocity by the azimuthal velocity 
$\varGamma / (2{\rm \pi} r)$ of a vortex filament of circulation 
$\varGamma$.
In figure 11 are shown contours plots of each velocity component in the rotor plane for each case. In figure 12 we compare, for each case and each component, the induced velocity on the blade (black dashed line), the mean (azimuthal average) induced velocity in the rotor plane (black solid line) and the mean induced velocity in the far wake (red solid line).
Induced velocity contours in the rotor plane (
$z=0$ plane) for 
$\eta = 0.05$, 
$\varepsilon = 0.01$, 
$N=2$ and various 
$\lambda$ (
$\lambda = -20, 11, 6.2, 4.3$); (a–d) 
$V_z^{ind} /V_{\infty }$, (e–h) 
$V_{\phi }^{ind} 2{\rm \pi} r/\varGamma$.
Induced velocity profiles. (a,c,e,g) Axial velocity 
$V_z^{{ind}}/V_{\infty }$, (b,d,f,h) azimuthal velocity 
$V_{\phi }^{{ind}} 2{\rm \pi} r/\varGamma$ for various values of 
$\lambda$ and 
$\eta =0.05$, 
$\varepsilon =0.01$ and 
$N=2$. In black, the profiles in the rotor plane (dashed: on the blade, solid: azimuthal average), in red, the azimuthal average in the far wake.
The first point to emphasize is the axisymmetric character of the induced flow (figure 11). In the rotor plane, both velocity components are almost axisymmetric. In the rotor plane, the difference between the azimuthally averaged field and the local field is mainly concentrated in the close neighbourhood of the vortices (see figure 11). This is also clearly visible in figure 12, where we can see that the velocity on the blade (black dashed line) departs from its azimuthal average (black solid line) at specific locations: at the blade tip, and at the location where one vortex crosses the rotor plane (
$r/R_b \approx 0.6$ for 
$\lambda =11$; 
$r/R_n\approx 2.1$ for 
$\lambda =6.2$).
The azimuthally averaged profiles of the induced velocity possess interesting features. As soon as we are outside the vortex structure, both the axial and the azimuthal components vanish. This is true in the far field and in the rotor plane as well (see figure 12). In the far field, both components are also constant inside the vortex structure. This is in agreement with known properties of the induced velocity of infinite helices (see Hardin Reference Hardin1982). In the rotor plane, we can note that both components are also constant inside the vortex structure for 
$\lambda =-20$ and 
$\lambda =4.3$. Moreover, for these cases, the constant is approximatively equal to the half of the induced axial velocity in the far field, as assumed in the momentum theory (e.g. Sørensen Reference Sørensen2016). For the two cases where the vortices cross the rotor plane, axial and azimuthal components are only almost constant by part: there is a jump at each radial location where a vortex crosses the rotor plane. For these cases, the relation between the mean induced axial flow in the rotor plane and in the far field is less obvious. It corresponds to the situations where momentum theory does not apply.
In figure 13 we compare the azimuthal average of each component in the rotor for the different values of 
$\lambda$. This figure is interesting to analyse the effect of 
$\lambda$ on the mean profile in the rotor plane. When normalized by 
$\varGamma /(2{\rm \pi} r)$, the azimuthal profile is particularly simple: in the central part of the rotor, it is always, +1 in the helicopter regimes and 
$-1$ in the wind turbine regimes. When no vortices cross the rotor plane, the azimuthal velocity just jumps to zero at the rotor edge (
$\lambda =-20$, 
$\lambda =4.3$). In the more complex configurations, where the vortices cross the rotor plane (
$\lambda =11$, 
$\lambda =6.2$), the normalized azimuthal velocity exhibits an addition jump at the radial location of crossing. For the helicopter regimes (
$\lambda =11$), it jumps from +1 to 
$-1$ at the radial location of crossing then back to 0 at the edge of the disk. For the wind turbine regimes (
$\lambda =6.2$), it jumps from 
$-1$ in the centre to 
$-2$ at the edge of the disk then to 0 at the radial location of crossing. Although complex, the azimuthal induced velocity is therefore easily modelled in all the configurations. The axial component has also a simple structure: it is approximatively constant by part. Except for 
$\lambda =11$, the induced axial flow is found to be fairly constant in the whole rotor disk, and therefore close to the mean value in the rotor disk. However, this value strongly varies with 
$\lambda$. For 
$\lambda =11$, that is in the helicopter descent flight regimes where the vortices cross the rotor disk, the induced flow is strong between the centre and the location where the vortices cross the rotor, and just compensates the mean advection velocity in the external part of the rotor. The mean induced axial flow on the rotor disk is therefore not expected to be a relevant quantity in that case as only a reduced inner part of the rotor contributes.
Azimuthal average of the induced velocity in the rotor plane for different values of 
$\lambda$ and 
$\eta =0.05$, 
$\varepsilon =0.01$ and 
$N=2$. (a) Axial velocity 
$V_z^{{ind}}/V_{\infty }$. (b) Azimuthal velocity 
$V_{\phi }^{{ind}} 2 {\rm \pi}r/\varGamma$.
In figure 14, the streamlines of the azimuthally averaged flow are displayed in the 
$(r,z)$ plane for the four typical regimes. In solid lines are plotted the streamlines of the induced flow, while the streamlines of the total flow are represented with dashed lines. In fast climbing helicopter regimes and in wind turbine regimes, the external flow is dominant, so the total flow remains unidirectional. When the external flow is weaker and opposed to the induced velocity, the flow can become bidirectional (see figure 14b,c); the vortex tube that goes through the rotor disk starts and ends on the same side of the rotor plane. For these configurations, as already mentioned above, momentum theory does not apply.
Streamlines of the azimuthally averaged flow in the 
$(r,z)$ plane for (a) 
$\lambda =-20$, (b) 
$\lambda =19$, (c) 
$\lambda =6.1$, (d) 
$\lambda =4.3$ with 
$\eta = 0.05$, 
$\varepsilon = 0.01$ and 
$N=2$. Black solid line: induced flow. Black dashed line: total flow. Red solid line: radial position of the tip vortex. Black dotted line: rotor disc.
3.3. Thrust and power coefficients
In both the helicopter and the wind turbine communities, it is common to characterize a rotor using the thrust and the power coefficient. These coefficients are non-dimensional forms of thrust and power. They can be obtained using the local lift force 
$\textrm {d} \boldsymbol {L}$ derived from Kutta–Joukowski theorem (e.g. Sørensen Reference Sørensen2016, chap. 9.2)
\begin{equation} \textrm{d} \boldsymbol{L}= \rho ({ \boldsymbol{V}_{\infty}} + { \boldsymbol{U}^{{ind}}} ) \times { \boldsymbol{\varGamma}_{{bound}} }\,\textrm{d}r, \end{equation}
where 
${ \boldsymbol {\varGamma }_{{bound}}} = \varGamma \boldsymbol {e_r}$ is the constant circulation of the bound vortex on the blade. The thrust 
$T$ and the torque 
$Q$ are then defined by integrating on the blade the local thrust 
$\textrm {d}T= \textrm {d} \boldsymbol {L} \boldsymbol {\cdot } \boldsymbol {e_z}$ and the local torque 
$\textrm {d}Q= ( \boldsymbol {r} \times \textrm {d} \boldsymbol {L}) \boldsymbol {\cdot } \boldsymbol {e_z}$
\begin{gather} T= \rho N \varGamma \int_a^{R-a} (r \varOmega _R - U_{\theta}^{{ind}})\,\textrm{d}r , \end{gather}
\begin{gather}Q= \rho N \varGamma \int_a^{R-a} (V_{\infty} + U_z^{{ind}})r\,\textrm{d}r. \end{gather}
Note that we neglect the contribution coming from the small intervals in the cores of hub and tip vortices where the circulation of the bound vortex should in principle be modified. The power 
$P$ delivered (or received) by the rotor is then given by
\begin{equation} P= \varOmega_R Q = \rho N \varGamma \varOmega_R \int_a^{R-a} (V_{\infty} + U_z^{{ind}}) r\,\textrm{d}r . \end{equation}
In agreement with our convention, 
$T$ is always positive, but 
$P$ may change sign. We expect 
$P<0$ in the helicopter regime, but 
$P>0$ in the wind turbine regime.
For helicopters, thrust and power coefficients are usually defined by
\begin{equation} C_T^{(Hel)} = \frac{T}{\frac{1}{2}\rho {\rm \pi}R^{4}\varOmega_R^{2} },\quad C_P^{(Hel)} = \frac{P}{\dfrac{1}{2}\rho {\rm \pi}R^{5}\varOmega_R^{3} } . \end{equation}
In figure 15 we display the contour levels of these coefficients in the 
$(1/\lambda , \eta )$ plane. Also indicated by red dashed lines is the approximation obtained when the induced velocity in the rotor plane is neglected
\begin{equation} C_{T0}^{(Hel)} \approx \frac{N \eta}{\rm \pi}\;;\quad C_{P0}^{(Hel)} \approx \frac{N \eta}{{\rm \pi} \lambda}. \end{equation}
As expected from the above formula, we observe that the thrust coefficient 
$C_T^{(Hel)}$ mainly depends on 
$\eta$. Note, however, the effect of the azimuthal induced flow; it decreases the thrust in the helicopter regime and increases it in the wind turbine regime.
(a) Thrust coefficient 
$C_T^{(Hel)}$ and (b) power coefficient 
$C_P^{(Hel)}$ for a Joukowski rotor with 
$N=2$, 
$\varepsilon =0.01$. The red dashed lines are the contours of the approximation obtained without taking into account the induced flow [
$(0.01, 0.02, 0.03, 0.04)$ for 
$C_T^{(Hel)}$, 
$(-0.004,-0.002,0,0.002,0.004)$ for 
$C_P^{(Hel)}$]. The solid yellow line in (b) is the contour 
$C_P^{(Hel)} =0$.
The mean induced velocity 
$V_i$ is an important quantity. It can be calculated directly by integrating 
$U_z^{{ind}}$ in the rotor disk or by computing the ratio 
$P/T$ that gives 
$V_i +V_{\infty }$. In figure 16(a), we have plotted the ratio 
$V_i/V_h$ versus 
$V_{\infty }/V_h$, where 
$V_h$ is the mean induced velocity in hover (that is when 
$V_{\infty }=0$). The calculation is performed for a fixed 
$C_T^{(Hel)}$. The results are compared to the momentum theory predictions which apply for 
$V_{\infty }/V_h < -2$ (wind turbine regime) or 
$V_{\infty }/V_h>0$ (helicopter climbing regime) and to the experimental data fit proposed by Leishman (Reference Leishman2006) for 
$-2 < V_{\infty }/V_h < 0$. The best agreement with the momentum theory is obtained for low loading data (
$C_T^{(Hel)}=0.01$, indicated by circles). A systematic departure is nevertheless observed when the mean induced velocity is computed from the power. We suspect that it comes from an overestimation of the contributions of the blade tip to the torque. Note also that the agreement is less good on the wind turbine side for large 
$C_T^{(Hel)}$. This is clearly an effect of the azimuthal induced velocity in the rotor plane. When this contribution is neglected, the departure from the momentum theory almost disappears, as seen in figure 16(b).
Mean induced velocity 
$V_i/V_h$ versus 
$V_{\infty }/V_h$, where 
$V_h$ is the mean induced velocity in hover for a fixed 
$C_T^{(Hel)}$ (
$C_T^{(Hel)}= 0.01$ with circles, 
$C_T^{(Hel)}=0.02,0.03,0.04$ with dots). The fixed parameters are 
$N=2$, 
$\varepsilon =0.01$. The other parameters are such that 
$-0.2<1/\lambda <0.3$ and 
$0.01<\eta <0.1$. The mean induced axial velocity 
$V_i$ is obtained from 
$V_i=P/T-V_{\infty }$ (red symbols) or by taking the mean of 
$U_z^{{ind}}$ in the rotor plane (blue symbols). The solid line is the momentum theory prediction (valid for 
$V_{\infty }/V_h$ smaller than 
$-2$ or larger than 0). The thick dashed line is a fit of experimental data provided by Leishman (Reference Leishman2006, (2.96), p. 87). (a) The data are computed using formula (3.5a) for 
$T$ and (3.6) for 
$P$. (b) The data are computed neglecting 
$U_{\theta }^{{ind}}$ in (3.5a) for 
$T$ and using 
$\langle U_z^{{ind}}\rangle _{\theta }$ instead of 
$U_z^{{ind}}$ in (3.6) for 
$P$.
For wind turbines, a different normalization is used for the thrust and power coefficients. They are defined by
\begin{gather} C_T^{(W)} = \frac{T}{\dfrac{1}{2}\rho {\rm \pi}R^{2} V_{\infty}^{2} } = \lambda ^{2} C_T^{(Hel)}, \end{gather}
\begin{gather}C_P^{(W)} = \frac{P}{\dfrac{1}{2}\rho {\rm \pi}R^{2} V_{\infty}^{3}} = \lambda^{3} C_P^{(Hel)}. \end{gather}
The simplest version of the momentum theory provides formulae for these coefficients in terms of the interference factor 
$a_i=V_i/V_{\infty }$ (e.g. Sørensen Reference Sørensen2016)
\begin{equation} C_T^{(W)} = 4 a_i (1-a_i),\quad C_P^{(W)} = 4 a_i (1-a_i)^{2}. \end{equation}
Our data are tested against these expressions in figure 17. This figure shows that momentum theory is recovered for small 
$a_i$. As expected, a departure is observed for 
$a_i>0.5$. Thrust continues to increase as 
$a_i$ gets larger than 0.5 contrarily to the momentum theory prediction. The induced field impact is clearly seen on these plots. It increases the thrust but strongly decreases the power. A difference between the data obtained with the mean induced velocity and those with the induced velocity on the blade is only visible for the power coefficient. When the induced velocity on the blade is used, the power is found to be smaller and further away from the momentum theory prediction.
Wind turbine thrust coefficient (a) and wind turbine power coefficient (b) versus the interference factor 
$a_i=V_i/ V_{\infty }$. The solid lines are the momentum predictions (3.10a,b). All the data are in the wind turbine regime for 
$N=2$, 
$\varepsilon =0.01$ and 
$0.01<\eta <0.1$. The green symbols are obtained using formula (3.5a) for 
$T$ and (3.6) for 
$P$. The blue symbols are obtained using the same formulae but replacing the induced velocity on the blade by its azimuthal average. The red symbols are when the induced velocity field is neglected.
It worth mentioning that in this section the local lift force on the blade has been calculated assuming a uniform circulation profile along the span (see formula (3.4)). When the circulation is not uniform on the blade, an equivalent Joukowski wake model can still be constructed, as mentioned in § 2. But, in that case, the local lift force can be computed differently, using the exact circulation profile on the blade and the azimuthal average of the induced velocity field, as in Durán Venegas et al. (Reference Durán Venegas, Le Dizès and Eloy2019) for instance.
4. Stability
In this section, we analyse the evolution of linear perturbations to the base flow obtained in the previous section.
4.1. Perturbation model
We now consider the full problem (2.2) and search for solutions of the form
\begin{equation} \boldsymbol{\xi}(\psi,\zeta) = \boldsymbol{\xi_0}(\zeta) + \boldsymbol{\xi'}(\psi,\zeta), \end{equation}
where 
$\boldsymbol {\xi _0} = (r_0(\zeta ),\phi _0(\zeta ),z_0(\zeta ))$ is a base flow satisfying (2.4a–c), and 
$\boldsymbol {\xi '}= (r',\phi ', z')$ a small perturbation. Linearizing (2.2), we obtain the following partial differential equations for the perturbation:
\begin{gather} \frac{\partial r'}{\partial \psi} = - \frac{\partial r'}{\partial \zeta} + \frac{1}{\varOmega_R} \left[ \left.\frac{\partial U_r}{\partial r}\right|_{ \boldsymbol{\xi_0}} r' + \left.\frac{\partial U_r}{\partial \phi}\right|_{ \boldsymbol{\xi_0}} \phi' + \left.\frac{\partial U_r}{\partial z}\right|_{ \boldsymbol{\xi_0}} z' \right], \end{gather}
\begin{gather}\frac{\partial \phi '}{ \partial \psi} = - \frac{\partial \phi'}{\partial \zeta} + \frac{1}{\varOmega_R} \left[\left. \frac{\partial \varOmega}{\partial r} \right|_{ \boldsymbol{\xi_0}}r' +\left. \frac{\partial \varOmega}{\partial \phi} \right|_{ \boldsymbol{\xi_0}}\phi' + \left.\frac{\partial \varOmega}{\partial z} \right|_{ \boldsymbol{\xi_0}}z' \right], \end{gather}
\begin{gather}\frac{\partial z'}{\partial \psi} = - \frac{\partial z'}{\partial \zeta} + \frac{1}{\varOmega_R} \left[ \left.\frac{\partial U_z}{\partial r}\right|_{ \boldsymbol{\xi_0}} r' + \left.\frac{\partial U_z}{\partial \phi} \right|_{ \boldsymbol{\xi_0}}\phi'+\left. \frac{\partial U_z}{\partial z} \right|_{ \boldsymbol{\xi_0}}z' \right]. \end{gather} The system is solved numerically. For this purpose, a finite difference discretization scheme has been implemented in 
$\psi$ and 
$\zeta$
\begin{equation} \frac{ \boldsymbol{\xi'}_{j+1} - \boldsymbol{\xi'}_{j} }{\Delta \psi} = \left[ -D_{\zeta} + \frac{1}{\varOmega_R} \boldsymbol{\nabla} \boldsymbol{U}|_{ \boldsymbol{\xi_0}} \right] \boldsymbol{\xi'}_j, \end{equation}
where the index 
$j$ represents the discretized position in 
$\psi$, 
$\Delta \psi$ the discretization interval in 
$\psi$, 
$D_{\zeta }$ the discretization matrix on 
$\zeta$ and 
$\boldsymbol {\nabla } \boldsymbol {U}|_{ \boldsymbol {\xi _0}}$ the Jacobian matrix of the induced velocity field of the base flow. To avoid numerical stability problems, a robust discretization scheme has to be chosen for the matrix 
$D_{\zeta }$. In our case, a four point backward discretization scheme has been used.
The perturbation is introduced as a displacement impulse in the axial direction at the azimuthal position 
$\zeta _p$
\begin{equation} \boldsymbol{\xi'}_{j=0}(\zeta) = \begin{pmatrix} 0 \\ 0 \\ A_p \delta(\zeta-\zeta_p) \end{pmatrix}, \end{equation}
where 
$A_p$ is the amplitude of the initial perturbation (that could be fixed to 1 since the problem is linear) and 
$\delta$ is the Dirac delta function. We keep the 
$2{\rm \pi} / N$ azimuthal symmetry of the base flow, so the Dirac impulse is applied to the 
$N$ tip vortices simultaneously.
In figure 18, the evolution of the axial displacement of the vortices induced by the Dirac perturbation is shown for a helicopter climbing flight regime. Here, the initial perturbation has been placed at 
$\zeta _p = 4{\rm \pi}$. We observe that the perturbation grows, extends and propagates downstream. Moreover, it forms alternate peaks that are associated with a particular interaction with the other parts of the helices.To understand this interaction, it is useful to look at the deformation of the helices induced by the perturbation. In figure 19(a) we plot in an axial–azimuthal view the undeformed helices as well as the helices deformed by the perturbation; we clearly see that the two vortices tend to get closer at specific locations. A similar behaviour has already been observed in experiments and shown to correspond to a local pairing instability (Quaranta et al. Reference Quaranta, Bolnot and Leweke2015, Reference Quaranta, Brynjell-Rahkola, Leweke and Henningson2018). A longitudinal cut across the helices shows that the displacement of the vortices towards each other is actually in phase opposition and makes a 
$45^{\circ }$ angle with respect to line connecting both vortices (see figure 19b), as for the pairing instability of an array of point vortices (Lamb Reference Lamb1932).
Axial displacement induced by the perturbation at three different instants (
$\psi =0, 2{\rm \pi} , 4 {\rm \pi}$). Parameters of the base flow are 
$\lambda = -10$, 
$\eta = 0.04$, 
$\varepsilon = 0.05$, 
$N=2$. Perturbation amplitude and initial position 
$A_p=10^{-4}R_{tip}$ and 
$\zeta _p = 4{\rm \pi}$.
Illustration of the instability. (a) View of the tip vortices in the axial–azimuthal plane. Solid line: perturbed solution. Dashed line: unperturbed solution. (b) View in longitudinal cross-cut of the tip vortices (on one side). Top: unperturbed; bottom: perturbed.
4.2. Temporal stability properties
In this section, we first consider the properties of the far field by introducing the perturbation sufficiently far away from the rotor (typically 
$\zeta _p=20{\rm \pi}$). As the flow is approximatively homogeneous at this location, each perturbation wavenumber is thus expected to have its own independent amplitude evolution. One of the advantages of using a Dirac perturbation is that it excites all the wavenumbers simultaneously. A single simulation can therefore be used to obtain the complete stability diagram of the flow. In practice, we proceed as follows.At each time step 
$\psi$, we perform a spatial Fourier transform in 
$\zeta$ of the function 
$z'(\psi ,\zeta )$ in order to obtain the amplitude 
$z'_k(\psi )$ associated with each wavenumber 
$k$. We then estimate the growth rate 
$\sigma (k)$ of each wavenumber, by looking at the slope of 
$\log (z'_k(\psi ))$ in terms of 
$\psi$. For four typical cases corresponding to the different regimes, we obtain the curves shown in figure 20(a). By construction, the growth rate 
$\sigma$ and the wavenumber 
$k$ that we get by this method are dimensionless. The growth rate has been normalized by the rotor rotation rate 
$\varOmega _R$ and the wavenumber corresponds to an azimuthal wavenumber (
$k=1$ means an oscillation for a rotation of 
$2{\rm \pi}$).
Growth rate as a function of the (azimuthal) wavenumber 
$k$ for 
$\varepsilon = 0.05$, 
$\eta =0.02$, 
$N=2$ and 
$\zeta _p = 20{\rm \pi}$ in different rotor regimes (solid line: 
$\lambda =-20$; dashed line: 
$\lambda =40$; dotted line: 
$\lambda =9.5$; dash-dotted line: 
$\lambda =8$). (a) Growth rate 
$\sigma ^{*}=\sigma /\varOmega _R$ (normalized by 
$\varOmega _R$). (b) Growth rate 
$\sigma 2 h_{FW}^{2}/(N^{2}\varGamma )$ (normalized by the characteristic advection time 
$2 h_{FW}^{2}/(N^{2}\varGamma )$). Dashed red line is the theoretical prediction by Gupta & Loewy (Reference Gupta and Loewy1974) for 
$h/R = 0.35$ and 
$\varepsilon =0.05$.
In order to compare the characteristics of the instability in the different regimes, the growth rate has to be non-dimensionalized by an advection time of the vortex structure. We use, as we did above, the advection time 
$2 h_{FW} ^{2}/(N^{2}\varGamma )$ obtained for a double array of 2-D point vortices (time needed to be advected by their separation distance 
$h_{FW}/N$). With such a rescaling, we observe in figure 20(b) that the stability curves obtained for the different regimes superimposed on a single curve corresponding to the theoretical prediction calculated by Gupta & Loewy (Reference Gupta and Loewy1974) for infinite helices. In principle, the theoretical prediction depends on 
$h/R$ and 
$\varepsilon$ but for the parameters we have considered this dependency is not significant, so a single curve has been plotted. The most unstable mode is obtained for 
$k=1$ with a value of 
$\sigma \simeq ({\rm \pi} /2) \varGamma /(2 (h_{FW}/N)^{2})$. As already noticed by Quaranta et al. (Reference Quaranta, Bolnot and Leweke2015), it is worth mentioning that this maximum growth rate also corresponds to the maximum growth rate of the pairing instability for an array of point vortices of circulation 
$\varGamma$ separated by a distance 
$h_{FW}/N$ (see also appendix B).
The good agreement with Gupta & Loewy (Reference Gupta and Loewy1974) for all the cases we have studied is a confirmation that the evolution of the perturbation is not affected by the presence of the rotor as soon as we are a few radii behind the rotor.
4.3. Spatio-temporal behaviour
Although the temporal stability properties do not change, the way the wave packet expands and is advected does vary according to the regime. In figure 21 we display the local growth rate contours in the 
$(\zeta ,\psi )$ plane of four characteristic cases. The perturbation has been introduced in the far field at 
$\zeta _p/(2 {\rm \pi}) = 10$. For both the climbing regime (figure 21(a) for 
$\lambda = -20$) and the wind turbine regime (figure 21(d) for 
$\lambda = 4.3$) the perturbation is advected away as it grows. This is characteristic of convectively unstable regimes. For the two other regimes, the wave packet does not completely move away as it expands. For 
$\lambda = 11$, a part of it clearly moves upstream (see figure 21b). The flow is then absolutely unstable in that case. For 
$\lambda =6.2$, the wave packet amplitude remains constant at the impulse location (see figure 21c). It corresponds to a marginally absolutely unstable configuration.
Spatio-temporal evolution of the perturbation for 
$\varepsilon = 0.01$, 
$\eta =0.05$, 
$N=2$. Contours of the amplitude envelope (in log scale) in the 
$(\zeta ,\psi )$ plane for a impulse placed at 
$\psi =0$ in 
$\zeta /2{\rm \pi} = 10$; (a) 
$\lambda = -20$, (b) 
$\lambda = 11$, (c) 
$\lambda =6.2$, (d) 
$\lambda =4.3$.
To analyse quantitatively the wave packet evolution, it is interesting to consider the growth rate 
$\sigma _V=\sigma (V_{\psi })\approx \log (|z'(\psi ,\zeta _p + V_{\psi } \psi )|)/\psi$ in the frame moving at a given (dimensionless) speed 
$V_{\psi }$ from the impulse location 
$\zeta _p$. The corresponding plot is shown in figure 22(a). Of particular interest are the velocities 
$v_{\psi }^{-}$, 
$v_{\psi }^{+}$ which delimit the frame velocity interval in which the perturbation grows and the frame velocity 
$v_{\psi }^{max}$ at which the growth is maximum. These velocities correspond to the slopes of the green and red lines indicated in figure 21(b). One can see in figure 22(a) that 
$v_{\psi }^{max}$ is close to 1 in all the regimes. This is not surprising as the maximum growth rate is expected to be reached in the frame moving with the vortex elements. These elements move at a negative angular velocity close to 
$- \Omega _R$, corresponding to the dimensionless velocity 
$v_{\psi }=1$. The difference comes from the self-inducted angular motion of the helices. As observed in figure 12, the azimuthal induced velocity is positive in helicopter regimes, and negative in wind turbine regimes. This explains why 
$v_{\psi }^{max}$ is slightly smaller than 1 for helicopters, and slightly larger than 1 for wind turbines.
Maximum temporal growth rate in a frame moving at the angular speed 
$V_{\psi }$ (or equivalently at the non-dimensionalized axial velocity 
$V_z 2 {\rm \pi}/ (\varOmega _R h_{FW})= V_{\psi }$) for 
$\varepsilon = 0.01$, 
$\eta =0.05$, 
$N=2$ and 
$\zeta _p = 20{\rm \pi}$ in different rotor regimes (solid line: 
$\lambda =-20$; dashed line: 
$\lambda =11$; dotted line: 
$\lambda =6.2$; dash-dotted line: 
$\lambda =4.3$). (a) Growth rate 
$\sigma _V/\varOmega _R$ (normalized by 
$\varOmega _R$) versus 
$V_{\psi }$. (b) Growth rate 
$\sigma _V 2 h_{FW}^{2}/(N^{2}\varGamma )$ (normalized by the characteristic advection time 
$2 h_{FW}^{2}/(N^{2}\varGamma )$) versus 
$(V_{\psi } -1)h_{FW}^{2} \varOmega _R/(N\varGamma {\rm \pi})$. Dashed red line is the theoretical prediction obtained from an array of 2-D point vortices 
$\sigma ^{*} = {\rm \pi}/2(1-X^{2})$.
In the far field, as the pitch is constant and equal to 
$h_{FW}$, a 
$2 {\rm \pi}$ increases of 
$\zeta$ corresponds to an axial displacement of a distance 
$h_{FW}$. It follows that 
$V_{\psi }$ can also be considered as an axial speed 
$V_z$ normalized by 
$h_{FW} \varOmega _R/(2{\rm \pi} )$.This can be used to compare the results with a 2-D point vortex model. As shown in appendix B, for an array of point vortices of circulation 
$\varGamma$ separated by 
$h_{FW}/N$, the pairing instability develops in a frame moving at the velocity 
$V_{rel}^{*}=V_{rel}2h_{FW}/(N \varGamma )$ with respect to the vortex frame with a growth rate given by
\begin{equation} \sigma_V 2h_{FW}^{2} /(N^{2} \varGamma) = \frac{\rm \pi}{2} \left(1 - (V_{rel}^{*})^{2}\right) . \end{equation}
For the far-wake helices, the relative velocity is given by 
$V_{rel}^{*}=(V_{\psi } -1)h_{FW}^{2}\varOmega _R/(N \varGamma {\rm \pi})$. We have compared our numerical results with formula (4.7) in figure 22(b). We can observe that the 2-D point vortex model works quite well. A similar agreement was already observed for the instability developing in an array of vortex rings (Bolnot et al. Reference Bolnot, Le Dizès and Leweke2014).
When the Dirac impulse is placed close to the rotor plane, the spatio-temporal development of the perturbation becomes more complex, especially for the configurations where the vortices are on both sides of the rotor. A few illustrations of the behaviour of the displacement amplitude are shown in figure 23. Whereas for the classical cases (
$\lambda =-20$ and 
$\lambda =4.3$) a well-defined wavepacket forms and is advected away from the rotor. For the two other cases, a dominant part of the wavepacket is still advected away as in the far field, but a small part of the wavepacket is now advected upstream. The wake has then become absolutely unstable. This part of the wavepacket is now trapped in the close neighbourhood of the rotor in the region where the vortices move on the other side of the rotor. No simple model is able to describe the evolution of the wavepacket in those cases.
Spatio-temporal evolution of the perturbation for 
$\varepsilon = 0.01$, 
$\eta =0.05$, 
$N=2$. Contours of the axial displacement 
$|z'|$ (in log scale) in the 
$(z,t)$ plane for a impulse placed at 
$\psi =0$ in 
$\zeta /2{\rm \pi} = 3$; (a) 
$\lambda = -20$, (b) 
$\lambda = 11$, (c) 
$\lambda =6.2$, (d) 
$\lambda =4.3$. The same colour map is used for all the figures. The red diamond corresponds to the impulse location. The rotor is located at 
$z=0$.
5. Conclusion
In this work, we have provided a vortex wake model for a rotor in any incident normal wind. It applies to all the vertical flight regimes of a helicopter (ascent, hover, descent) as well as all the regimes of a wind turbine. This model, which was first introduced a century ago by Joukowski, is based on a simple description of the wake by tip and root vortex filaments for each blade. Contrarily to the classical studies where a helical geometry was in general prescribed for the wake, here we have solved the governing equations to obtain the geometry of the vortex structure. This has allowed us to account for the contraction of the wake in the helicopter regimes and for its expansion in the wind turbines regimes. In the helicopter regimes, we have been able to obtain solutions up to the regimes of rapid descent, for which the vortices move above the rotor before being advected below it. Similarly, in the wind turbine regimes, we have shown that for large tip ratios the vortices could also cross the rotor plane. The properties of these new types of solutions have been analysed in detail, by providing their geometric characteristics, their induced velocity field in the rotor plane and in the far wake as well as their thrust and power coefficients. The stability of the solutions has also been considered by analysing the linear impulse response. All the solutions have been found to be unstable with respect to a local pairing instability, and well described in the far wake by the theoretical predictions for uniform helices (Gupta & Loewy Reference Gupta and Loewy1974) and for an array of 2-D point vortices. However, the more complex solutions have been found to be more unstable, with perturbations propagating upstream (absolute instability) and exhibiting a non-trivial spatio-temporal behaviour close to the rotor.
It is clear that for those cases, the vortex structure is expected to be globally unstable. We suspect that it is replaced by a time-dependent solution corresponding to the so-called vortex ring state. This state is characterized by the presence of a ring structure in the rotor plane that is periodically shed (Drees & Hendal Reference Drees and Hendal1951; Green, Gillies & Brown Reference Green, Gillies and Brown2005; Stack, Caradonna & Savaş Reference Stack, Caradonna and Savaş2005).
Acknowledgements
The authors are grateful to T. Leweke and H. Quaranta for having provided the images shown in figure 10.
Funding
This work received support from the French Agence Nationale de la Recherche under the A*MIDEX grant ANR-11-IDEX-0001-02, the LABEX MEC project ANR-11-LABX-0092 and the ANR TWIN-HELIX project ANR-17-CE06-0018-01. E.D.V. also acknowledges DGA and Aix Marseille Université for his PhD grant.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Numerical convergence study
In this appendix, we analyse the convergence of the results with respect to discretization parameters that are the number of discretization points per turn, the number of turns in the computation domain and the number of turns used to describe the far-wake structure. We consider the highly deformed case obtained for the parameters 
$\lambda = -40$, 
$\varepsilon = 0.05$, 
$\eta = 0.02$ and 
$N = 2$ which corresponds to one of the most challenging configuration. Both the helical structure and the stability results are analysed.
In figure 24, we show the effect of the numerical parameters on the far-wake and the near-wake geometries. The effect on the far wake is analysed by considering the radius 
$R_{FW}$ and the pitch 
$h_{FW}$. For the near wake, we consider the position of the points of the first two turns. More specifically, we define the maximum distance 
$d_{NW}^{max}$ between the points of the first two turns of the helix compared to a well-converged case. For this analysis, the converged reference case has 60 points per turn, 40 turns for the computation domain and 40 turns for the far wake.
Effect of numerical parameters on the far-wake and near-wake geometries for 
$\lambda = -40$, 
$\varepsilon = 0.05$, 
$\eta = 0.02$ and 
$N = 2$. The default parameters are 30 discretization points per turn, 30 turns for the computation domain and 30 turns for the far-wake domain. (a,d,g) Radius 
$R_{FW}$ in the far wake; (b,e,h) pitch 
$h_{FW}$ in the far wake; (c,f,i) near-wake maximum distance 
$d_{NW}^{max}$. (a–c) Number of discretization points per turn; (d–f) number of turns in the computation domain; (g–i) number of turns in the far wake.
The effect of the number of the discretization points per turn was already analysed in Durán Venegas & Le Dizès (Reference Durán Venegas and Le Dizès2019) for periodic configurations. Similar results are obtained for semi-infinite configurations: when one changes the number of points per turn from 30 to 60, we observe a difference smaller than 
$0.1\,\%$ for 
$R_{FW}$ and 
$0.8\,\%$ for 
$h_{FW}$. For the near-wake structure, 
$d_{NW}^{max}/R_b$ is found to be 
$5\times 10^{-3}$.
Concerning the number of turns in the computation domain, 
$R_{FW}$ and 
$h_{FW}$ are found not to vary much above 10 turns. Between 20 turns and 40 turns, we obtain variations of 
$0.05\,\%$ for 
$R_{FW}$, of 
$0.2\,\%$ for 
$h_{FW}$ while 
$d_{NW}^{max}/R_b$ is 
$2\times 10^{-4}$. Concerning the number of turns in the far wake, the difference between taking 10 and 40 turns is 
$0.02\,\%$ for 
$R_{FW}$ and 
$0.08\,\%$ for 
$h_{FW}$. This difference is reduced by one order of magnitude if one compares 30 turns with 40 turns. In the near wake, 
$d_{NW}^{max}/R_b$ is found to be 
$10^{-3}$ when 30 turns are compared to 40 turns.
The effect of the numerical parameters on the instability growth rate is shown in figure 25. We observe that larger numerical parameters are needed to reach convergence for first two instability bands. For instance, if one takes 20 points per turn, the wake structure is well described, as well as the first instability band (see figure 25a), but the second instability peak is not properly reproduced. One needs 60 points per turn to obtain converged results for this instability band. Concerning the number of turns in the computation domain (figure 25b), no significant differences are observed when the number of turns varies from 15 to 40. The limitation actually comes from the advection of the perturbation; the domain has to be sufficiently large such that the wave packet has not reached the far wake before being completely developed. Finally, for the number of turns in the far wake (figure 25c), a quite good convergence is already obtained as soon as it is larger than 10. Above 20 turns, the difference becomes imperceptible (see figure 25c).
Effect of the numerical parameters on the growth rate as a function of the (azimuthal) wavenumber 
$k$ for 
$\lambda = -40$, 
$\varepsilon = 0.05$, 
$\eta = 0.02$, 
$N = 2$ and 
$\zeta _p = 16{\rm \pi}$. The default parameters are 60 points per turn, 30 near-wake turns and 30 far-wake turns. (a) Variation of the number of points per turn: 20 (dotted), 40 (dashed dotted), 60 (solid) and 80 (grey thick solid). (b) Variation of the number of turns in the near wake domain: 15 (dotted), 20 (dashed dotted), 30 (solid) and 40 (grey thick solid). (c) Variation of the number of turns in the far-wake domain: 5 (dotted), 10 (dashed dotted), 20 (solid) and 30 (grey thick solid).
Appendix B. Spatio-temporal development of the pairing instability in an array of 2-D point vortices
In this appendix, we review the stability properties of an infinite array of point vortices of circulation 
$\varGamma$ located at 
$z_p = p h$, where 
$p$ covers all the integers and 
$h$ is the separation distance. As shown by Lamb (Reference Lamb1932, § 156), this array is linearly unstable in the frame of the vortices with respect to a displacement of the form 
$z_p'=\delta \exp ({\textrm {i} k h p - \textrm {i} \omega t})$ with a complex frequency 
$\omega = \omega ^{*} \varGamma /(2h^{2})$ that is related to the wavenumber 
$k=k^{*}/h$ by
\begin{equation} \omega^{*} = \textrm{i} \frac{\rm \pi}{2} \left(1 - \left(\frac{k^{*}}{\rm \pi} -1\right)^{2}\right). \end{equation}
In a frame moving at the velocity 
$V_{rel}= V_{rel}^{*} \varGamma /(2 h)$ with respect to the vortices, the dispersion relation is
\begin{equation} \omega_V^{*}= \omega^{*} - k^{*} V_{rel}^{*} = \textrm{i} \frac{\rm \pi}{2}\left(1 - \left(\frac{k^{*}}{\rm \pi} -1\right)^{2}\right) - k^{*} V_{rel}^{*}. \end{equation}
As explained in Huerre & Monkewitz (Reference Huerre and Monkewitz1990), in order to find the maximum growth rate of a wave packet in this moving frame, one should consider the wavenumber 
$k_V^{*}$ that cancels the group velocity
\begin{equation} \frac{\partial \omega_V^{*}}{\partial k^{*}}(k_V^{*})= 0, \end{equation}which gives
\begin{equation} k_V^{*} ={\rm \pi}\left(1 + \textrm{i} V_{rel}^{*} \right). \end{equation}
The maximum growth rate of the wave packet in the frame moving at the velocity 
$V_{rel}^{*}$ is thus given by 
$\sigma _V^{*} =\Im m( \omega _V^{*}(k_V^{*}))$ with
\begin{equation} \omega_V^{*}(k_V^{*}) = \textrm{i} \frac{\rm \pi}{2}\left( 1 - (V_{rel}^{*})^{2}\right) - {\rm \pi}V^{*}. \end{equation}In dimensional form, we then obtain
\begin{equation} \sigma_V = \frac{{\rm \pi}\varGamma}{4h^{2}} \left( 1 - \left(\frac{V_{rel} 2 h}{\varGamma} \right)^{2}\right) \end{equation}
for the growth rate of a wavepacket moving at a velocity 
$V_{rel}$ with respect to the vortices.
