Numerical simulation of the flow around oscillating wind turbine airfoils Part 2 : Free vibrating airfoil

ABSTRACT The present paper is Part 2 of a two parts paper on flow around vibrating wind turbine airfoils. The first part of the paper dealt with a forced oscillating airfoil. Part 2 focuses on free vibrating airfoils. The flow induced vibrations on two airfoils used for wind turbine blades are investigated by applying a fluid structure interaction approach. A commercial Computational Fluid Dynamics (CFD) code is coupled to a computational structural program that solves the dynamic equations of the airfoil oscillations. The fluid governing equations are described in the Arbitrary Lagrangian Eulerian coordinates and solved with a moving mesh. A straightforward meshing technique is implemented in a subroutine called by the CFD code at each time step for updating the grid. The method is applied to a free pitch oscillating airfoil and to combined pitch and vertical oscillations known as the flutter instability.


INTRODUCTION
Wind turbines aeroelastic stability problems begun to appear since 1990 with the development of larger blades.Most aeroelastic stability analyses devoted to wind turbines were performed using traditional aeroelastic design tools where the aeroelastic codes contain an aerodynamic part and a structural part: (i) the blade force and torque being first computed for given parameters (wind speed, rotational velocity …) with engineering-type dynamic stall models, where more often, semi-empirical models initially developed for helicopter rotor blades have been modified and adapted for airfoil sections used on wind turbines [1] and then (ii) the dynamic rotor equations being solved to describe the response of the wind turbine [2,3].However, with the progress of computational facilities, improved analysis of the mechanical behaviour of a rotating wind turbine is now possible by a Fluid Structure Interaction (FSI) approach.In these techniques, the aerodynamic forces are computed from the solution of the time-accurate Navier-Stokes equations.The calculated forces are then used for the solution of the rotor dynamic equations to determine the response of the structure at each time step.FSI techniques have been widely used in many industrial problems [4,5] and for aerospace applications where most papers are applied to compressible flows [6][7][8].Also, it was shown that compressibility effect plays significant role in aeroelastic stability limits.As wind turbines operate in an incompressible environment, specific studies have to be performed for wind turbine blades.To the author's knowledge, the main contribution with an FSI approach applied to wind turbine blades was made under the European project KnowBlade where Ellipsys3D, an in house computational Fluid Dynamics (CFD) code, and a structural code were used to simulate flap-lead/lag vibrations of a wind turbine blade [9].Recently, Svacek et al [10] presented the study of the classical Flutter by a FSI approach where the incompressible fluid equations were solved by the finite element method, using an ALE formulation of the Navier Stokes equations.
In this study, flow induced vibrations on airfoils are simulated with FSI computations using a commercial CFD code based on the finite volume method.In part one of this paper was implemented a straightforward technique for moving the mesh with a user subroutine.Then, time-accurate computations of the airfoil forces and torque based on the solution of the Navier Stokes equations expressed in ALE were applied to a forced oscillating airfoil.In this second part of the paper, the case of a free vibrating airfoil is considered.The aerodynamic forces, solution of the Navier-Stokes equations, are used to determine the response of the structure with the solution at each time step, of the airfoil dynamic equations.The coupling method of the CFD code with dynamic equations by this FSI approach is presented in the next section.Then the method is applied to a NACA 0012 airfoil in free pitch oscillations.Thereafter, the issue of the classical Flutter or combined pitch and vertical oscillations of a NACA 63 2 415 airfoil is considered.

NUMERICAL APPROACH
It is well known that wind turbines operate in a turbulent environments and it was shown in [11] that the lift and drag forces computed with the solution of the Reynolds Averaged Navier Stokes equations are sensitive to the free stream turbulence intensity.This aspect of the problem will be considered subsequently and there computations are carried for laminar flow as in the first part of the paper [12].The unsteady incompressible Navier-Stokes equations are described in ALE coordinates [13,14].The PISO algorithm is applied for the solution of the coupled pressure velocity equations with the first order differencing UPWIND scheme for the discretisation of the convection-diffusion terms.The temporal discretisation is performed with the implicit θ -scheme.The fluid forces, computed from the solution of the Navier Stokes equations, are used to determine the response of the structure with the solution, in a time accurate sequence, of the airfoil dynamic equations.The airfoil displacements are applied to the fluid mesh for the solution of the fluid equations and the determination of the blade forces at the next iteration.The dynamic equations of the oscillating airfoil, the coupling approach and the equations describing the moving mesh technique are summarized in the following.

