Flap/lead-lag computational investigations on NREL S809 airfoil

– This paper aims at getting some insight into aeroelastic instabilities on NREL S809 airfoil. The study is based on the solution of the Unsteady Reynolds Averaged Navier-Stokes (URANS) equations where turbulence is modelled using a two-equation SST k - ω model. Preliminary ﬂuid ﬂow computations are performed for a ﬁxed airfoil set at various angles of attack. The results showed good agreement with experimental data found in the literature. Flow induced-vibrations on the airfoil are then investigated using a strong coupling ﬂuid-structure interaction approach. The URANS equations are then deﬁned in Arbitrary Lagrangian-Eulerian (ALE) coordinates. Simulations are carried out for ﬂap/lead-lag free oscillations. All computations are performed in 2D domain in which the governing equations are solved considering an incompressible ﬂow and a Reynolds number based on airfoil chord Re = 3 . 33 × 10 5 to Re = 2 × 10 6 .


Introduction
In recent decades, wind power has grown significantly.Research firm Global Data reported that installed capacity increased from 2006 to 2012 at a compound annual growth rate of 25%.This means a jump from 73.6 GW in 2006 to 280.6 GW in 2012.This huge development has been possible because of the tremendous works which has been done to reduce wind turbine accident and mainly, wind turbine blade damage.Indeed, wind turbine blades receive the biggest part of the wind load and the unexpected aerodynamics loads, which create large vibrations of the blades, are one of the major causes of damage [1,2].Since the manufacturing cost of wind turbine blades is about 15% to 20% of wind turbine production cost [3], ensuring wind turbine blade stability is a crucial issue.
Aerodynamics is a necessary tool for predicting and modelling blade loads [4].Due to the great development of numerical tools, Navier-Stokes equation solvers i.e.CFD codes are often used for that purpose.To study the blade aeroelastic stability, CFD codes are coupled to a structural code that solves the airfoil dynamic equations to predict the airfoil behaviour.Among works that have been done to study the aeroelastic behaviour of wind turbine blades, we can mention Johansen [5] who used both EllipSys2D, an in-house 2D CFD code, and the Beddoes-Leishman model, to perform an aeroelastic study of a wind turbine airfoil with two and three degrees of a Corresponding author: a.bekhti@cder.dzfreedom.In CFD computations, turbulence was represented by the SST k-ω model.For stable conditions good agreement was obtained with both calculations, but for unstable conditions, only the CFD code predicted the correct response.
Chaviaropoulos et al. [2] carried out an aeroelastic stability analysis of a NACA 0015 airfoil under combined pitch/flap and flap/lead-lag motions.Different cases were simulated and all computations were performed for a Reynolds number based on the airfoil chord equal to 10 6 .Results highlighted the airfoil oscillatory behaviour and allowed to bring out parameters corresponding to stable cases and to unstable cases.
Numerical investigations are reported on aeroelastic stability of a wind turbine rotor blade that were carried out in the frame of KNOW-BLADE (Task-4) European project [6].The computations were performed in two and three dimensions.A bad agreement is shown at high wind speed.
Baxevanou et al. [7] presented an aeroelastic model to study the aeroelastic behaviour of wind turbine blades.They concluded that, for predicting NACA 0015 airfoil response, the turbulence model is as important than the coupling scheme or the space discretization scheme.
In reference [8], flow-induced vibrations have been investigated for two airfoils used on wind turbine blades, NACA 0012 and NACA 63 2 415.The method was based on a strong coupling approach for Fluid-Structure Interaction (FSI).The dynamic blade response due to the fluid forces was determined in time accurate sequence.Their computations were performed for laminar and turbulent flow.The resulting aerodynamic forces were different but similar airfoil response was obtained.
Further fluid-structure interaction (FSI) studies can be found but they were applied to helicopter blades, aircraft [9] or cylinders.Chen et al. [10] performed simulations for incompressible flow in which they calculate flow-induced vibrations of a cylinder and an airfoil.This paper focuses on the applied E-CUSP upwind scheme.Comparison of their numerical results with experimental ones attested of the robustness of the used scheme.This study was extended to 3-D transonic flutter prediction of a flexible wing, using Roe scheme [11].
Finite element method was also employed to study flow-induced vibrations on NACA 63 2 -415 airfoil: Computations were performed for laminar flow [12].Turbulence was also implemented using k-ω model [13].
As far as experimental investigation is concerned, we can mention those performed by Poirel et al. [14] and Poirel and Yuan [15].Using an aeroelastic apparatus composed of a wing moving with two degrees of freedom, the authors studied the effects of blade stiffness and aeroelastic axis position for a NACA 0012 airfoil.Their experiments were conducted at Reynolds number equal to 10 4 .
To the authors' knowledge, there is no wind turbine blade aeroelasticity study applied to a free vibrating NREL S809 airfoil.The aim of the present paper is to simulate flow-induced vibrations on the NREL S809 wind turbine airfoil using a strong coupling fluid-structure interaction approach.Then the dynamic response of the wind turbine blade is determined in time accurate sequence.Two degrees of freedom are considered, flap and lead-lag motions.The numerical simulations are performed using Code Saturne, an open source CFD code in which the computational structural program for the solution of the blade dynamic equations is implemented in user subroutines for both calculations coupling and grid updating.The viscous flow-induced vibrations on the airfoil are then simulated.Beforehand computations are applied to a fixed airfoil set at different angles of attack to have confidence in our modelling approach.

