E Durán Venegas

Stéphane Le Dizès
email: ledizes@irphe.univ-mrs.fr

Christophe Eloy
email: eloy@irphe.univ-mrs.fr

A strongly

Keywords: Rotors, Helicopters, Vortices, Flexible Blade Model, Fluid-Structure Interaction

published or not. The documents may come

Introduction

Rotors have been widely studied, often empirically, with practical applications in mind. Many models have also been proposed to estimate rotor performance in different operational conditions. General reviews of these models can be found in Leishman's book [START_REF] Leishman | Principles of Helicopter Aerodynamics[END_REF] for helicopters and in Hansen's book [START_REF] Hansen | Aerodynamics of wind turbines[END_REF] for wind turbines.

In several applications, rotor blades are long and flexible structures. These characteristics make blades prone to deformation when subject to aerodynamics loads. For helicopters and large wind turbines, significant deformations are indeed observed close to the blade tip, even if the blade material is not particularly compliant. By contrast, for some drones, the material is so elastic that the blade shape drastically changes when the rotor is in operation.

The external flow properties and, specially, the rotation of the blades introduce an heterogeneous airload distribution on the blades. These heterogeneous loads deform the blades, and consequently change both the blade position and the induced wake. As the wake also affects the incident velocity seen by the blade, the fluid-structure problem is always strongly coupled. In this paper, instead of trying to solve this nonlinear problem by direct numerical simulation [START_REF] Bazilevs | 3d simulation of wind turbine rotors at full scale. Part I: Geometry modeling and aerodynamics[END_REF], our objective is to introduce few simplifications that will allow us to perform a parametric analysis at low computation cost.

The first simplification concerns the blade geometry. Due to the slenderness of the blades, it is natural to use beam theory [START_REF] Timoshenko | On the transverse vibration of bars of uniform cross-section[END_REF] that exploits this property. Different one-dimensional equations for the bending and twisting angles are obtained along the blade centerline depending on the asymptotic limit considered [START_REF] Carrera | Beam structure: classical and advanced theories[END_REF]. The most used models are the simple Euler-Bernouilli rod model [START_REF] Dias | A non-linear rod model for folded elastic strips[END_REF] and the nonlinear ribbon model [START_REF] Volovoi | Assessment of beam modeling methods for rotor blade applications[END_REF][START_REF] Dias | Wunderlich meet Kirchhoff: A general and unified description of elastic ribbons and thin rods[END_REF] with an intrinsic non-linear coupling between bending and twisting. Several authors have studied the effect of this coupling on the adaptability of the rotors to non-optimal operational conditions [START_REF] Lobitz | Load mitigation with bending/twist coupled blades on rotors using modern control strategies[END_REF][START_REF] Maheri | Performance prediction of wind turbines utilizing passive smart blades: approaches and evaluation[END_REF][START_REF] Bottasso | Optimization based study of bend twist coupled rotor blades for passive and integrated passive/active load alleviation[END_REF]. In the present work, we shall use the simple rod model.

One of the main difficulties of rotor modelling is to have a good description of the flow that goes through the rotor plane. The most popular strategy is to use a momentum balance to relate the thrust and torque felt by the rotor to changes of axial and angular momentum in the flow [START_REF] Glauert | Airplane propellers. Division L[END_REF]. This so-called Blade Element Momentum Theory is widely used for the design of wind turbines [START_REF] Hansen | State of the art in wind turbine aerodynamics and aeroelasticity[END_REF]. It has been progressively improved to account for the vortical wake, multiple blades or three-dimensional effects [START_REF] Sørensen | General momentum theory for horizontal axis wind turbines[END_REF]. However, it assumes a particular topology of the flow streamlines that is not always satisfied, especially in the context of helicopters. In particular, in hovering and in weakly descending regimes, Blade Element Momentum theory does not apply.

