**3. 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 ana‐ lytical models. Few references related to analytic developments of aeroelastic phenomena are [10], [11] and [12].

#### **3.1. 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 shed‐ ding to travel from the leading edge backward along the airfoil [14]. The analytical develop‐ ment 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:

$$L = c\_L \, \frac{1}{2} \rho V^2 c \,\tag{7}$$

Where


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

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 linear‐ ly with α. Stall happens at a critical angle of attack whereby the lift reaches a maximum val‐ ue and flow separation on the suction side occurs.

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

For unsteady conditions, a delay exists prior to reaching stability and is an essential condi‐ tion 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 sig‐ nificant 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 var‐ iations of the AoA between 0o and 15o,, the onset of stall is delayed and the lift is considera‐ bly 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 forma‐ tion 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.

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

**Figure 7.** Aerodynamic stall mechanism- Vortex creation at the leading edge

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 linear‐ ly with α. Stall happens at a critical angle of attack whereby the lift reaches a maximum val‐

**Figure 5.** Lift coefficient under static and dynamic stall conditions (dashed line for steady conditions, plain line for un‐

For unsteady conditions, a delay exists prior to reaching stability and is an essential condi‐ tion 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 sig‐ nificant 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 var‐

bly 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 forma‐ tion 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.

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

and 15o,, the onset of stall is delayed and the lift is considera‐

ue and flow separation on the suction side occurs.

steady conditions)

92 Advances in Wind Power

iations of the AoA between 0o

**Figure 8.** Aerodynamic stall mechanism- Vortex separation at the leading edge and creation of vortices at the trailing edge

**Figure 9.** Aerodynamic stall mechanism – Vortex shedding at the trailing edge

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 pa‐ rameters, study of flow transport, boundary layer analysis (shape factor and thickness), vor‐ tex 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 aer‐ oelastic consideration. The proper modelling of transition from laminar to turbulent flow is also essential for accurate prediction of stall parameters.

#### **3.2. 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 air‐ speed 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.

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

The airspeed at which the elastic twist increases rapidly to the point of failure is called the divergence airspeed, *UD*.. This phenomenon, being highly dangerous and prejudicial for wind blades, makes the accurate calculation of *UD* 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.

$$
\alpha = \alpha\_r + \Theta \tag{8}
$$

The elastic twist angle is proportional to the moment at the elastic axis:

$$
\Theta = \mathbb{C}^{\prime \theta \theta} T \tag{9}
$$

where *C θθ* is the flexibility coefficient of the spring. The total aerodynamic moment with re‐ spect to the elastic axis is given by:

$$T = \{\mathbf{C}\_l e + \mathbf{C}\_m c\} qS \tag{10}$$

where

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 pa‐ rameters, study of flow transport, boundary layer analysis (shape factor and thickness), vor‐ tex 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 aer‐ oelastic consideration. The proper modelling of transition from laminar to turbulent flow is

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 air‐ speed 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

The airspeed at which the elastic twist increases rapidly to the point of failure is called the divergence airspeed, *UD*.. This phenomenon, being highly dangerous and prejudicial for wind blades, makes the accurate calculation of *UD* 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

*α* = *α<sup>r</sup>* + θ (8)

also essential for accurate prediction of stall parameters.

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

The elastic twist angle is proportional to the moment at the elastic axis:

**3.2. Divergence**

94 Advances in Wind Power

to structural failure.

as the elastic twist angle.


The lift coefficient is related to the angle of attack measured with respect to a zero lift condi‐ tion as follows:

$$C\_l = \frac{\partial C\_l}{\partial \alpha} \left(\alpha\_r + \Theta\right) \tag{11}$$

Here ∂*Cl* <sup>∂</sup> *<sup>α</sup>* represents the slope of the lift curve. The elastic twist angle *θ*, can be obtained by simple mathematical manipulations of the three previous equations:

$$T = \begin{bmatrix} \frac{\partial \mathcal{L}\_l}{\partial \alpha} (\alpha\_r + \Theta) . \mathcal{e} + \mathcal{C}\_m \mathcal{c} \end{bmatrix} q \mathbf{S} \tag{12}$$

Hence,

$$\mathcal{Q} = \mathbb{C}^{\partial \mathcal{O}} \Big[ \frac{\partial \mathcal{C}\_l}{\partial a} (\alpha\_r + \mathcal{O}) . e + \mathcal{C}\_m c \Big] \mathcal{q} S \tag{13}$$

$$\Theta = \mathbb{C}^{\theta \theta} \Big[ \underbrace{^{\partial \mathcal{L}\_l}}\_{\partial \alpha} \alpha\_r eqS + \frac{\partial \mathcal{L}\_l}{\partial \alpha} \Theta eqS + \mathcal{C}\_m cqS \Big] \tag{14}$$

Regrouping *θ* :

$$\left[\Theta\right]\mathbf{1} - \frac{\partial \mathbf{C}\_l}{\partial a}\mathbf{C}^{\partial \theta} \epsilon q \mathbf{S} \mathbf{J} = \mathbf{C}^{\partial \theta} \left[\frac{\partial \mathbf{C}\_l}{\partial a} \alpha\_r \epsilon q \mathbf{S} + \mathbf{C}\_m \epsilon q \mathbf{S}\right] \tag{15}$$

This leads to:

$$\Theta = \mathbb{C}^{\Theta \Theta} \frac{\stackrel{\circ \, \mathcal{C}\_l}{\stackrel{\circ \, \mathcal{C}\_l}{\circ \,}} \stackrel{\circ}{\alpha\_r} \epsilon q \mathbb{S} + \mathcal{C}\_m \epsilon q \mathbb{S}}{\stackrel{\circ \, \mathcal{C}\_l}{1 - \stackrel{\circ \, \mathcal{C}\_l}{\circ \,}} \mathcal{C}^{\{\theta\}} \epsilon q \mathbb{S}} \tag{16}$$

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 aer‐ odynamic divergence. When the denominator tends to zero:

$$\mathbf{1} - \frac{\partial \mathcal{C}\_{\parallel}}{\partial a} \mathbf{C}^{\prime \theta \theta} \mathbf{e} \mathbf{g} \mathbf{S} = \mathbf{0} \tag{17}$$

The dynamic pressure is given by:

$$q = \frac{1}{2}\rho v^2\tag{18}$$

Thus, we come up with:

$$\mathbf{1} \cdot \frac{\partial \mathcal{C}\_l}{\partial a} \mathbf{C}^{\theta \theta} e \frac{1}{2} \rho v^2 \mathbf{S} = \mathbf{0} \tag{19}$$

Hence, the divergence velocity can be expressed as:

$$\mathcal{U}\mathcal{U}\_D = \sqrt{\frac{1}{C^{\alpha\theta} \overbrace{\stackrel{\partial C\_\eta}{\rightleftharpoons}} \epsilon \frac{1}{\varDelta t} \rho S}}\tag{20}$$

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 ∂*Cl* <sup>∂</sup> *<sup>α</sup>* depends both on the shape of the airfoil and flow conditions [23]. We note that as divergence velocity is approached, the elastic twist an‐ gle 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:

$$
\partial\_r \Theta\_r = \mathbb{C}^{\partial \theta} T = \mathbb{C}^{\partial \theta} (\frac{\partial \mathbb{C}\_l}{\partial \alpha} e \; \alpha\_r + \mathbb{C}\_m \mathbb{C}) q \mathbb{S} \tag{21}
$$

Hence:

$$
\Theta = \frac{\theta\_r}{1 - \overbrace{^{\text{d'C}\_l}\text{C}^{\text{e@\theta}}\text{eqS}}^{\theta\_r}} \tag{22}
$$

which leads to:

Aeroelasticity of Wind Turbines Blades Using Numerical Simulation http://dx.doi.org/10.5772/52281 97

$$\Theta = \frac{\theta\_r}{1 \cdot \frac{q}{q\nu}} = \frac{\theta\_r}{1 \cdot \left(\frac{l\ell}{l\ell\_D}\right)^2} \tag{23}$$

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 aero‐ dynamic solicitation. The latter is calculated by solving for the moment applied on the pro‐ file 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 di‐ vergence speed are detailed in [28]:

$$\mathcal{U}\mathcal{U}\_D = \sqrt{\frac{1}{C^{\circ 0} \frac{\text{\reflectbox{ $\circ$ }}{\text{\reflectbox{ $\circ$ }}}{\text{\reflectbox{ $\circ$ }}} \frac{\text{\reflectbox{ $\circ$ }}}{\text{\reflectbox{ $\circ$ }}} eS}}\tag{24}$$

Detailed eigen values and eigenvectors analysis related to divergence phenomenon is pre‐ sented in [27].

#### **3.3. Aerodynamic flutter**

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 aer‐

<sup>∂</sup> *<sup>α</sup> <sup>C</sup> θθeqS* =0 (17)

<sup>2</sup> *<sup>ρ</sup><sup>v</sup>* <sup>2</sup> (18)

*S* =0 (19)

<sup>∂</sup> *<sup>α</sup>* depends both on the shape of the airfoil and

<sup>∂</sup> *<sup>α</sup> e α<sup>r</sup>* + *Cmc*)*qS* (21)

(20)

(22)

odynamic divergence. When the denominator tends to zero:

The dynamic pressure is given by:

Thus, we come up with:

96 Advances in Wind Power

1 - ∂*Cl*

1 - ∂*Cl* <sup>∂</sup> *<sup>α</sup> <sup>C</sup> θθ<sup>e</sup>*

Hence, the divergence velocity can be expressed as:

airfoil, *ρ* depends upon the used fluid and

solving for the elastic twist angle, we have:

Hence:

which leads to:

*<sup>q</sup>* <sup>=</sup> <sup>1</sup>

*UD* <sup>=</sup> <sup>1</sup> *<sup>C</sup> θθ* <sup>∂</sup> *Cl* <sup>∂</sup> *<sup>α</sup> e* 1 <sup>2</sup> *ρS*

*<sup>θ</sup><sup>r</sup>* <sup>=</sup>*<sup>C</sup> θθ<sup>T</sup>* <sup>=</sup>*<sup>C</sup> θθ*( <sup>∂</sup>*Cl*

*<sup>θ</sup>* <sup>=</sup> *<sup>θ</sup><sup>r</sup>* 1 - ∂ *Cl* <sup>∂</sup> *<sup>α</sup> C θθeqS*

1 <sup>2</sup> *<sup>ρ</sup><sup>v</sup>* <sup>2</sup>

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

∂*Cl*

flow conditions [23]. We note that as divergence velocity is approached, the elastic twist an‐ gle 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‐

As previously mentioned, flutter is caused by the superposition of two structural modes – pitch 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-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 posi‐ tioned 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, d1 and d2 are geometrical quantities as shown in Figure 11.

**Figure 11.** Model defining parameters

$$L = 2\pi\rho L\,^2b\Big|\frac{i\omega C(k)h\_0}{L} + C\{k\}a\_0 + \left[1 + C\{k\}(1 - 2a)\right]\frac{i\omega ba\_0}{2L} - \frac{\omega^2 b h\_0}{2L^2} + \frac{\omega^2 b^2 a a\_0}{2L^2}\Big|\tag{25}$$

$$\Delta M = 2\pi \rho \lambda I^2 b \left| d\_1 \right| \xrightarrow{i\omega C(k)h\_0} + \mathcal{C}(k)a\_0 + \left[1 + \mathcal{C}(k)(1 - 2a)\right] \xrightarrow{i\omega h\_0} d\_2 \xrightarrow{i\omega h\_0} -\frac{\omega^2 b^2 a}{2L^2} h\_0 + \left(\frac{1}{\mathcal{S}} + a^2\right) \frac{\omega^2 b^2 \kappa\_0}{2L^2} \tag{26}$$

Theodorsen's equation can be rewritten in a form that can be used and analyzed in Matlab Simulink as follows:

$$L = 2\pi \rho L^2 b \left[ \frac{\mathbb{C}(k)}{\mathcal{U}} \dot{h} + \mathbb{C}(k) \infty + \left[ 1 + \mathbb{C}(k)(1 - 2a) \right] \right] \frac{b}{2\mathcal{U}} \dot{\infty} \\ \stackrel{b}{\to} + \frac{b}{2\mathcal{U}^2} \ddot{h} + \frac{b^2 a}{2\mathcal{U}^2} \overset{\cdots}{\infty} \tag{27}$$

$$M = 2\pi \rho L^2 b \left| d\_i \right| \frac{C(k)h^\cdot}{L^2} + C(k) \approx \left[ 1 + C(k) \left( 1 - 2a \right) \frac{b}{2L} \dot{\mathbf{c}} \times \dot{\mathbf{k}} \right] + d\_2 \frac{b}{2L} \dot{\mathbf{c}} \dot{\mathbf{c}} \times \frac{ab^2}{2L^2} \ddot{\mathbf{h}} \left| - \left( \frac{1}{8} + a^2 \right) \frac{b^2 \dot{\mathbf{c}}}{2L^2} \right| \tag{28}$$

#### *3.3.1. Flutter movement*

The occurrence of the flutter has been illustrated in Section 2 (Figure 3). To better under‐ stand 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.

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

#### *3.3.2. 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:

$$\alpha(\mathbf{x}, \ y, t) = \Theta\_T + \Theta(t) + \frac{\mathbb{H}(t)}{\mathcal{U}\_0} + \frac{l(\mathbf{x})\theta(t)}{\mathcal{U}\_0} \text{ - } \frac{w\_\mathbf{g}(\mathbf{x}, \ y, t)}{\mathcal{U}\_0} \tag{29}$$

The Lagrangian form equations are constructed for the mechanical system. The first one cor‐ responds to the vertical displacement z and the other is for the angle of attack α:

$$J\_0 \ddot{\alpha} + m d \cos(\alpha) \ddot{\hat{z}} + c(\alpha \cdot \alpha\_0) = M\_0 \tag{30}$$

$$m\ddot{z} + md\cos(\alpha)\ddot{\alpha} \text{ - } m\sin(\alpha)\dot{\alpha}^2 + kz = F\_Z \tag{31}$$

Numerical solution of these equations requires expressing Fz and Mo as polynomials of α. Moreover, *Fz* (*α*)= <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>SV* <sup>2</sup> *Cz* (*α*) and *Mo* (*α*)= <sup>1</sup> <sup>2</sup> *<sup>ρ</sup>LSV* <sup>2</sup> *Cm*<sup>0</sup> (*α*) for S being the surface of the blade, Cz, the lift coefficient, Cm0 being the pitch coefficient, Fz being the lift, Mo, the pitch moment. Cz and Cm values are extracted from NACA 4412 data. Third degree interpolations for Cz and Cm with respect to the AoA are given below:

$$C\_z = -0.0000983 \text{ } \alpha^3 \text{ - } 0.0003562 \alpha^2 + 0.1312 \alpha + 0.4162$$

$$C\_{m0} = -0.00006375 \alpha^3 + 0.00149 \alpha^2 - 0.001185 \text{ } \alpha \text{ - } 0.9312$$

These equations will be used in the modeling of a lumped representation of flutter present‐ ed in the last section of this chapter.