Numerical approach
All computations are carried out in a 2D domain for an incompressible and turbulent flow with the solution of the Unsteady Reynolds Averaged Navier-Stokes (URANS) equations.The fluid flow equations, the oscillating airfoil dynamic equations, the coupling approach and the moving mesh technique equations are summarised in the following.

Fluid equations
The fluid equations are described in ALE coordinates in which the grid is considered as a referential frame moving at an arbitrary velocity.Let Ω f ∈ R 2 a 2D spatial domain occupied by the fluid.The incompressible Reynolds Averaged Navier-Stokes equations in ALE form write as, in tensor notation: where t, is the time, x i the Cartesian coordinate of a point of Ω f , u i , the averaged velocity component in the i direction, u cj , the moving grid velocity, p, the pressure, ρ, the fluid density, and τ ij , the stress tensor components.

Turbulence modelling
The Reynolds stresses due to turbulence are described by the SST k-ω turbulence model.It was shown that this model is efficient for the simulation of wall bounded flow with pressure gradient [16][17][18].Detailed equations of this turbulent model can be found in [19].When the RANS equations are written in ALE formulation, the terms ∂(u i k)/∂x i and ∂(u i ω)/∂x i in the transport equations for k, the turbulent energy and ω, the turbulent specific dissipation write as: For turbulence equations, the applied wall function is a two equation model that involves the friction velocity of the fluid u * and a turbulent friction velocity u k (see description in Appendix A).

Computational domain
The computational domain is of type O-H, the airfoil being located inside a cylinder of R1 radius.The domain extends upstream to a distance of 10c and downstream to a distance of 30c, where c is the airfoil chord.The North and South boundaries are located at +10c and −10c respectively.The domain thickness is equal to c (front and back boundaries are located at z = 0 and z = c respectively (Fig. 1)).The computational domain dimensions are chosen so that free stream conditions are achieved at outflow, north and south boundaries and no incoming flow can occur at outflow boundary.

Boundary conditions
Inflow condition is applied at the West boundary.Given the free-stream velocity U ∞ (Fig. 1), the turbulent kinetic energy k ∞ and the turbulent dissipation rate ∞ are defined as:  where D h is the hydraulic diameter (D h = c = 1), I is the free stream turbulent intensity (I ∞ = 2%), C µ and κ are constant coefficients (C µ = 0.09 and κ = 0.42).At the airfoil surface, fluid velocity is equal to the airfoil speed: where u is the fluid velocity on the airfoil surface and v b is the airfoil velocity.

