Aeroelasticity of Wind Turbines Blades Using Numerical Simulation

With roller coaster traditional fuel prices and ever increasing energy demand, wind energy has known significant growth over the last years. To pave the way for higher efficiency and profitability of wind turbines, advances have been made in different aspects related to this technology. One of these has been the increasing size of wind turbines, thus rendering the wind blades gigantic, lighter and more flexible whilst reducing material requirements and cost. This trend towards gigantism increases risks of aeroelastic effects including dire phe‐ nomena like dynamic stall, divergence and flutter. These phenomena are the result of the combined effects of aerodynamic, inertial and elastic forces. In this chapter, we are present‐ ing a qualitative overview followed by analytical and numerical models of these phenomena and their impacts on wind turbine blades with special emphasize on Computational Fluid Dynamics (CFD) methods. As definition suggests, modeling of aeroelastic effects require the simultaneous analysis of aerodynamic solicitations of the wind flow over the blades, their dynamic behavior and the effects on the structure. Transient modeling of each of these char‐ acteristics including fluid-structure interaction requires high level computational capacities. The use of CFD codes in the preprocessing, solving and post processing of aeroelastic prob‐ lems is the most appropriate method to merge the theory with direct aeroelastic applications and achieve required accuracy. The conservation laws of fluid motion and boundary condi‐ tions used in aeroelastic modeling will be tackled from a CFD point of view. To do so, the chapter will focus on the application of finite volume methods to solve Navier-Stokes equa‐ tions with special attention to turbulence closure and boundary condition implementation. Three aeroelastic phenomena with direct application to wind turbine blades are then stud‐ ied using the proposed methods. First, dynamic stall will be used as case study to illustrate


Introduction
With roller coaster traditional fuel prices and ever increasing energy demand, wind energy has known significant growth over the last years.To pave the way for higher efficiency and profitability of wind turbines, advances have been made in different aspects related to this technology.One of these has been the increasing size of wind turbines, thus rendering the wind blades gigantic, lighter and more flexible whilst reducing material requirements and cost.This trend towards gigantism increases risks of aeroelastic effects including dire phenomena like dynamic stall, divergence and flutter.These phenomena are the result of the combined effects of aerodynamic, inertial and elastic forces.In this chapter, we are presenting a qualitative overview followed by analytical and numerical models of these phenomena and their impacts on wind turbine blades with special emphasize on Computational Fluid Dynamics (CFD) methods.As definition suggests, modeling of aeroelastic effects require the simultaneous analysis of aerodynamic solicitations of the wind flow over the blades, their dynamic behavior and the effects on the structure.Transient modeling of each of these characteristics including fluid-structure interaction requires high level computational capacities.The use of CFD codes in the preprocessing, solving and post processing of aeroelastic problems is the most appropriate method to merge the theory with direct aeroelastic applications and achieve required accuracy.The conservation laws of fluid motion and boundary conditions used in aeroelastic modeling will be tackled from a CFD point of view.To do so, the chapter will focus on the application of finite volume methods to solve Navier-Stokes equations with special attention to turbulence closure and boundary condition implementation.Three aeroelastic phenomena with direct application to wind turbine blades are then studied using the proposed methods.First, dynamic stall will be used as case study to illustrate the traditional methodology of CFD aeroelastic modeling: mathematical analysis of the phenomenon, choice of software, computational domain calibration, mesh optimization and turbulence and transition model validation.An S809 airfoil will be used to illustrate the phenomenon and the obtained results will be compared to experimental ones.The divergence will be then studied both analytically and numerically to emphasize CFD capacity to model such a complex phenomenon.To illustrate divergence and related study of eigenvalues, an experimental study conducted at NASA Langley will be analyzed and used for comparison with our numerical modeling.In addition to domain, mesh, turbulence and transition model calibration, this case will be used to illustrate fluid-structure interaction and the way it can be tackled in numerical models.Divergence analysis requires the modeling of flow parameters on one side and the inertial and structural behavior of the blade on the other side.These two models should be simultaneously solved and continuous exchange of data is essential as the fluid behavior affects the structure and vice-versa.This chapter will conclude with one of the most dangerous and destructive aeroelastic phenomena -the flutter.Analytical models and CFD tools are applied to model flutter and the results are validated with experiments.This example is used to illustrate the application of aeroelastic modeling to predictive control.The computational requirements for accurate aeroelastic modeling are so important that the calculation time is too large to be applied for real time predictive control.Hence, flutter will be used as an example to show how we can use CFD based offline results to build Laplacian based faster models that can be used for predictive control.The results of this model will be compared to experiments.