THE DYNAMIC EQUATIONS
It is assumed that the airfoil is an elastic body and the airfoil elasticity is described by a simple two degrees of freedom model.The dynamic equations of the airfoil are defined from the Lagrange equations.In the case of an airfoil in pitch oscillations (or torsion around the elastic axis EO) and in flap-wise (oscillations in the vertical direction) as shown on Fig. 1, the non linear equations that describe the elastic airfoil motion write as: (1) m y C y S k y F y y y where m is the airfoil mass, α, y, , , and are the airfoil rotational and vertical displacement, velocity and acceleration respectively, F y is the vertical component of the aerodynamic force, M, the torsion moment, a, the distance between the point of application of F y and the elastic axis EO, with r G the distance between the centre of gravity and the elastic axis EO, k α and k y are the torsion and bending stiffness and C y , C α are defined as C α = 2J α ω α ζ α and C y = 2mω y ζ y , with J α inertia moment around the elastic axis EO, ζ α and ζ y , dampening coefficients and are the airfoil natural frequencies.It is assumed here that the aerodynamic force acts at the aerodynamic centre located at a quarter chords from the leading edge.
For the time integration of the dynamics equations, three algorithms are compared: a linear scheme, the Newmark algorithm and the Crank Nicholson scheme.Expressed for the solution of one degree of freedom equation: (3) these algorithms write as: • the linear scheme [15], The airfoil scheme. (4c) • the Newmark algorithm • the Crank Nicholson method [16] (6a) with: where is the time step and β 1 and β 2 , are constants of integration.The 3 algorithms are compared to the analytic solution of the following equation: given by the relation: (8) where the constants A 1 and A 2 are determined from the initial conditions and .It is shown on Fig. 2 that the solution obtained with the Crank Nicholson scheme coincides with the analytic solution.This algorithm is then used for the solution of the equations ( 1) and ( 2) that describe the airfoil oscillations.

THE COUPLING METHOD
The coupling method is based on the explicit scheme of Euler where the position of the computational domain at time t n+1 is determined with the relation: where q n is the displacement of the structure at time t n , and are the velocity of the airfoil at time t n and t n-1 respectively, is the time step and and are the coupling scheme coefficients.
Explicit methods are first order accurate and they do not conserve energy at the moving fluid -solid interface and it is well known that implicit algorithms are more suitable [13].However, these techniques are more computationally expensive and heavy to implement and to integrate in a commercial CFD code.Therefore they are not used in this study.

THE MOVING MESH TECHNIQUE
The moving mesh method is analogous to the meshing technique applied for the airfoil in forced oscillation [12].As previously, the airfoil is located in the centre of an O-H computational domain and the meshing method is based on algebraic interpolations: Given the cylindrical coordinates of the vertices in the 2D computational grid at the time t n , the Cartesian coordinates and and the airfoil displacements calculated with the solution of the dynamic equations, the displacement of the mesh is performed as follow:

•
The vertices of the circular sub-domain of radius R1 move at the airfoil velocity: • For the vertices of the annular sub-domain delimited by the circles of radius R1 and R2, the applied relations are: where , and NR 12 is the number of cells according to the radial direction, between R 2 and R 1 .

•
The vertices of the outer domain are stationary.
This meshing technique is implemented in a user subroutine called by the CFD code at the beginning of each time step.Our computations are applied to the test case corresponding to Re = 76 064, ω n = 57 Hz and ζ = 0.05, with the initial angle of attack equal to 10°.The initial flow field is determined by simulations that are carried out for an airfoil at a fixed 10°incidence until the time t = 0.60 s (corresponding to the dimensionless time t * = U ∞ ؒ t/c ≈ 48), when the time wise variations of the force coefficients become periodic.The airfoil is then released in the fluid and computations are performed in Fluid Structure Interaction.

Instantaneous velocity magnitude
For this airfoil, the static separation occurs around α = 10°.Therefore when the airfoil is at a fixed incidence α = 10°, the flow is separated on the extrados (Fig. 3a).In the rest of the domain, the flow-field is steady (Fig. 4a).As soon as the profile is free oscillating, some vortices detach from the airfoil surface and are spread of in the wake of the airfoil (Fig. 3b to 3d and 4b to 4d).This phenomenon persists even when the airfoil oscillations are low (Figs.3e and 4e).At the time t = 1.62 s from the beginning of the airfoil oscillations, the stationary flow is nearly established (Figs.3f and 4f).

Moment and force coefficients
The lift coefficient decreases quickly as soon as the profile is free oscillating, as far as reaching a zero lift coefficient.This result is reflected in the Fig. 5 which shows the variation of the torsion moment as a function of the angle of attack.

Airfoil oscillations
The time-wise variations of the airfoil position are depicted on Fig. 6.The curve show that the airfoil oscillations are attenuated and that the position tends to a zero incidence.This result was expected as the elastic axis is located at the aerodynamic centre.Moreover, the obtained results are merged with the analytic solution expressed by the relation (7).This is due to the high value of the airfoil natural frequency.