For these cases, a more precise description of the wake based on vortex methods has often been used [START_REF] Gupta | Comparison of momentum and vortex methods for the aerodynamic analysis of wind turbine[END_REF]. The wake is then composed of vortex elements (such as tubes or sheets) that can have a prescribed form or move according to the Biot-Savart law [START_REF] Leishman | Free-vortex filament methods for the analysis of helicopter rotor wakes[END_REF]. Free-vortex methods can be very precise but they are unstable and costly when the number of vortex elements is large. A good alternative, that we adopt in this work, is to use a reduced vortex model involving only a tip and a hub vortex from each blade [START_REF] Joukowski | Théorie tourbillonnaire de l'hélice propulsive[END_REF]. This so-called Joukowsky model offers a simple but relatively efficient way to describe the wake [START_REF] Gupta | Comparison of momentum and vortex methods for the aerodynamic analysis of wind turbine[END_REF].

The paper is organized as follows. In section 2 we present the model, introducing all its building blocks. The model is then applied to two different flexible rotors in section 3. The effects of blade flexibility and external flow conditions on the blade deformation and wake properties are then analysed. Finally, some conclusions are drawn in section 5.

The model

In this section, we provide a description of the model. The model can be divided into five parts: the calculation of the aerodynamic loads, a model for the deformation of the blades, a model for the near wake, the calculation of the wake structure and the fluid-structure coupling.

Aerodynamic loads

The aerodynamic forces acting on the blade are calculated using the two-dimensional lifting-line theory [START_REF] Phillips | Modern adaptation of prandtl's classic lifting-line theory[END_REF]. The effective velocity felt by a blade element at the radial coordinate r b , when the rotor is rotating at the angular velocity Ω in an external axial flow V ∞ , is given by

U = V ∞ + V ind z 2 + Ω + Ωind 2 r 2 b , (1)
where V ind z and Ωind are the axial and angular velocities induced by the wake. The inflow angle with respect to the rotor plane is given by

ϕ = arctan V ∞ + V ind z Ωr b + Ωind r b . (2)
If the pitch angle of the blade element is β (figure 1), we can then define the angle of attack by

α = β -ϕ. (3)
The aerodynamic lift f L and drag f D forces per unit length on the blade element are given by where C L (α) and C D (α) are the lift and drag coefficients, ρ is the density of the fluid and c is the chord of the element. The dependence of C L and C D on the Reynolds number is neglected, as well as any hysteretic behavior. The forces f L and f D are exerted perpendicularly and parallel to the velocity direction (figure 1). The aerodynamic force components perpendicular and parallel to the blade element direction are then given by

f L = 1 2 ρU 2 c C L (α), f D = 1 2 ρU 2 c C D (α). (4) (a) (b)
f ⊥ = f L cos α + f D sin α, f = f D cos α -f L sin α. (5)
In the present model, we shall consider azimuthally averaged quantities for the induced velocities. This averaging procedure is necessary to smooth out the strong velocity variations that are artificially introduced by our wake model (see §2.3).

Deformation model

The blade deformation is obtained using a beam model for the blade. Such a 1D model provides a simple way to compute bending and torsion from the forces exerted on the blade. We are looking for equilibrium configurations. These configurations satisfy the Kirchhoff equations [START_REF] Dias | Wunderlich meet Kirchhoff: A general and unified description of elastic ribbons and thin rods[END_REF], expressing the balance of forces and moments on the blade element located at the curvilinear coordinate s ∂T ∂s

+ f = 0, (6)
∂M ∂s + d 3 × T + m = 0. ( 7
)
where f and m are the external forces and moments per unit length, T and M the internal forces and moments and d 3 the vector tangent to the blade centerline. These equations are conveniently written in the local orthonormal frame (d 1 , d 2 , d 3 ) (figure 2a). These vectors d i , i = 1, 2, 3 satisfy

∂d i ∂s = ω(s) × d i (s), (8)
where the Darboux vector ω(s) is defined by

ω = κ 1 (s)d 1 (s) + κ 2 (s)d 2 (s) + τ (s)d 3 (s). ( 9
)
The components of ω in the local frame are the normal curvature κ 1 (s), the geodesic curvature κ 2 (s) and the torsion τ (s). The normal curvature κ 1 and the torsion τ can also be defined in terms of a bending angle θ and a torsion angle γ (figure 2b)

∂θ ∂s = κ 1 , ∂γ ∂s = τ. (10)
Equations ( 6) and [START_REF] Volovoi | Assessment of beam modeling methods for rotor blade applications[END_REF] have to be completed by the beam constitutive law