Characteristics of aeroelastic phenomena
Aeroelasticity refers to the science of the interaction between aerodynamic, inertial and elastic effects.Aeroelastic effects occur everywhere but are more or less critical.Any phenomenon that involves a structural response to a fluid action requires aeroelastic consideration.In many cases, when a large and flexible structure is submitted to a high intensity variable flow, the deformations can be very important and become dangerous.Most people are familiar with the "auto-destruction" of the Tacoma Bridge.This bridge, built in Washington State, USA, 1.9 km long, was one of the longest suspended bridges of its time.The bridge connecting the Tacoma Narrows channel collapsed in a dramatic way on Thursday, November 7,1940.With winds as high as 65-75 km/h, the oscillations increased as a result of fluidstructure interaction, the base of Aeroelasticity, until the bridge collapsed.Recorded videos of the event showed an initial torsional motion of the structure combined very turbulent winds.The superposition of these two effects, added to insufficient structural dumping, amplified the oscillations.Figure 1 below illustrates the visual response of a bridge subject to aeroelastic effects due to variable wind regimes.The simulation was performed using multiphysics simulation on ANSYS-CFX software.Some more details on similar aeroelastic modelling can be viewed from [1], [2], [3] and [4].In an attempt to increase power production and reduce material consumption, wind turbines' blades are becoming increasingly large yet, paradoxically, thinner and more flexible.The risk of occurrence of damaging aeroelastic effects increases significantly and justifies the efforts to better understand the phenomena and develop adequate design tools and mitigation techniques.Divergence and flutter on an airfoil will be used as introduction to aeroelastic phenomena.When a flexible structure is subject to a stationary flow, equilibrium is established between the aerodynamic and elastic forces (inertial effects are negligible due to static condition).However, when a certain critical speed is exceeded, this equilibrium is disrupted and destructive oscillations can occur.This is illustrated with Figure 2 where α is the angle of attack due to a torsional movement as a result of aerodynamic solicitations.If we consider an angle of attack sufficiently small such that cosα ≈ 1 and sinα≈ α, and writing the equilibrium of the moments, M, with respect to the centre of the rotational spring, we have: Where the lift L is: S, surface area of the profile, C l is the lift coefficient, M 0 is the moment coefficient.This leads to an angle of attack at equilibrium corresponding to: For a zero flow condition, the angle of attack α z , is such that: Divergence occurs when denominator in equation (3) becomes 0 and this corresponds to a dynamic pressure, q D expressed as: When velocity increases such that dynamic pressure q approaches critical dynamic pressure q D , the angle of attack dangerously increases until a critical failure value -divergence.This is solely a structural response due to increased aerodynamic solicitation due to fluid-structure interaction.This is an example of a static aeroelastic phenomenon as it involves no vibration of the airfoil.Flutter is an example of a dynamic aeroelastic phenomenon as it occurs when structure vibration interacts with fluid flow.It arises when structural damping becomes insufficient to damp aerodynamic induced vibrations.Flutter can appear on any flexible vibrating object submitted to a strong flow with positive retroaction between flow fluctuations and structural response.When the energy transferred to the blade by aerodynamic excitation becomes larger than the normal dynamic dissipation, the vibration amplitude increases dangerously.Flutter can be illustrated as a superposition of two structural modes -the angle of attack (pitch) torsional motion and the plunge motion which characterises the vertical flexion of the tip of the blade.Pitch is defined as a rotational movement of the profile with respect to its elastic center.As velocity increases, the frequencies of these oscillatory modes coalesce leading to flutter phenomenon.This may start with a rotation of the blade section (at t=0 s in Figure 3).The increased angle amplifies the lift such that the section undertakes an upward vertical motion.Simultaneously, the torsional rigidity of the structure recoils the profile to its zero-pitch condition (at t=T/4 in Figure 3).The flexion rigidity of the structure tends to retain the neutral position of the profile but the latter then tends to a negative angle of attack (at t = T/2 in Figure 3).Once again, the increased aerodynamic force imposes a downward vertical motion on the profile and the torsional rigidity of the latter tends to a zero angle of attack.The cycle ends when the profile retains a neutral position with a positive angle of attack.With time, the vertical movement tends to damp out whereas the rotational movement diverges.If freedom is given to the motion to repeat, the rotational forces will lead to blade failure.

Mathematical analytical models
Several examples of aeroelastic phenomena are described in the scientific literature.When it comes to aeroelastic effects related to wind turbines, three of the most common and dire ones are dynamic stall, aerodynamic divergence and flutter.In this section, we will provide a summarized definition of these aeroelastic phenomena with associated mathematical analytical models.Few references related to analytic developments of aeroelastic phenomena are [10], [11] and [12].

Dynamic stall
In fluid dynamics, the stall is a lift coefficient reduction generated by flow separation on an airfoil as the angle of attack increases.Dynamic stall is a nonlinear unsteady aerodynamic effect that occurs when there is rapid change in the angle of attack that leads vortex shedding to travel from the leading edge backward along the airfoil [14].The analytical development of equations characterizing stall will be performed using illustrations of Figures 4 and  5.The lift per unit length, expressed as L is given by: Where • ρ is the air density • c is the chord length of the airfoil We will, first, present static stall as described in [13].During stationary flow conditions, no flow separation occurs and the lift, L, acts approximately at the quarter cord distance from the leading edge at the pressure (aerodynamic) centre.For small values of α, L varies linearly with α.Stall happens at a critical angle of attack whereby the lift reaches a maximum value and flow separation on the suction side occurs.For unsteady conditions, a delay exists prior to reaching stability and is an essential condition for building analytical stall models [15].In such case, we can observe a smaller lift for an increasing angle of attack (AoA) and a larger one for decreasing AoA when compared with a virtually static condition.In a flow separation condition, we can observe a more significant delay which expresses itself with harmonic movements in the flow which affects the aerodynamic stall phenomenon.Figure 5, an excerpt from [16], shows that for harmonic variations of the AoA between 0 o and 15 o, , the onset of stall is delayed and the lift is considerably smaller for the decreasing AoA trend than for the ascending one.Hence, as expressed by Larsen et al. [17], dynamic stall includes harmonic motion separated flows, including formation of vortices in the vicinity of the leading edge and their transport to the trailing edge along the airfoil.Figures 6-9, which are excerpts from [18], illustrate these phenomena.Stall phenomenon is strongly non-linear such that a clear cut analytical solution model is impossible to achieve.This complex phenomenon requires consideration of numerous parameters, study of flow transport, boundary layer analysis (shape factor and thickness), vortex creation and shedding as well as friction coefficient consideration in the boundary layer.The latter helps in the evaluation of separation at the leading edge and is important for aeroelastic consideration.The proper modelling of transition from laminar to turbulent flow is also essential for accurate prediction of stall parameters.