Aerodynamic force computation
Aerodynamic force resultant is obtained from the Navier-Stokes equations solution by integrating the pressure p and the shear stress over the blade surface.The lift and the drag coefficients C L and C D are computed using the following relations: and where F y and F x are the y and x components of the aerodynamic force, respectively.

Algorithm and schemes
All simulations are performed using Code Saturne.The equations are solved using the finite volume method with a fractional time step integration scheme.Variables are collocated at the cell centres.Rhie and Chow interpolation is used to stabilise pressure oscillations [20].The Centred scheme is used for the discretization of the convection-diffusion terms and temporal discretization is performed with the implicit θ-scheme.

Dynamic equations
It is assumed that the airfoil is a free vibrating body whose elasticity is described by a two degree freedom model.The airfoil moves in horizontal and vertical directions (Fig. 3) and airfoil displacements are caused by the aerodynamic lift (L) and drag (D) forces.The dynamic equations of the airfoil are defined from the Lagrange equations: where x, y, ẋ, ẏ, ẍ and ÿ are the airfoil lateral and vertical displacement, velocity and acceleration respectively, m is the airfoil mass, k x and k y are the bending stiffness, ζ x and ζ y are the damping coefficients ratio, and ω y = (k y /m) 1/2 are the airfoil natural frequencies.
The solution algorithm of the dynamic equations is based on the Newmark scheme (see Appendix B).As it is well known, first order explicit scheme does not conserve energy at the moving fluid-solid interface, so implicit algorithms are more suitable.Thus, an implicit internal coupling technique is applied in these simulations.

The moving mesh equations
The O-H block structured grid allows a meshing technique easy to implement.As the airfoil oscillates, the grid nodes move with a velocity u c derived from the solution of a Laplace smoothing equation [21]: and the new position of a node is: where x new and x old are the new and old node positions and Δt is the time step.λ is a virtual mesh viscosity defined by: λ = 10 10 for r ≤ R 1 (10) λ = 10 10 e −2(r−R1) + 1 for r > R 1 (11) r is the distance from cell centre to centre of the moving airfoil.It is assumed that mesh behaves as a viscous fluid; near the airfoil surface, high value of mesh fluid viscosity is imposed to a circular area within a radius R 1 (Eq.( 10)).Beyond the circle, mesh fluid viscosity varies according to an exponential law (Eq.( 11)).The outer domain vertices are stationary.This technique avoids large displacements for cells and nodes.Given an appropriate choice of R 1 , the resulting mesh distortions are small and the mesh quality is preserved.The mesh displacement is represented in Figure 4.

Results and discussion
To check our numerical model (used turbulence model, computational mesh and numerical schemes), preliminary computations are carried out for the simulation of the turbulent flow around a fixed S809 airfoil set at different angles of attack.Then the airfoil subjected to flow-induced vibrations behaviour is investigated.

Mesh resolution test
Simulations are carried out for three types of grid that differ mainly by the total cell number and cell number on the airfoil surface.For the grid "mesh 1", the number of cells at the airfoil surface is 140; for the grid "mesh 2", it is 220 and for the grid "mesh 3", N s is 320.
The resulting total numbers of cells are 23 850, 30 250 and 44 300 for the first, second and third grid respectively.The height of the first cell row around the airfoil is the same for the three meshes so that y + remains between 15 and 45.The test case is performed for a stationary airfoil set at a fixed angle of attack of 13 • and for a Reynolds number Re = 10 6 .
The mean lift and drag coefficients are computed with the obtained values from which the solution begins to stabilise.According to Figure 5 where transitional variations of lift and drag coefficients are depicted, mean values of C L and C D are computed from t = 1 to 3 s.
Computed results are compared to experimental data obtained at Ohio State University and DUT wind tunnels where tests were performed with transition free [22,23] (Table 1).
It is noticed that experimental data of DU T were obtained at a fixed value of Reynolds number Re = 10 6 (equal to the one of the present study) and in a low turbulence level wind tunnel.Whereas at OSU , the Reynolds number varies from 0.99 × 10 6 to 1.04 × 10 6 .
It appears that the S809 geometry is not correctly represented by 140 elements (results of mesh 1).Both lift and drag coefficients computed with meshes 2 and 3 are slightly different but the drag coefficient obtained with mesh 3 is somewhat closer to experimental data of DUT.However, computational time with mesh 2 is lower than that with mesh 3. Accordingly, mesh 2 was chosen for the following calculations as a compromise between accuracy of results and computational time.we see a slight drop in lift coefficient.That could be due to laminar separated bubbles observed at the trailing edge.C L still increases with α up to 15 • and then it suddenly decreases while the drag coefficient C D increases roughly.