NACA 63 2 415 IN FLUTTER
Simulations of the classical flutter are performed for a NACA 63 2 415 airfoil similar to that of Svacek and al [10].The model has a chord length C = 0.30 m, a span of 0.50 m, with the following elastic parameters: The elastic axis EO is localized at 40% from the leading edge and the centre of gravity G is at 37% of the leading edge.Computations are started with initial displacements and zero initial vibration velocity.The two following configurations are considered: where y init and α init are the initial airfoil positions and U ∞ is the free stream velocity.
As previously, the initial flow-field is determined with computations carried out for a fixed airfoil at the above initial positions.The airfoil is then released in the fluid and the computations are performed in FSI.The obtained results are depicted according to the initial position and the speed index, a dimensionless parameter defined as with and , the natural frequency of the airfoil.The aim was to verify if the flutter point is attained or not.With the applied values of the airfoil elastic parameters, the free stream velocities U ∞ = 2, 26 and 45 m/s correspond to the speed indexes V * = 0.15, 1.99 and 3.45, respectively.

Results obtained with y init = +0.050 m and α init = 6°T
he contours of velocity magnitudes obtained with the different values of U ∞ are depicted on Figs. 7 to 9, where U max is the maximal velocity magnitude over the airfoil surface.When U ∞ = 2 m/s (the corresponding Reynolds number being Re = 4 × 10 4 ), it is found that the  vortex streets are more important in the near wake of the airfoil than with the higher free stream velocities.In addition, when the airfoil is at rest, the ratio increases slightly with U ∞ .With U ∞ =26 and 45 m/s, as soon as the airfoil is released in the fluid, the ratio decreases as reaching the value 1.30 at the time t ≈ 1s.However, when sensitive to the free stream velocity value.However, it is found that the airfoil responses are not influenced by the value of the fluid velocity inlet.The phenomenon of divergence is not observed (Fig. 11).Then the flutter point is not reached.Moreover, the results have similar behaviour to that of an analytic solution obtained with F Y =0 and M = 0.This is due to the airfoil elasticity parameters and dampening and to the low values of the speed index.In addition, our computed airfoil displacements are different to those of Ref. [10].The discrepancies are attributed to different initial airfoil displacements.

3.2.2.
Results obtained with y init = -0.050m, α init = 12°and v * = 1.99 Due to higher initial angle of incidence, the flow around the fixed airfoil is unsteady and it separates on the extrados; the ratio is more important.When the airfoil is free vibrating, the flow stabilizes and reattaches at the airfoil surface whereas the airfoil tends toward a position of balance with an incidence of 0°.
In comparison to the previous computations with y init = +0.050m, α init = 6°and the same speed index V * = 1.99, it is found that: (i) At the time t ≈ 1s from the beginning of the oscillations, the ratio (Fig. 12).(ii) The computed force coefficients are lower compared to the prior configuration with the same value of the free stream velocity (Fig. 13).(iii) The initial amplitudes of the airfoil pitch oscillations are higher; however, the airfoil vibrations are rapidly damped.The vertical oscillations are virtually opposite (Fig. 14).
Other simulations whose results are presented in a companion paper have been performed for an airfoil without dampening and with different value of the airfoil natural frequency.The influence of the free-stream velocity on the airfoil oscillations is then more important.The method was successfully applied for the computation of the viscous laminar flows around a free pitch oscillating airfoil and around an airfoil in combined pitch and vertical oscillations (classic Flutter).As wind turbines operate in turbulent environment, further computations will be performed with the of the incompressible Reynolds averaged Navier-Stokes equations.
This method can easily be extended to 3D computational domain to study aeroelastic stability of a parked wind turbine blade at high wind speeds.However, although the numerical tools performance are more and more increasing, a FSI approach is hard to apply for the aeroelasticity modelling of the entire wind turbine i.e. including the tower and the generator.Then future challenges lies in developing the Reduced Order Modelling (ROM) technique that is efficient to model larger systems.In particular, the Proper Orthogonal Decomposition (POD) technique, that performances have been confirmed in fluid mechanics, has been successfully tested in FSI applications [18,19].

Figure 2
Figure 2 Comparison of the discretisation schemes.

Figure 3
Figure 3 Contours of the velocity magnitude in the neighbourhood of the NACA 0012 airfoil.

1 .Figure 4
Figure 4 Contours of the velocity magnitude around the pitch oscillating NACA 0012 in the whole computational domain.

396 Numerical simulation of the flow around oscillating wind turbine airfoils 10 Figure 5
Figure 5 Variation of the moment with the pitch angle for the oscillating NACA 0012 airfoil.

Figure 6
Figure 6 Time wise variation of the angle of attack of the pitch oscillating NACA 0012 airfoil (t = 0s corresponds to the time at which the airfoil is released in the fluid).