Divergence
We consider a simplified aeroelastic system of the NACA0012 profile to better understand the divergence phenomenon and derive the analytical equation for the divergence speed.Figure 10 illustrates a simplified aeroelastic system, the rigid NACA0012 profile mounted on a torsional spring attached to a wind tunnel wall.The airflow over the airfoil is from left to right.The main interest in using this model is the rotation of the airfoil (and consequent twisting of the spring), α, as a function of airspeed.If the spring were very stiff and/or airspeed very slow, the rotation would be rather small; however, for flexible spring and/or high flow velocities, the rotation may twist the spring beyond its ultimate strength and lead to structural failure.The airspeed at which the elastic twist increases rapidly to the point of failure is called the divergence airspeed, U D .. This phenomenon, being highly dangerous and prejudicial for wind blades, makes the accurate calculation of U D very important.For such, we define C as the chord length and S as the rigid surface.The increase in the angle of attack is controlled by a spring of linear rotation attached to the elastic axis localized at a distance e behind the aerodynamic centre.The total angle of attack measured with respect to a zero lift position equals the sum of the initial angle α r and an angle due to the elastic deformation θ, known as the elastic twist angle.
The elastic twist angle is proportional to the moment at the elastic axis: where C θθ is the flexibility coefficient of the spring.The total aerodynamic moment with re- spect to the elastic axis is given by: where • S is the rigid surface area of the blade section The lift coefficient is related to the angle of attack measured with respect to a zero lift condition as follows: Here ∂ C l ∂ α represents the slope of the lift curve.The elastic twist angle θ, can be obtained by simple mathematical manipulations of the three previous equations: Regrouping θ : This leads to: We can note that for a given value of the dynamic pressure q, the denominator tends to zero such that the elastic twist angle will then tend to infinity.This condition is referred to as aerodynamic divergence.When the denominator tends to zero: The dynamic pressure is given by: Thus, we come up with: Hence, the divergence velocity can be expressed as: To calculate the theoretical value of the divergence velocity, certain parameters need to be found.These are C θθ , which is specific to the modeled spring, S and e being inherent to the airfoil, ρ depends upon the used fluid and ∂ C l ∂ α depends both on the shape of the airfoil and flow conditions [23].We note that as divergence velocity is approached, the elastic twist angle will increase in a very significant manner towards infinity [24].However, computing is finite and cannot model infinite parameters.Therefore, the value of the analytical elastic twist angle is compared with the value found by the coupling.In the case wherein the elastic twist angle introduces no further aerodynamic solicitations, by introducing α = α r , and re- solving for the elastic twist angle, we have: Hence: which leads to: Hence, we note that the theoretical elastic twist angle depends on the divergence speed and the elastic twist angle calculated whilst considering that it triggers no supplementary aerodynamic solicitation.The latter is calculated by solving for the moment applied on the profile at the elastic axis (T) during trials in steady mode.These trials are conducted using the k-ω SST intermittency transitional turbulence model with a 0.94 intermittency value [25].In this section, we will present only the expression used to calculate the divergence speed.The development of this expression and the analytical calculation of a numerical value of the divergence speed are detailed in [28]: Detailed eigen values and eigenvectors analysis related to divergence phenomenon is presented in [27].

Aerodynamic flutter
As previously mentioned, flutter is caused by the superposition of two structural modespitch and plunge.The pitch mode is described by a rotational movement around the elastic centre of the airfoil whereas the plunge mode is a vertical up and down motion at the blade tip.Theodorsen [16][17][18] developed a method to analyze aeroelastic stability.The technique is described by equations ( 61) and (62).α is the angle of attack (AoA), α 0 is the static AoA, C(k) is the Theodorsen complex valued function, h the plunge height, L is the lift vector positioned at 0.25 of the chord length, M is the pitching moment about the elastic axis, U is the free velocity, ω is the angular velocity and a, b, d 1 and d 2 are geometrical quantities as shown in Figure 11.
Theodorsen's equation can be rewritten in a form that can be used and analyzed in Matlab Simulink as follows:

Flutter movement
The occurrence of the flutter has been illustrated in Section 2 (Figure 3).To better understand this complex phenomenon, we describe flutter as follows: aerodynamic forces excite the mass -spring system illustrated in Figure 12.The plunge spring represents the flexion rigidity of the structure whereas the rotation spring represents the rotation rigidity.

Flutter equations
The flutter equations originate in the relation between the generalized coordinates and the angle of attack of the model that can be written as: The Lagrangian form equations are constructed for the mechanical system.The first one corresponds to the vertical displacement z and the other is for the angle of attack α: Numerical solution of these equations requires expressing F z and M o as polynomials of α.C z = -0.0000983α 3 -0.0003562α 2 + 0.1312α + 0.4162 C m0 = -0.00006375α 3 + 0.00149α 2 -0.001185 α -0.9312These equations will be used in the modeling of a lumped representation of flutter presented in the last section of this chapter.

Computational fluid dynamics (CFD) methods in aeroelastic modeling
Aeroelastic modeling of wind blades require complex representation of both fluid flows, including turbulence, and structural response.Fluid mechanics aims at modeling fluid flow and its effects.When the geometry gets complex (flow becomes unsteady with turbulence intensity increasing), it is impossible to solve analytically the flow equations.With the advent of high efficiency computers, and improvement in numerical techniques, computational fluid dynamics (CFD), which is the use of numerical techniques on a computer to resolve transport, momentum and energy equations of a fluid flow has become more and more popular and the accuracy of the technique has been an a constant upgrading trend.Aeroelastic modeling of wind blades includes fluid-structure interaction and is, in fact, a science which studies the interaction between elastic, inertial and aerodynamic forces.The aeroelastic analysis is based on modeling using ANSYS and CFX software.CFX uses a finite volume method to calculate the aerodynamic solicitations which are transmitted to the structural module of ANSYS.Within CFX, several parameters need to be defined such as the turbulence model, the reduced frequency, the solver type, etc. and the assumptions and limitations of each model need to be well understood in order to validate the quality and pertinence of any aeroelastic model.These calibration considerations will be illustrated in the stall modeling section as an example.