M = EIκ 1 d 1 + EJ 2(1 + ν) τ d 3 , (11)
where ν is the Poisson ratio, E the Young modulus, I and J the second moments of area in the direction d 1 and d 3 respectively. We have assumed uniform elastic properties of the blades and neglected deformations in the direction d 2 (κ 2 = 0). Moreover, in the following, we assume that the airfoil profile is the same all along the blade, meaning that we can write I = I * c 4 and J = J * c 4 , where I * and J * are dimensionless quantities. The cross section area of the airfoil can also be written as A = A * c 2 where A * is a dimensionless quantity.

The blade deformation problem is then defined by the parameters linked to the blade profile (C L , C D , I * , J * and A * ), one parameter characterizing the blade aspect ratio c * = c/R b and three parameters associated with the blade material

ν, E * = E ρ b gR b , ρ * b = ρ b ρ , (12)
where R b is the length of the blade, ρ b and ρ the density of the blade and of the fluid respectively, and g the gravitational acceleration. The parameter E * compares the elastic forces with gravity. It is also useful to introduce the Froude number Fr that compares gravitational forces with centrifugal forces:

Fr 2 = R b Ω 2 g . (13)
In order to obtain two differential equations for θ and γ, we manipulate equations ( 6)- [START_REF] Bottasso | Optimization based study of bend twist coupled rotor blades for passive and integrated passive/active load alleviation[END_REF]. For the equation on θ, we differentiate equation ( 7) with respect to s and we take its projection on the local direction d 1 . After imposing κ 2 = 0, we obtain for the first term:

d 1 • ∂ 2 M ∂s 2 = ∂ 2 ∂s 2 (M • d 1 ) + τ 2 (M • d 1 ) ∂ 2 ∂s 2 (M • d 1 ) , (14)
where we neglect the term implying τ 2 , as the torsion curvature is small compared to the deflection, which will be verified a posteriori. So the final expression will be:

∂ 2 (M • d 1 ) ∂s 2 + ∂T ∂s • d 2 = 0, ( 15
)
where κ 2 = 0 has been used again to obtain the second term and neglect the term coming from the external 70 moments.

From the constitutive law [START_REF] Bottasso | Optimization based study of bend twist coupled rotor blades for passive and integrated passive/active load alleviation[END_REF] and using the definition [START_REF] Maheri | Performance prediction of wind turbines utilizing passive smart blades: approaches and evaluation[END_REF], the first term of equation ( 15) can be written as

∂ 2 (M • d 1 ) ∂s 2 = E I ∂ 3 θ ∂s 3 + 2 ∂I ∂s ∂ 2 θ ∂s 2 + ∂ 2 I ∂s 2 ∂θ ∂s . (16)
The second term of ( 15) is the projection on the direction d 2 of the forces exerted on the blade (see equation ( 6)). They include the aerodynamic force perpendicular to the blade, f ⊥ , the centrifugal force and the weight of the blade per unit length. In the following, the rotor axis will be assumed to be aligned with the direction of gravity, as for a helicopter in vertical flight. The projection on d 2 of the centrifugal force and weight can then be written:

f cent • d 2 = Aρ b r b (s)e r • d 2 = Aρ b Ω r b (s) sin θ = Aρ b Ω R b s 0 cos θds sin θ, (17)
f w • d 2 = Aρ b ge z • d 2 = Aρ b g cos θ sin(β + γ). (18)
Developing equation ( 15) with ( 5), ( 16)-( 18) we obtain the final equation for θ. The equation for γ is obtained by projecting equation ( 7) on d 3 with ( 10) and [START_REF] Bottasso | Optimization based study of bend twist coupled rotor blades for passive and integrated passive/active load alleviation[END_REF]. Finally, the two dimensionless differential equations for the bending and torsion angles θ and γ are:

ρ * b E * I * Fr 2 c * 4 θ + 8c * c * 3 θ + 4(3c * c * 2 + c * c * 3 )θ - 1 2 c * U * 2 (C L (α) cos α + C D (α) sin α) +ρ * b A * c * 2 s 0 cos θds sin θ + ρ * b A * c * 2 Fr 2 cos θ sin(β + γ) = 0, (19)
ρ * b E * J * 0 2(1 + ν)Fr 2 c * 4 γ + 4c * 3 c * γ - 1 2 c * 2 U * 2 [(C L (α) cos α + C D (α) sin α)δ ac cm + C m,ac (α)] = 0, ( 20
)
where the prime denotes differentiation with respect to s. The last two terms of the bending equation ( 19) correspond to the centrifugal force and the weight of the blade. The second term of the torsion equation ( 20) corresponds to the aerodynamic moment. The quantity δ ac cm is the distance (non-dimensionalized by c) from the center of mass to the aerodynamic center of the airfoil, typically located at a distance of c/4 from the leading edge for subsonic flows. The coefficient C m,ac is the pitching moment coefficient of the airfoil, which is in general constant for small angles of attack and equal to zero for symmetric airfoils. The angle of attack α(s) and the normalized incident velocity U * (s) = U (s)/(ΩR b ) are obtained from the wake model discussed in the next two subsections. These equations have to be integrated with the following boundary conditions

θ(0) = θ (1) = θ (1) = 0 , γ(0) = γ (1) = 0, (21)
that correspond to a "clamped" condition at the blade root, and a "free" condition at the blade tip.

In order to analyze the blade deformations in a more intuitive way, we also introduce the bending deflection function f θ (s) that measures the distance of each blade element from its original undeformed position.

Joukowski model

To describe the wake, we use a Joukowski model [START_REF] Joukowski | Théorie tourbillonnaire de l'hélice propulsive[END_REF]. In this model, the wake of each blade is composed of a bound vortex on the blade and two free vortices of opposite circulation Γ, emitted from two positions r t and r h close to the tip and hub of the blade. The points of emission and the vortex circulation Γ are computed from the circulation profile Γ(s) of the blade which is obtained from the local 2D Kutta-Joukowski formula

Γ(s) = 1 2 U (s)c(s)C L (α(s)). ( 22
)
The vortex circulation Γ is defined as the maximum value of Γ(s), while r tip = 0, i.e. the emission is on the rotor axis. This hypothesis is justified if we consider that the hub vortices emitted by all the blades are expected to merge away from the rotor on the rotation axis. It actually corresponds to the initial model introduced by Joukowski [START_REF] Joukowski | Théorie tourbillonnaire de l'hélice propulsive[END_REF]. For the other point of emission, it corresponds to the tip if the circulation reaches its maximum at this point. In general, we shall be in this situation, which is illustrated in figure 3

Wake structure

The wake structure is computed using a free-vortex method [START_REF] Leishman | Principles of Helicopter Aerodynamics[END_REF]. Each vortex is discretized in small segments of vorticity that move as material lines in the fluid according to

dξ dt = U (ξ) = U ∞ + U ind (ξ), ( 23
)
where ξ is the position vector of the vortex filament, U the velocity field, composed of the external velocity U ∞ and the induced velocity U ind (ξ) generated by the vortex filaments. When the rotor has N blades, the induced velocity deduced from the Biot-Savart law [START_REF] Saffman | Vortex dynamics[END_REF] is given by

U ind (ξ) = Γ 4π N j=1 (ξ j -ξ) × dτ j |ξ j -ξ| 2 , ( 24
)
where the integrals cover each vortex filament defined by its position vector ξ j and tangent vector τ j . This expression is also used to obtain the mean induced velocity on the blade defined in equation [START_REF] Leishman | Principles of Helicopter Aerodynamics[END_REF]. On the vortex filament, the Biot-Savart integral is singular, the self-induced velocity diverges. To avoid this singularity, a small but finite core size a is introduced. This allows us to obtain the self-induced velocity by the so-called cut-off method [START_REF] Saffman | Vortex dynamics[END_REF]. In the present study, we assume an identical and constant core size for all the vortices. In the present work, we do not consider any temporal evolution of the rotor system. We focus on regimes where the rotor is equilibrium with its wake. We therefore look for helical wake structures that are stationary in the rotor frame. This means that the vortices should be steady in the frame rotating at the angular velocity Ω of the rotor. This condition of steadiness implies a condition for the velocity field with respect to the vortex structures: it should everywhere be tangent to the structure. This condition reads