Lift and drag forces
The static stall occurs.These results are in good agreement with experimental data of [22] that were conducted for Reynolds number Re = 10 6 .Similar results were obtained by [24].Using ANSYS-Fluent code, the later studied the S809 airfoil in similar flow conditions and their results were compared to Delft University Experiments.Thus the efficiency of our applied numerical approach is confirmed and could be carried out for the flap/lead-lag investigations.

Velocity contours
The velocity magnitude iso contours are represented in Figure 7 for the airfoil set at different angles of attack.As expected, flow separation position on the upper surface moves toward the leading edge as the angle of attack α increases (Fig. 7).For α = 12 • detachment covers 40% of the airfoil upper surface and for α = 18 • it covers more then 90%.The vortex shedding frequency f v obtained from FFT analysis of the lift coefficient times signal is described by the corresponding non-dimensional Strouhal number: Pellegrino et al. [24] reported that for the airfoil set at α = 0 to 20 • the flow is largely attached and so no periodic vortex shedding would be expected.Indeed, the Strouhal number computed for the airfoil set at α = 13 • is about 0.001 whereas for the airfoil set at α = 40 • and 50 • , the computed Strouhal numbers are 0.25 and 0.18 respectively (it is noticed that Pellegino et al. [24] found Strouhal numbers of 0.23 and 0.17 respectively for the airfoil set at the same angles of attack).

Stall-induced flap/lead-lag flutter
These numerical simulations are applied to an airfoil of a parked wind turbine.Two degrees of freedom are considered, vertical and horizontal displacements.Computations are performed using the following structural parameters: -Case 1: m = 10 kg, k x = 15 000, k y = 7500.-Case 2: m = 10 kg, k x = 1 500, k y = 750.These values correspond to typical wind turbine airfoil, where stiffness in edgewise direction is higher than that in flapwise direction.It is noticed that parameter values vary with the blade span.
The influence of initial airfoil attack angle, free-stream velocity and airfoil stiffness on airfoil displacements is investigated.Results are given as a function of dimensionless time (t * ) defined as: where U ∞ is the free-stream velocity, and c is the airfoil chord.

Effect of initial angle of attack
The effect of initial airfoil angle of attack is considered for α = −6 • , 0 • , 8  Δy(= y f − y 0 ) increases with the increase of the angle of attack from −6 • to 13 • and it decreases for α = 18 • (Fig. 13).This difference is attributed to the aerodynamic forces induced by fluid flow.
Figure 14 shows the turbulent energy k and the turbulent specific dissipation rate ω contours.These contours are depicted in the vicinity of the airfoil set at α = 8 • , 13 • and 18 • .It shows the massive turbulent separation on the body upper side when the airfoil is set at the high angle of attack located in the post stall region.It is shown also that both turbulent quantities increase at high incidence angles.
Non dimensional vortex shedding frequency is calculated from FFT analysis of the time accurate lift curves.The related non dimensional Strouhal number varies from St = 0.31 for α = 0 • to St = 0.21 for α = 18 • (Table 2).Thus comparing these results to those of the static airfoil at low angle of attack, shows that the vortex shedding frequency increases.However, it is found that the frequency of the vortex shedding is well above the natural frequencies of the blade.So resonance does not occur.Moreover, low amplitude.This reveals that small vortices are continously shed in the wake of greater frequency.