Dynamic stall
In this section of the chapter, we will illustrate aeroelastic modeling of dynamic stall on a S809 airfoil for a wind blade.The aim of this section, apart from illustrating this aeroelastic phenomenon is to emphasize on the need of parameter calibration (domain size, mesh size, turbulence and transition model) in CFD analysis.The different behaviour of lift as the AoA increases or decreases leads to significant hysteresis in the air loads and reduced aerodynamic damping, particularly in torsion.This can cause torsional aeroelastic instabilities on the blades.Therefore, the consideration of dynamic stall is important to predict the unsteady blade loads and, also, to define the operational boundaries of a wind turbine.In all the following examples, the CFD based aeroelastic models are run on the commercial ANSYS-CFX software.CFX, the fluid module of the software, models all the aerodynamic parameters of the wind flow.ANSYS structural module defines all the inertial and structural parameters of the airfoil and calculates the response and stresses on the structure according to given solicitations.The MFX module allows fluid-structure modelling, i.e., the results of the aerodynamic model are imported as solicitations in the structural module.The continuous exchange of information allows a multi-physics model that, at all time, computes the action of the fluid on the structure and the corresponding impact of the airfoil motion on the fluid flow.

Model and convergence studies 4.1.1.1. Model and experimental results
In an attempt to calibrate the domain size, mesh size, turbulence model and transition model, an S809 profile, designed by NREL, was used.This airfoil has been chosen as experimental results and results from other sources are available for comparison.The experimental results have been obtained at the Low Speed Laboratory of the Delft University [31] and at the Aeronautical and Astronautical Research Laboratory of the Ohio State University [32].The first work [31], performed by Somers, used a 0.6 meters chord model at Reynolds numbers of 1 to 3 million and provides the characteristics of the S809 profile for angles of incidence from -20 0 to 20 0 .The second study [32], realized by Ramsey, gives the characteristics of the airfoil for angles of incidence ranging from -20 0 to 40 0 .The experiments were conducted on a 0.457 meters chord lenght for Reynolds numbers of 0.75 to 1.5 million.Moreover, this study provides experimental results for the study of the dynamic stall for incidence angles of (8 0 , 14 0 and 20 0 ) oscillating at (±5.5 0 and ± 10 0 ) at different frequencies for Reynolds numbers between 0.75 and 1.4 million.

Convergence studies
In this section, we will focus on the definition of a calculation domain and an adapted mesh for the flow modelling around the mentioned airfoil.This research is realized by the study of the influence of the distance between the boundaries and the airfoil, the influence of the size of the chord for the same Reynolds number and finally, the influence of the number of elements in the mesh and computational time.

Computational domain
The computational domain is defined by a semi-disc of radius I 1 ×c around the airfoil and two rectangles in the wake, of length I 2 ×c.This was inspired from works conducted by Bhaskaran presented in Fluent tutorial.As the objective of this study was to observe how the distance between the domain boundary and the airfoil affects the results, only I 1 and I 2 were varied with other values constant.As these two parameters will vary, the number of elements will also vary.To define the optimum calculation domain, we created different domains linked to a preliminary arbitrary one by a homothetic transformation with respect to the centre G and a factor b. Figure 14 gives us an idea of the different parameters and the outline of the computational domain whereas table 1 presents a comparison of the different meshed domains.The drag coefficient diminishes as the homothetic factor increases but tends to stabilize.This stabilization is faster for low angles of attack (AoA) and seems to be delayed for larger homothetic factors and increasing AoA.The trend for the lift coefficients as a function of the homothetic factor is quite similar for the different angles of attack except for an angle of 8.2 0 .The evolution of the coefficients towards stabilization illustrate an important physical phenomenon: the further are the boundary limits from the airfoil, this allows more space for the turbulence in the wake to damp before reaching the boundary conditions imposed on the boundaries.Finally, a domain having a radius of semi disc 5.7125 m, length of rectangle 9,597 m and width 4.799 m was used.

Meshing
Unstructured meshes were used and were realized using the CFX-Mesh.These meshes are defined by the different values in table 2. We kept the previously mentioned domain.

Table 2. Mesh parameters
Figure 16 gives us an appreciation of the mesh we have used of in our simulations: Several trials were performed with different values of the parameters describing the mesh in order to have the best possible mesh.Lift and drag coefficient distributions have been computed according to different AoA for a given Reynolds number and the results were compared with experiments.The mesh option that provided results which fitted the best with the experimental results was used.The final parameters of the mesh were 66772 nodes and 48016 elements.[24] presents these comparisons which lead us to the following conclusions: the k-ω SST model is the only one to have a relatively good prediction of the large separated flows for high angles of attack.So, the transport of the turbulent shear stress really improves the simulation results.The consideration of the transport of the turbulent shear stress is the main asset of the k-ω SST model.However, probably a laminar-turbulent transition added to the model will help it to better predict the lift coefficient between 6° and 10°, and to have a better prediction of the pressure coefficient along the airfoil for 20°.This assumption will be studied in the next section where the relative performance of adding a particular transitional model is studied.

Transition model
ANSYS-CFX proposes in the advanced turbulence control options several transitional models namely: the fully turbulent k-ω SST model, the k-ω SST intermittency model, the gamma theta model and the gamma model.As the gamma theta model uses two parameters to define the onset of turbulence, referring to [33], we have only assessed the relative performance of the first three transitional models.The optimum value of the intermittency parameter was evaluated.A transient flow analysis was conducted on the S809 airfoil for the same Reynolds number at different AoA and using three different values of the intermittency parameter: 0.92, 0.94 and 0.96. Figure 17 illustrates the drag and lift coefficients obtained for these different models at different AoA as compared to DUT and OSU experimental data.We note that for the drag coefficient, the computed results are quite similar and only dif-fer in transient mode exhibiting different oscillations.For large intermittency values, the oscillations are larger.Figure 18 shows that for the lift coefficients, the results from CFX differ from 8.2 0 .For the linear growth zone, the different results are close to each other.The difference starts to appear near maximum lift.The k-ω SST intermittency model with γ=0.92, under predicts the lift coefficients as compared with the experimental results.The results with γ=0.94 predicts virtually identical results as compared to the OSU results.The model with γ= 0.96 predicts results that are sandwiched between the two experimental ones.Analysis of the two figures brings us to the conclusion that the model with γ=0.94 provides results very close to the DUT results.Therefore, we will compare the intermittency model with γ=0.94 with the other transitional models.Figure 18 shows that the drag coefficients for the three models are very close until 18 0 after which the results become clearly distinguishable.As from 20 0 , the γ-θ model over predicts the experimental values whereas such phenomena appear only after 22.1 0 for the two other models.For the lift coefficients, Figure 18 shows that the k-ω SST intermittency models provide results closest to the experimental values for angles smaller than 14 0 .The k-ω SST model under predicts the lift coefficients for angles ranging from 6 0 to 14 0 .For angles exceeding 20 0 , the intermittency model does not provide good results.Hence, we conclude that the transitional model helps in obtaining better results for AoA smaller than 14 0 .However, for AoA greater than 20 0 , a purely turbulent model needs to be used.