U ∞ (ξ m j ) + U ind (ξ m j ) × τ m j = 0, ( 25
)
where the external field U ∞ is the sum of the external axial velocity V ∞ e z and the azimuthal velocity -rΩe ψ . This condition can also be written for each vortex filament as

dr w dψ = V ind r Ω ind -Ω , dz w dψ = V ind z + V ∞ Ω ind -Ω ( 26
)
where r w (ψ) and z w (ψ) are the radial and axial positions of the filament. Concerning the boundary conditions, it is imposed that the free vortices are attached to their point of emission (the axis center and in are the radial and axial coordinates of the tip emission point. Far from the rotor, the wake is assumed to become homogeneous. In other words, we assume that each filament becomes a uniform helix after a certain number of turns, typically of order 15. The impact of this far wake is then computed by considering 8 additional turns of perfect helices. As shown by Ali & Abid [START_REF] Ali | Helical vortex wake: How far is the infinity?[END_REF], this is sufficient to have a good approximation of the effect of infinite helices.

It is worth mentioning that the wake structure deduced from the Joukowski model only depends on the number of blades N and three dimensionless parameters:

λ = r tip b Ω V ∞ , η = Γ (r tip b ) 2 Ω , ε = a r tip b , (27)
where λ is known as the tip-speed-ratio and η and ε represent the vortex strength and the vortex core size respectively.

Fluid-structure interaction

To obtain the circulation profile of the blade, we must determine the flow in the rotor plane. But this flow, especially the induced flow, also depends on the circulation and emission point of the vortices. The wake/rotor problem is thus always strongly coupled. The solution to this problem is here obtained through an iterative procedure following a fluid-structure interaction loop (see figure 4).

In practice, the problem is first solved for a rigid rotor. The first loop is performed by neglecting the induced velocity, so that first estimates for the circulation profile and the resulting wake can be computed from the external velocity only. This first loop also gives an estimate of the induced velocity, so that a second loop can be performed with the full velocity field. The process is repeated until convergence.

The flexible case is treated by considering the rigid rotor configuration as guess value. The difference with the respect to the rigid case is that the circulation profile is now obtained after having computed the new position of each blade element using the beam model (section 2.2). As for the rigid case, the process is repeated until convergence.

For both rigid and flexible rotors, convergence is typically obtained in 5 or 6 loops.

Applications

Description of the rotors

In this section, we apply the model to two different two-bladed rotors, named here rotor A and rotor B. For both rotors, we consider a same NACA0012 profile for the blade cross-section. The characteristic constants for this profile are I * = 0.0033, J * = 3.28 and A * = 0.0822. For the aerodynamic coefficients C L and C D , we assume the functions given in [START_REF] Chritzos | Aerodynamic characteristics of NACA 0012 airfoil section at angles of attack from 0 to 180[END_REF] and shown in figure 5 In the following section, we shall vary the operational conditions of the rotors (i.e. the angular rotation Ω and the external wind speed V ∞ ) that affect the tip speed ratio λ and the Froude number F r. The material properties of the blade will also be varied such that the effect of the dimensionless parameters ρ * b and E * will be considered. However, the Poisson ratio will be kept fixed and equal to ν = 0.5.

Except for the twist angle of the rotor A blade, the geometry of the two rotors without external forces will not be varied. We shall also not analyse the effect of the vortex core size that is fixed to ε = 0.01.

Results

In this section, we first analyse the effect of the tip-speed ratio λ on the wake characteristics and blade shape of rotor A. We fix the Young modulus E * , the density ratio ρ * b , the Froude number Fr, and vary λ by changing the external fluid velocity.

In figure 6, we have illustrated the radial trajectories of the vortices in the wake when λ is varied. In the rotor frame, the external wind is going downwards when λ is negative and upwards when it is positive. For a helicopter, this means that negative values of λ correspond to climbing flight, and positive values of λ to descending flight. Hover is associated with an infinite value of λ. Normal flight situations of a helicopter are shown in figure 6(a). For this case, the wake is contracting and moves downwards. The contraction increases when λ becomes more negative, that is when the climbing speed decreases. This contraction process continues when the climbing speed vanishes and changes sign, i.e. when we move to a slow descent flight regime corresponding to large positive values of λ. However, for this regime, the vortices are emitted upwards and therefore cross the rotor plane before going downwards. For smaller values of λ, a downward wake ceases to exist, and we jump to another type of solution shown in figure 6(b): the wake expands and goes upwards. This situation corresponds to the so-called windmill brake regime of helicopters. If gravity was not aligned with the rotor axis, it would correspond to the wind turbine regime. As expected, the stronger the external wind (i.e. the smaller λ), the less expanding is the wake.

The intervals of λ where downward and upward wakes exist do not seem to overlap. For all the parameters that we have considered, we have found an interval of λ where both solutions cease to exist. The limits of this interval are different for each case. We suspect that this is reminiscent of the Vortex Ring State (VRS) [START_REF] Drees | The field of flow through a helicopter rotor obtained from wind tunnel smoke tests[END_REF]. In figure 7, we have plotted the circulation profile obtained on the blade for the configurations shown in figure 6. For normal flight conditions (figure 7(a)), the circulation profile is found not to vary much. It is mainly associated with the blade rotation that prescribes a linear dependence of the circulation on the radial coordinate. The small bump observed for small values of s/R b is associated with the lift coefficient crisis obtained for small angles of attack [see figure 5(b)]. The small overshoot of circulation observed when λ = 15 corresponds to an induction effect of the vortices that are above the rotor plane for s/R b > 0.7. In the windmill brake regimes (figure 7(b)), a larger effect of variation is observed, with an increase of 40 % of the total circulation when λ changes from 3 to 5.7. This can be understood by the larger relative contribution of the axial wind in the total velocity for these cases.

The effect of λ on the blade deformation is shown in figure 8. The material chosen for the blade is weakly flexible so the deformation in terms of torsion angles (left plots) and bending angles (central plots) remains small. For normal flight conditions (upper plots), the deformation increases as the climbing velocity decreases. The largest deformation is reached for the slow descending regime. In the windmill brake regime (lower plots), the opposite behavior is observed : the deformation diminishes when the descent speed decreases. Both behaviors are in agreement with the variation of circulation with respect to λ. Larger circulation leads to larger deformation. Note however that the variations of γ and θ along the blade are much smoother and do not exhibit the bumps and jumps observed in the circulation profile.

The effect of the other parameters on blade deformation for a typical climbing regime of rotor A (λ = -20) is analyzed in figure 9. The torsion angle γ remains always much smaller than the deflection angle θ. However, both are similarly affected by variations of E * (upper plots), ρ * b (central plots) and Fr (lower plots): twisting and bending increase when E * or ρ * b decreases or when Fr increases. These similarities can be understood by going back to equation [START_REF] Saffman | Vortex dynamics[END_REF]. For the parameters of figure 9, weight and centrifugal forces remains small.

If the corresponding terms (third and fourth terms in equation ( 19)) are neglected, the equations for θ and γ become dependent on the parameters E * , ρ * b and Fr 2 through the combination ρ * b E * /Fr 2 only. It is then obvious that decreasing ρ b by a factor of 10 is equivalent to decreasing E * by the same factor, or to increasing Fr 2 by a factor of 10.

In figure 9, the deformation is relatively small. Highly deformed cases can nevertheless be calculated in the same way. In figure 10 we illustrate such a configuration by considering a very compliant material for rotor A. We observe in figure 10(c), that the deflection reaches 50% of the blade length when we consider the flow conditions λ = -10, ρ * b = 1 and Fr 2 = 100. The same rotor at rest is however almost undeformed as seen on this figure. We see also that, even for this very deformed case, the torsion remains very small compared to the deflection. This agrees with the assumption made in equation [START_REF] Sørensen | General momentum theory for horizontal axis wind turbines[END_REF], where the terms involving τ 2 were neglected.

So far, we have considered the simple geometry of rotor A. There is no difficulty to consider the more complex geometry of rotor B. In figure 11, we compare the circulation profile and blade deformation of both rotors for the same flow conditions, and the same material. We have chosen a value of λ close to conditions for which the rotor B has been designed. In figure 11(a), we observe that the circulation profile associated with rotor B is almost constant for s/R b > 0.4 in agreement with the design properties [START_REF] Quaranta | Long-wave instability of a helical vortex[END_REF]. The circulation profile obtained by rotor A is completely different. It gives a smaller vortex circulation but it gets more deformed by the flow (see figure 11(b-d)).

Comparison with experimental data

In this section, we compare our model with experimental data for a wind turbine rotor and a helicopter rotor in hover. For both cases, the blades are rigid. The data for the wind turbine rotor are taken from the MEXICO project [START_REF] Schepers | Analysis of mexico wind tunnel measurements: Final report of iea task 29[END_REF]. This rotor has variable chord and pitch angle and resembles rotor B. In the left column of figure 12, we have plotted the wake geometry, the lift force along the blade, and the predicted circulation profile for a tip-speed ratio λ = 6.67. The experimental data are indicated by symbols. We observe that there is a very good agreement for both the wake geometry and the lift force. This is not surprising as the circulation profile is mainly flat (lower left plot). The Joukowski model is expected to work for such a case. The good agreement therefore validates the numerical procedure.

The data for the helicopter rotor are extracted from a two-bladed rotor experiment from NASA technical rapport TM81232 [START_REF] Caradonna | Experimental and analytical studies of a model helicopter rotor in hover (TM 81232)[END_REF]. The rotor has in this case an untwisted constant-chord blade with a uniform pitch angle β = 12 • as rotor A. Similar quantities as for the wind turbine rotor are plotted in the right column of figure 12. For this rotor, the agreement with the model is not as good as for the wind turbine case. In the model, the wake is more contracted and the lift is underestimated in the inner part of the blade and weakly overestimated near the tip. As in the wind turbine case, a better agreement near the tip would be obtained if a tip correction was implemented in the model (see Sørensen's book [START_REF] Sørensen | General momentum theory for horizontal axis wind turbines[END_REF]). The circulation profile obtained from the model is shown in figure 12 (lower right plot). It strongly varies with the radial position. Note in particular that it becomes negative for r/R < 0.4. In principle, with such a circulation profile, we expect, from Prandtl lifting-line theory, vortex sheet shedding all along the blade, and the formation for each blade of three vortices: an inner vortex at r/R ≈ 0.2, a tip vortex at r/R ≈ 1 and third opposite sign vortex close to r/R = 0.6 resulting from the roll-up of the vortex sheet. It is the merging of the root and middle vortices of all the blades that is supposed to give the hub vortex of Joukowski model. This complicated roll-up and merging process is not considered. It may naturally affect the induced flow in the rotor plane. For the rotor in hover, the induced flow becomes the dominant part of the flow close to the rotor center. It is therefore in this region that the largest inaccuracy can be expected. The discrepancy observed in the wake geometry (upper right plot of figure 12) may not be associated with the model. Indeed, the experimental rotor is fixed on a big shaft aligned with the wake axis. We suspect that the shaft limits the contraction of the wake. It is not excluded that it could also influence the lift forces. For these reasons, it is difficult to have a definite opinion on the errors associated with the model for this case.

Conclusion

We have presented a model providing the equilibrium state of a flexible rotor in vertical flight conditions. The solution is stationary (in the frame rotating at the rotor angular speed), but takes into account the strong coupling between the rotor geometry and its wake. A simple but complete description of the blade deformation and of the wake generated by the rotor has been provided. A parametric study has also been performed for a simple rotor geometry in order to quantify the effect of external flow conditions and material properties. We have shown that the model is able to describe rigid as well as very flexible rotors. It also works in operational conditions where classical momentum theory does not apply. In particular, we have shown that we can describe rotors in slow descending regimes of helicopter flight, where the vortices created by the rotor move above the rotor plane before being advected downwards.

We have compared the rigid version of the model to experimental data for two rotors. We have found a good agreement for the lift distribution and the wake geometry for a wind turbine rotor exhibiting an almost constant circulation distribution along its blades. But further comparisons are needed to fully assess the validity of the model in hover or for other rotors.

It is worth emphasizing that we have only consider steady solutions. Strongly unsteady regimes such as the Vortex Ring State [START_REF] Stack | Flow visualizations and extended thrust time histories of rotor vortex wakes in descent[END_REF][START_REF] Quaranta | Instabilities in a swirling rotor wake[END_REF] are outside the scope of the present work. Furthermore, we have no information on the stability of the solutions. We know that the wake is a priori unstable with respect to long-wavelength instability [START_REF] Widnall | The stability of a helical vortex filament[END_REF][START_REF] Gupta | Theoretical analysis of the aerodynamic stability of multiple, interdigitated helical vortices[END_REF][START_REF] Quaranta | Long-wave instability of a helical vortex[END_REF]. But other instabilities associated with the blade flexibility, such as flutter, may also be present [START_REF] Eloy | Flutter of a rectangular plate[END_REF][START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF].

Figure 1 :

 1 Figure 1: Blade element airfoil in generic climbing (a) and descending (b) flight cases.

Figure 2 :

 2 Figure 2: Blade schemes showing the local orthonormal frame (a) and bending and torsion angles (b).

b

  and r hub b are the centers of mass on either side of the point of maximum of Γ(s) of the quantity ∂Γ/∂s. For the point of emission close to the hub, we shall assume r hub b

  (a).

Figure 3 :

 3 Figure 3: (a) The Joukowski model. (b) Illustration of the discretization into straight segments of the tip vortex.

Figure 4 :

 4 Figure 4: Fluid-structure loop scheme.

Figure 5 :

 5 Figure 5: (a) Blade geometry of rotor A (upper plot) and rotor B (lower plot). (b) Lift coefficient C L (dashed line) and drag coefficient C D (solid line) with respect to the angle of attack α for a NACA0012 profile, from [21]. (c) Chord c * and twist angle β (in degree) along the rotor B blade. Adapted from [22].

Figure 6 :

 6 Figure 6: Radial position of the wake in the axial direction for different tip-speed-ratio values of rotor A. (a) Helicopter regimes for β = 40 • . Climbing flight (λ = -10, λ = -20), hover (λ = ∞), and weakly descending flight (λ = 15). (b) Wind turbine (or windmill brake) regimes for β = 15 • . Other parameters are E * = 10 7 , ρ * b = 100, Fr 2 = 1000.

Figure 7 :

 7 Figure 7: Circulation profile obtained along the blade of rotor A for different tip-speed-ratio values. (a) Helicopter regimes (β = 40 • ); (b) Wind turbine regimes ( β = 15 • ). Other parameters are the same as in figure 6.

7 Figure 8 :

 78 Figure 8: Blade deformation of rotor A for the helicopter regimes with β = 40 • (upper plots) and for the wind turbine regimes with β = 15 • (lower plots). Left: Torsion angle γ; Center: Bending angle θ; Right: Bending deflection function f θ . The other parameters are the same as in figure 6.

Figure 9 :

 9 Figure 9: Effect of the other parameters on blade deformation of rotor A for β = 30 • and λ = -20. The default parameters are E * = 10 6 , ρ * b = 100 and Fr 2 = 100 Left column: Twist angle γ; Central column: Bending angle θ; Right column: Bending deflection function f θ . Upper line: Young Modulus effect. Middle line: Blade density effect. Lower line: Froude number effect.

Figure 10 :

 10 Figure 10: Highly deformable case. Rotor A with a twist angle β = 30 • for E * = 10 6 , ρ * b = 1 and λ = -10. Solid line: Fr 2 = 100; Dashed line: Fr = 0 (rotor at rest). (a) Torsion angle γ; (b) Bending angle θ; (c) Bending deflection function f θ . (d) Three-dimensional visualization of the wake and the deformed rotor.

Figure 11 :

 11 Figure 11: Effect of the blade geometry. Solid lines: rotor A with β = 20 • ; Dashed lines: rotor B. (a) Circulation Γ along the blade; (b) Torsion angle γ; (c) Bending angle θ; (d) Bending deflection function f θ . The other parameters are λ = 3.5, E * = 10 7 , ρ b = 100, Fr 2 = 100.

Figure 12 :

 12 Figure 12: Comparison of the model (solid line) with two experimental cases (•): A wind turbine configuration from MEXICO project [24] at λ = 6.67 (left column) and and a rotor in hover with an untwisted constant-chord blade at a constant pitch angle 12 • from NASA [25] (right column). Upper line: Radial position of the tip vortices in the axial direction. Central line: Lift force along the blade. Lower line: Circulation profile along the blade.

Acknowledgments

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. EDV also acknowledges support from DGA and Aix Marseille Université for his PhD grant.