Reynolds number effect
Simulations are carried out for the airfoil set at various initial positions within the range of α between 8 • and 13 • , with free-stream velocity U ∞ varying from 5 m.s −1 to 30 m.s −1 .Stable airfoil oscillations are found for all flapwise cases.However the stabilization time (t s ) decreases when Reynolds number increases (Fig. 15a).This suggests that high free-stream velocity induces high aerodynamic damping.
As for airfoil y position, it is found, as previously, that final stable airfoil position (y f ) is different from initial airfoil position and that final stable airfoil position increases with increasing Reynolds number.Moreover, final airfoil position varies with airfoil stiffness.Figure 15b shows that Δy increases when stiffness increases.Amplitude of horizontal displacement is smaller than vertical displacement.It is shown on graphs presenting edgewise displacement that vibration amplitudes have a fixed value for the initial position of −6 • , 0 • , 8 • and 13 • and since the drag force suffered an abrupt increase for α = 18 • , the vibration amplitude increases.
Figures 16 and 17 represent the variation of the horizontal and vertical displacements versus horizontal and vertical velocities, respectively.It appears that the amplitude of horizontal displacement is lower than that of vertical displacement.This is due to higher value of k x , the airfoil stiffness, in stream-wise direction.Furthermore, vertical displacement and airfoil velocity increase with increasing Reynolds number (Fig. 16), whereas horizontal displacement and velocity decrease with increasing Reynolds number (Fig. 17).Thereafter, airfoil flow induced vibrations are simulated for a parked wind turbine airfoil.Two degrees of freedom are considered, horizontal and vertical displacements.The influence of various parameters such as freestream velocity, airfoil initial position, angle of attack and airfoil physical parameters on airfoil response is studied.It is found that final airfoil position is different from initial position i.e. the distance between both positions (final and initial) increases with the increase of the Reynolds number.On the other hand increasing the airfoil stiffness leads to the decrease of this distance.

Conclusion
Results of the present work allow us to distinguish stable cases from unstable ones that might induce structural damages to the wind turbine blades and that should be avoided by blade designers.Another phenomenon to avoid is related to the Strouhal number that should not approach the natural frequency of the blade.In perspective, as blade stiffness varies with the span, it is expected to perform numerical simulations of the flow around a 3D blade in flap/lead-lag motion and normal flutters to analyse the 3D effects and flow at the tip of the free vibrating blade.
u k q f ν and ν is the molecular kinematic viscosity.

Appendix B: Newmark method
The Newmark algorithm is given as: Then x i+1 , the new airfoil position, is calculated using the following equation: x i+1 = a 0 M +a 1 C +K F i+1 +(a 1 x i +a 4 ẋi +a 5 ẍi )C +(a 0 x i + a 2 ẋi +a 3 ẍi M (B.7) where:

Figure 6
Figure6shows the variation of lift and drag coefficients versus the angle of attack α.As expected, it appears that the lift coefficient increases with increasing α until 10 • when a first light stall occurs.At about α = 10 • ,

Force coefficients Mesh 1 Fig. 6 .
Fig. 6.Lift and drag coefficients with respect to the angle of attack-fixed airfoil.
Numerical simulations are carried out to investigate flow-induced vibrations on NREL S809 airfoil.These computations are based on the solution of URANS equations written in ALE formulation with turbulence represented by SST k-ω model.Preliminary computations are performed for the simulation of the fluid flow around a fixed airfoil set at fixed angles of attack.The resulting lift and drag coefficients are in good agreement with experimental data.When Re = 10 6 , static stall occurs at α = 15 • .

Table 1 .
Comparison of CL and CD for different meshes.
• , 13 • and 18 • .The computations are performed for a Reynolds number Re = 10 6 .Figures 8 to 12 show that the new stable airfoil position (y f ) is different from initial airfoil position (y 0 ).Moreover,

Table 2 .
St numbers for the free vibrating airfoil.