Results
In order to validate the quality of stall results, the latter are compared with OSU experimental values and with Leishman-Beddoes model.Moreover, modelling of aeroelastic phenomena is computationally very demanding such that we have opted for an oscillation of 5.5 0 around 8 0 , 14 0 and 20 0 for a reduced frequency of k = ωc 2U ∞ = 0.026, where c is the length of the chord of the airfoil and U ∞ is the unperturbed flow velocity.From a structural point of view, the 0.457 m length profile will be submitted to an oscillation about an axis located at 25% of the chord.The results which follow illustrates the quality of our aeroelastic stall modelling at three different angles, all with a variation of 5.5 sin(w)*t.α = 8 0 ± 5.5sin(ωt) 0 Figure 19 illustrates the evolution of the aerodynamic coefficients with oscillation of the AoA around 8 0 with amplitude 5.5 0 and a reduced frequency of 0.026.When the angle is less than 14 0 , the transitional k-ω SST intermittency model was used.For the drag coefficient, the results are very close to experimental ones and limited hysteresis appears.As for the lift coefficient, we note that the «k-ω SST intermittency» transitional turbulent model underestimates the hysteresis phenomenon.Furthermore, this model provides results with inferior values as compared to experimental ones for increasing angle of attack and superior values for decreasing angle of attack.We, also, notice that the onset of the stall phenomenon is earlier for the «k-ω SST intermittency» model.α = 14 0 ± 5.5sin(ωt) 0 Figure 20 illustrates the evolution of the aerodynamic coefficients with oscillation of the AoA around 14 0 with amplitude of 5.5 0 and a reduced frequency of 0.026.As the angle exceeds 14 0 , the purely turbulent k-ω SST model was used.
We note that before 17°, the model overestimates the drag coefficient, both for increasing and decreasing AoA.For angles exceeding 17°, the model approaches experimental results.As for the drag coefficient, the model provides better results for the lift coefficient when the angle of attack exceeds 17°.Furthermore, we note that the model underestimates the lift coefficient for both increasing and decreasing angles of attack.Moreover, the predicted lift coefficients are closer to experimental results for AoA less than 13°. Figure 21 illustrates the evolution of the aerodynamic coefficients with oscillation of the Ao Aaround 20 0 with amplitude of 5.5 0 and reduced frequency of 0.026.As the angle exceeds 14 0 , the purely turbulent k-ω SST model was used.For the drag coefficients, the results provided by the k-ω SST model are quite close to the experimental results but predict premature stall for increasing AoA and reattachment for decreasing AoA.Furthermore, we note oscillations of the drag coefficient for the experimental data showing higher levels of turbulence.The k-ω SST model under predicts the value of the lift coefficient for increasing AoA.Furthermore, due to vortex shedding, we note oscillations occurring in the results.Finally, the model prematurely predicts stall compared to the experiments.

Conclusions on stallmodelling
In this section, we presented an example of aeroelastic phenomenon, the dynamic stall.We have seen through this section the different steps to build and validate the model.Better results are obtained for low AoA but as turbulence intensity gets very large, the results diverge from experimental values or show oscillatory behaviour.We note that, though, the CFD models show better results than the relatively simple indicial methods found in literature, refinements should be brought to the models.Moreover, this study allows us to appreciate the complexity of fluid structure interaction and the calibration work required upstream.It should be emphasized that the coupling were limited both by the structural and aerodynamic models and refinements and better understanding of all the parameters that can help achieve better results.This study allows us to have a very good evaluation of the different turbulence models offered by CFX and their relative performances.

Aerodynamic divergence
In this section we will illustrate the different steps in modeling another aeroelastic phenomenon, the divergence and whilst using this example to lay emphasis on the ability of CFX-ANSYS software to solve fluid-structure interaction problems.As from the 1980s, national and international standards concerning wind turbine design have been enforced.With the refinement and growth of the state of knowledge the "Regulation for the Certification of Wind Energy Conversion Systems" was published in 1993 and further amended and refined in 1994 and 1998.Other standards aiming at improving security for wind turbines have been published over the years.To abide to such standards, modelling of the aeroelastic phenomena is important to correctly calibrate the damping parameters and the operation conditions.For instance, Nweland [34] makes a proper and complete analysis of the critical divergence velocity and frequencies.These studies allow operating the machines in secure zones and avoid divergence to occur.Such studies' importance is not only restrained to divergence but also apply to other general dynamic response cases of wind turbines.Wind fluctuations at frequencies close to the first flapwise mode blade natural frequency excite resonant blade oscillations and result in additional, inertial loadings over and above the quasi-static loads that would be experienced by a completely rigid blade.Knowledge of the domain of such frequencies allows us to correctly design and operate the machines within IEC and other norms.We here present a case where stall can be avoided by proper knowledge of its parameters and imposing specific damping.As the oscillations result from fluctuations of the wind speed about the mean value, the standard deviation of resonant tip displacement can be expressed in terms of the wind turbulence intensity and the normalized power spectral density at the resonant frequency, R u (n1) [34]: where: • x 1 is the first mode component of the steady tip displacement.
• U is the mean velocity (usually averaged over 10 minutes) • δ is the logarithmic decrement damping • K sx (n1) is the size of the reduction factor which is present due to lack of correlation of wind along the blade at the relevant frequency.
It is clear from equation ( 32) that a key determinant of resonant tip response is the value of damping present.If we consider for instance a vibrating blade flat in the wind, the fluctuating aerodynamic force acting on it per unit length is given by: μ 1 (r) is the first mode shape and m 1 is the generalized mass given by: Here, ω 1 is the first mode natural frequency given in radian per second.The logarithmic decrement is obtained by multiplying the damping ratio by 2 π.To properly estimate operating conditions and damping parameters, knowledge of the vibration frequencies and shape modes are important.The need to know these limits is again justified by the fact that when maximum lift is theoretically achieved toward maximum power when stall and other aeroelastic phenomena are also approached.

ANSYS-CFX coupling
To achieve the fluid structure coupling study, we make use of the ANSYS multi-domain (MFX).This module was primarily developed for fluid-structure interaction studies.On one side, the structural part is solved using ANSYS Multiphysics and on the other side, the fluid part is solved using CFX.The study needs to be conducted on a 3D geometry.If the geometries used by ANSYS and CFX need to have common surfaces (interfaces), the meshes of these surfaces can be different.The ANSYS code acts as the master code and reads all the multi-domain commands.It recuperates the interface meshes of the CFX code, creates the mapping and communicates the parameters that control the timescale and coupling loops to the CFX code.The ANSYS generated mapping interpolates the solicitations between the different meshes on each side of the coupling.Each solver realizes a sequence of multi-domain, time marching and coupling iterations between each time steps.For each iteration, each solver recuperates its required solicitation from the other domain and then solves it in the physical domain.Each element of interface is initially divided into n interpolation faces (IP) where n is the number of nodes on that face.The 3D IP faces are transformed into 2D poly-gons.We, then, create the intersection between these polygons, on one hand, the solver diffusing solicitations and on the other hand, the solver receiving the solicitations.This intersection creates a large number of surfaces called control surfaces as illustrated in Figure 22.These surfaces are used in order to transfer the solicitation between the structural and fluid domains.We can make use of different types of resolutions, either using a simultaneous scheme or using a sequential scheme, in which case we need to choose which domain to solve first.For lightly coupled domains, CFX literature recommends the use of the simultaneous scheme.As for our case, the domains are strongly coupled and for such reasons, we make use of the sequential scheme.This scheme has as advantage to ensure that the most recent result or solicitation of a domain solver is applied to the other solver.In most simulations; the physics of one domain imposes the requirements of the other domain.Hence, it is essential to ade-quately choose the code to solve first in the sequential scheme.In the case of the divergence, it is the fluid that imposes the solicitations on the solid such that the CFX code will be the first to be solved followed by the ANSYS code.The ANSYS workbench flow-charts that illustrates such interaction is illustrated in Figure 24 below:

Overview of the experimental results
An aeroelastic experiment was conducted at the Duke University Engineering wind tunnel facility [35].The goals of this test were to validate the analytical calculations of non-critical mode characteristics and to explicitly examine the aerodynamic divergence phenomenon.

Configuration description
The divergence assessment testbed (dat) wind tunnel model consists of a typical section airfoil with a flexible mount system providing a single degree of freedom structural dynamic mode.The only structural dynamic mode of this model is torsional rotation, or angle of attack.The airfoil section is a NACA 0012 with an 8-inch chord and a span of 21 inches.The ratio of the trailing edge mass to the total mass is 0.01.This spans the entire test section from the floor to ceiling.The structural dynamic parameters for this model are illustrated in table 3: However, some parameters were unavailable in [35] such that an iterative design process was used to build the model in ANSYS.Using parameters specified in [35], a preliminary model was built and its natural frequencies verified using ANSYS.The model was successively modified until a model as close as possible to the model in the experiment was obtained.The aims of the studies conducted in [35] were to: 1) find the divergence dynamic pressure;2)examine the modal characteristics of non-critical modes, both sub-critically and at the divergence condition; 3) examine the eigenvector behaviour.Heeg [35] obtained several interesting results among which the following graphic showing the variation of the angle of attack with time.The aim of our simulations was to determine how the numerical AN-SYS-CFX model will compare with experiments.The test was conducted by setting as close as possible to zero the rigid angle of attack, α 0 , for a zero airspeed.The divergence dynamic pressure was determined by gradually increasing the velocity and measuring the system response until it became unstable.The dynamic pressure was being slowly increased until the angle of attack increased dramatically and suddenly.This was declared as the divergence dynamic pressure, 5.1 psf (244 N/m2).The time history shows that the model oscillates around a new angle of attack position, which is not at the hard stop of the spring.It is speculated that the airfoil has reached an angle of attack where flow has separated and stall has occurred [35].

The ANSYS-CFX model
The model used in the experiment was simulated using a reduced span-wise numerical domain (quasi 2D).The span of the airfoil was reduced 262.5 times, from 21 inches to 0.08 inches or 2.032 mm, while the chord of the airfoil was maintained at 8 inch or 203.2 mm.We used a cylinder to simulate the torsion spring used in the experiment.

Results
In [23], the authors have derived the analytical mathematical equation to calculate the divergence velocity, U D .The expression was: (37) In order to calculate the theoretical value of the divergence velocity, certain parameters need to be found first.These are C θθ , which is specific to the modeled spring, S being inherent to the profile, e, which depends both on the profile (elastic axis) and on the aerodynamic model, ρ, which is dependent upon the used fluid and ∂ C l ∂ α which depends both on the shape of the profile but, also, on the turbulent model [23].We note that, as divergence velocity is approached, the elastic twist angle will increase in a very significant manner and tend to infinity [24].However, numerical values are finite and cannot model infinite parameters.We will, therefore, formulate the value of the analytical elastic twist angle in order to compare it with the value found by the coupling.In the case wherein the elastic twist angle introduces no further aerodynamic solicitations, by introducing α = α r , and resolving for the elastic twist angle, we have: Algebraic manipulations of the expressions lead us to the following formulation: This leads to: (40)

Advances in Wind Power
Hence, we can note that the theoretical elastic twist angle depends on the divergence speed and the elastic twist angle calculated whilst considering that it triggers no supplementary aerodynamic solicitation.To calculate the latter, we will solve for the moment applied on the profile at the elastic axis (T) during trials in steady mode.These trials are conducted using the k-ω SST intermittency transitional turbulence model with a 0.94 intermittency value [24].To model the flexibility coefficient of the rotational spring C θθ , used in the NASA experiments we used a cylinder as a torsion spring.The constant of the spring used in the experiment is K α = 5.8262 N•m/rad and since we used a reduced model, with an span 262.5 times smaller than the original, the dimensions and properties of the cylinder are such that: and the flexibility coefficient is: The slope of the lift ∂ C l ∂ α , can be calculated for an angle α = 5 0 in the following way: We have calculated the lift coefficient at 4.0 0 and 6.0 0 such that: The distance e, between the elastic axis and the aerodynamic centre for the model is 0.375•b.The rigid area is calculated to be S, being the product of the chord and the span and is calculated as follows: Hence the divergence velocity is calculated as: The theoretical divergence speed given in Table 4 of the NASA experiment [35] is 19.15 m/s.This slight difference is due to the value of slope of the lift profile ∂ C l ∂ α taken into consideration, which in the NASA work was 2π, or 6.283 rad -1 , whereas we used a value of 6.532 rad -1 .Furthermore, a difference between our calculated speed and that presented in [35] might also be explained by the size of the used tunnel and the possible wall turbulence interaction.Furthermore, using the model, domain and mesh parameters detailed in the previous sections of this article, divergence was modelled as follows: the airfoil used in [35] was fixed and exempted from all rotational degrees of liberty and subjected to a constant flow of velocity 15 m s -1 .Suddenly, the fixing is removed and the constant flow can be then compared to a shock wave on the profile.The profile then oscillates with damped amplitude due to the aerodynamic damping imposed.Figure 27 illustrates the response portrayed by AN-SYS-CFX software.We can extract the amplitude and frequency of oscillation of around 8 Hz which is close to the 7.9 Hz frequency presented in [35].

Aerodynamic flutter
In this section, we illustrate a CFD approach of modeling the most complex and the most dangerous type of aeroelastic phenomenon to which wind turbine blades are subjected.While illustrating stall phenomenon, we calibrated the CFD parameters for aeroelastic modeling.In the divergence section, the example was used to reinforce the notion of multiphysics modeling, more precisely, emphasis was laid on fluid structure interaction modeling within ANSYS-CFX MFX.Flutter example will be used to illustrate the importance of using lumped method.

Computational requirement and Lumped model
Aeroelastic modeling requires enormous computational capacity.The most recent quad core 16 GB processor takes some 216 hours to simulate flutter on a small scale model and that for a 12 second real time frame.The aim of simulating and predicting aeroelastic effects on wind blades has as primary purpose to apply predictive control.However, with such enormous computational time, this is impossible.The need for simplified lumped (2D Matlab based) models is important.The CFD model is ran preliminarily and the lumped model is built according to simulated scenarios.In this section we will illustrate flutter modeling both from a CFD and lumped method point of view.

Matlab-Simulink and Ansys-CFX tools
For flutter modelling, again, ANSYS-CFX model was used to simulate the complex fluidstructure interaction.However, due to excessively important computational time that rendered the potential of using the predictive results for the application of mitigation control impossible, the results of the CFD model was used to build a less time demanding lumped model based on Simulink.Reference [36] describes the Matlab included tool Simulink as an environment for multi-domain simulation and Model-Based Design for dynamic and embedded systems.It provides an interactive graphical environment and a customizable set of block libraries that let you design, simulate, implement, and test a variety of time-varying systems.For the flutter modelling project, the aerospace blockset of Simulink has been used.The Aerospace Toolbox product provides tools like reference standards, environment models, and aerospace analysis pre-programmed tools as well as aerodynamic coefficient importing options.Among others, the wind library has been used to calculate wind shears and Dryden and Von Karman turbulence.The Von Karman Wind Turbulence model uses the Von Karman spectral representation to add turbulence to the aerospace model through preestablished filters.Turbulence is represented in this blockset as a stochastic process defined by velocity spectra.For a blade in an airspeed V, through a frozen turbulence field, with a spatial frequency of Ω radians per meter, the circular speed ω is calculated by multiplying V by Ω.For the longitudinal speed, the turbulence spectrum is defined as follows: Here,L ω represents the turbulence scale length and σis the turbulence intensity.The corresponding transfer function used in Simulink is: For the lateral speed, the turbulence spectrum is defined as: and the corresponding transfer function can be expressed as : We notice that at the beginning there is a non-established instability, followed by a recurrent oscillation.The peak to peak distance corresponds to around 2.5 seconds that is a frequency of 0.4 Hz.The oscillation can be defined approximately by an amplitude of 0 0 ± 17 0 . .The same modelling was performed using the Simulink model and the result for the AoA variation and the plunge displacement is shown below: We note that for the AoA variation, the aerospace blockset based model provides very similar results with Heeg's results [35].The amplitude is, also, around 0 0 ± 17 0 and the frequency is 0.45 Hz.Furthermore, we notice that the variation is very similar.We can conclude that the aerospace model does represent the flutter in a proper manner.It is important to note that this is a special type of flutter.The frequency of the beat is zero and, hence, represents divergence of "zero frequency flutter".Using Simulink, we will vary the angular velocity of the blade until the eigenmode tends to a negative damping coefficient.The damping coefficient, ζ is obtained as: ζ = c 2mω , ω is measured as the Laplace integral in Simulink, c is the viscous damping and ω= k m .Figure 30 illustrates the results obtained for the variation of the damping coefficient against rotational speed and flutter frequency against rotor speed.We can note that as the rotation speed increases, the damping becomes negative, such that the aerodynamic instability which contributes to an oscillation of the airfoil is amplified.We also notice that the frequency diminishes and becomes closer to the natural frequency of the system.This explains the reason for which flutter is usually very similar to resonance as it occurs due to a coalescing of dynamic modes close to the natural vibrating mode of the system.We present here the results obtained for the same case study using ANSYS-CFX.The frequency of the movement using Matlab is 6.5 Hz while that using the ANSYS-CFX model is 6.325 Hz compared with the experimental value of 7.1Hz [35].Furthermore, the amplitudes of vibration are very close as well as the trend of the oscillations.For the points identified as 1, 2 and 3 on the flutter illustration, we illustrate the relevant flow over the airfoil.The maximum air speed at moment noted 1 is 26.95 m/s.We note such a velocity difference over the airfoil that an anticlockwise moment will be created which will cause an increase in the angle of attack.Since the velocity, hence, pressure difference, is very large, we note from the flutter curve, that we have an overshoot.The velocity profile at moment 2, i.e., at 1.88822s shows a similar velocity disparity, but of lower intensity.This is visible as a reduction in the gradient of the flutter curve as the moment on the airfoil is reduced.Finally at moment 3, we note that the velocity profile is, more or less, symmetric over the airfoil such that the moment is momentarily zero.This corresponds to a maximum stationary point on the flutter curve.After this point, the velocity disparity will change position such that angle of attack will again increase and the flutter oscillation trend maintained, but in opposite direction.This cyclic condition repeats and intensifies as we have previously proved that the damping coefficient tends to a negative value.

Figure 1 .
Figure 1.Aeroelastic response of a bridge

Figure 4 .
Figure 4. Illustration of an airfoil used for analytical development of stall related equations

Figure 5 .
Figure 5. Lift coefficient under static and dynamic stall conditions (dashed line for steady conditions, plain line for unsteady conditions)

Figure 6 .
Figure 6.Aerodynamic stall mechanism-Onset of separation on the leading edge

Figure 10 .
Figure 10.Simplified aeroelastic model to illustrate divergence phenomenon

Figure 12 .
Figure 12.Illustration of both pitch and plunge

1 2
ρLSV 2 C m0 (α) for S being the surface of the blade, C z , the lift coefficient, C m0 being the pitch coefficient, F z being the lift, M o , the pitch moment.C z and C m values are extracted from NACA 4412 data.Third degree interpolations for C z and C m with respect to the AoA are given below:

Figure 14 .
Figure 14.Shape of the calculation domainFigure 15 below respectively illustrates the drag and lift coefficients as a function of the homothetic factor b for different angles of attack.

Figure 16 .
Figure 16.Unstructured mesh along airfoil, boundary layer at leading edge and boundary layer at trailing edge

4. 1 . 1 . 5 .
Turbulence model calibrationCFX proposes several turbulence models for resolution of flow over airfoils.Scientific literature makes it clear that different turbulence models perform differently in different applications.CFX documentation recommends the use of one of three models for such kind of applications, namely the k-ω model, the k-ω BSL model and the k-ω SST model.The Wilcox k-ω model is reputed to be more accurate than k-ε model near wall layers.It has been successfully used for flows with moderate adverse pressure gradients, but does not succeed well for separated flows.The k-ω BSL model (Baseline) combines the advantages of the Wilcox k-ω model and the k-ε model but does not correctly predict the separation flow for smooth surfaces.The k-ω SST model accounts for the transport of the turbulent shear stress and overcomes the problems of k-ω BSL model.To evaluate the best turbulence model for our simulations, steady flow analyses at Reynolds number of 1 million were conducted on the S809 airfoil using the defined domain and mesh.The different values of lift and drag obtained with the different models were compared with the experimental OSU and DUT results.D'Hamonville et al.

Figure 17 .
Figure 17.Drag and lift coefficients for different AoA using different intermittency valuesFigure 18 illustrates the drag and lift coefficients for different AoA using different transitional models.

Figure 18 .
Figure 18.Drag and lift coefficients for different AoA using different transitional models

Figure 19 .
Figure 19.Drag and lift coefficients vs. angle of attack for stall modelling

Figure 21 .
Figure 21.Drag and lift coefficients vs. angle of attack for stall modelling where x ˙ is the blade flatwise velocity, C d is the drag coefficient and C(r) is the local blade chord.Hence the aerodynamic damping per unit length, C a ^(r) = ρU -C d C(r) and the first aerodynamic damping mode is:

Figure 22 .
Figure 22.Transfer Surfaces The respective MFX simultaneous and sequential resolution schemes are presented in figure 23.

Figure 23 .
Figure 23.Simultaneous or sequential resolution of CFX and ANSYS

Figure 27 .
Figure 27.Oscillatory response to sudden subject to a constant flow of 15m/s

Figure 29 .
Figure 29.Flutter response obtained from Matlab Aerospace blockset

Figure 30 .Figure 31 .
Figure 30.Damping coefficient against rotational speed and flutter frequency against rotor speed

Table 1 .
Description of the trials through homothetic transformation Aeroelasticity of Wind Turbines Blades Using Numerical Simulation http://dx.doi.org/10.5772/52281 Figure 15.Drag and lift coefficients vs. homothetic factor for different angles of attack

Table 3 .
[35]rpt from Table5in "Jennifer Heeg"[35]: Structural dynamic parameters associated with wind tunnel model configurations Table4lists the analytical calculations for divergence conditions for the considered model presented in[35].

Table 4 .
Analytical calculations for divergence conditions for the considered model presented