*α* =20<sup>0</sup> ± 5.5**sin**(**ωt**)<sup>0</sup>

Figure 21 illustrates the evolution of the aerodynamic coefficients with oscillation of the Ao Aaround 200 with amplitude of 5.50 and reduced frequency of 0.026. As the angle exceeds 140 , the purely turbulent k-ω SST model was used. For the drag coefficients, the results pro‐ vided by the k-ω SST model are quite close to the experimental results but predict prema‐ ture 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 turbu‐ lence. 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.

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

#### *4.1.3. 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 re‐ sults are obtained for low AoA but as turbulence intensity gets very large, the results di‐ verge 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 litera‐ ture, refinements should be brought to the models. Moreover, this study allows us to appre‐ ciate 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.

#### **4.2. Aerodynamic divergence**

In this section we will illustrate the different steps in modeling another aeroelastic phenom‐ enon, 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 phenom‐ ena 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 pa‐ rameters 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, *Ru*(*n*1) [34]:

$$\frac{\sigma\_{\text{x1}}}{\sigma\_{\text{x1}}} = \frac{\sigma\_{\text{u}}}{\sigma} \frac{\pi}{\sqrt{2\delta}} \sqrt{\mathbf{R}\_{\text{u}}(\mathbf{n1})} \sqrt{\mathbf{K}\_{\text{sc}}(\mathbf{n1})} \tag{32}$$

where:

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

the model prematurely predicts stall compared to the experiments.

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

*4.1.3. Conclusions on stallmodelling*

with amplitude of 5.50

Figure 21 illustrates the evolution of the aerodynamic coefficients with oscillation of the Ao

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 re‐ sults are obtained for low AoA but as turbulence intensity gets very large, the results di‐ verge 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 litera‐ ture, refinements should be brought to the models. Moreover, this study allows us to appre‐ ciate 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.

, the purely turbulent k-ω SST model was used. For the drag coefficients, the results pro‐ vided by the k-ω SST model are quite close to the experimental results but predict prema‐ ture 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 turbu‐ lence. 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,

and reduced frequency of 0.026. As the angle exceeds

*α* =20<sup>0</sup> ± 5.5**sin**(**ωt**)<sup>0</sup>

Aaround 200

106 Advances in Wind Power

140

$$R\_u(\text{n1}) = \frac{\, ^{n.S\_u} \text{(n1)}}{\sigma^2 \, ^{2} \text{(} \, ^{2}\text{)}} \tag{33}$$

**•** *x*1 is the first mode component of the steady tip displacement.


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 fluctuat‐ ing aerodynamic force acting on it per unit length is given by:

$$\frac{1}{2}\rho \begin{pmatrix} \mathbf{\dot{\cdot}} & \mathbf{\dot{\cdot}} \end{pmatrix}^{2} \mathbf{C}\_{d} \cdot \mathbf{C}(\mathbf{r}) \cdot \frac{1}{2}\rho \mathbf{\dot{\cdot}} \mathbf{\dot{I}}^{2} \mathbf{C}\_{d} \mathbf{C}(\mathbf{r}) \equiv \rho \mathbf{\dot{\cdot}} \mathbf{I} \, \dot{\mathbf{\dot{\cdot}}} \mathbf{C}\_{d} \mathbf{C} \Big|\_{\mathbf{r}} \tag{34}$$

where *x*˙ is the blade flatwise velocity, *Cd* is the drag coefficient and C(r) is the local blade chord. Hence the aerodynamic damping per unit length, *Ca* ^ (*r*)=*ρ<sup>U</sup>* - *CdC*(*r*) and the first aero‐ dynamic damping mode is:

$$\varepsilon\_{a1} = \frac{\mathcal{C}\_{a1}}{2m\_1\omega\_1} = \frac{\int\_0^R \hat{\mathcal{C}}\_a(r)\mu\_1^2(r)dr}{2m\_1\omega\_1} = \frac{\int\_0^r \mathcal{L}\mathcal{L}\_a\mathbb{I}\_0^R \mathcal{C}(r)\mu\_1^2(r)dr}{2m\_1\omega\_1} \tag{35}$$

*μ*1 (*r*) is the first mode shape and *m*1 is the generalized mass given by:

$$
\mu\_1 = \int\_0^R m(r) \mu\_1^2(r) dr \tag{36}
$$

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 oper‐ ating 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.

#### *4.2.1. 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 geome‐ tries 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 dif‐ ferent 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 dif‐ fusing 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.

**Figure 22.** Transfer Surfaces

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

where *x*˙ is the blade flatwise velocity, *Cd* is the drag coefficient and C(r) is the local blade


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 oper‐ ating 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

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 geome‐ tries 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 dif‐ ferent 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‐

2*m*1*ω*<sup>1</sup>

*Cd* .*C*(*r*)≅*ρU*


^ (*r*)=*ρ<sup>U</sup>* -

2(*r*)*dr* (36)

*<sup>x</sup>*˙*CdC*(*r*) (34)

*CdC*(*r*) and the first aero‐

(35)

<sup>2</sup> *ρU* - 2

ing aerodynamic force acting on it per unit length is given by:

chord. Hence the aerodynamic damping per unit length, *Ca*

*<sup>ε</sup>a*<sup>1</sup> <sup>=</sup> *Ca*<sup>1</sup> <sup>2</sup>*m*1*ω*<sup>1</sup> <sup>=</sup> *<sup>∫</sup>* 0 *RCa* ^ (*r*)*μ*<sup>1</sup> 2(*r*)*dr* <sup>2</sup>*m*1*ω*<sup>1</sup> <sup>=</sup> *<sup>ρ</sup><sup>U</sup>*

aeroelastic phenomena are also approached.

*4.2.1. ANSYS-CFX coupling*

*Cd* - *<sup>C</sup>*(*r*) - <sup>1</sup>

(*r*) is the first mode shape and *m*1 is the generalized mass given by:

*m*<sup>1</sup> = *∫* 0 *Rm*(*r*)*μ*<sup>1</sup>

1 <sup>2</sup> *<sup>ρ</sup>*(*<sup>U</sup>* - - *x*˙ ) 2

dynamic damping mode is:

108 Advances in Wind Power

*μ*1

The respective MFX simultaneous and sequential resolution schemes are presented in fig‐ ure 23.

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

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 so‐ licitation 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 il‐ lustrates such interaction is illustrated in Figure 24 below:

**Figure 24.** ANSYS workbench divergence flow-chart

#### *4.2.2. Comparison with experimental results*

#### *4.2.2.1. 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.

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


**Table 3.** Excerpt from Table 5 in "Jennifer Heeg" [35]: Structural dynamic parameters associated with wind tunnel model configurations

Table 4 lists 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

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 succes‐ sively modified until a model as close as possible to the model in the experiment was ob‐ tained.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 sever‐ al 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.

**Figure 25.** Divergence of wind tunnel model configuration #2

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 pres‐ sure was being slowly increased until the angle of attack increased dramatically and sud‐ denly. 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].

#### *4.2.2.3. The ANSYS-CFX model*

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

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.

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

> **α (Hz)**

**Table 3.** Excerpt from Table 5 in "Jennifer Heeg" [35]: Structural dynamic parameters associated with wind tunnel

Table 4 lists the analytical calculations for divergence conditions for the considered model

**ζ**

ceiling. The structural dynamic parameters for this model are illustrated in table 3:

5.8262 49.5 7.88 0.053

**ωα (rads/sec)**

lustrates such interaction is illustrated in Figure 24 below:

**Figure 24.** ANSYS workbench divergence flow-chart

110 Advances in Wind Power

*4.2.2. Comparison with experimental results*

*4.2.2.1. Overview of the experimental results*

*4.2.2.2. Configuration description*

**Kα (N∙m/rad)**

model configurations

presented in [35].

The model used in the experiment was simulated using a reduced span-wise numerical do‐ main (quasi 2D). The span of the airfoil was reduced 262.5 times, from 21 inches to 0.08 in‐ ches 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.

**Figure 26.** ANSYS built geometry with meshing

#### *4.2.2.4. Results*

In [23], the authors have derived the analytical mathematical equation to calculate the diver‐ gence velocity, *UD*. The expression was:

$$
\Delta U\_D = \sqrt{\frac{1}{C^{\alpha 0} \frac{\stackrel{\text{J} \cdot C}{\delta \cdot a}}{\delta \cdot a} e^{\frac{1}{\Delta} \rho S}}} \tag{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 ∂*Cl* <sup>∂</sup> *<sup>α</sup>* which depends both on the shape of the pro‐ file 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]. How‐ ever, numerical values are finite and cannot model infinite parameters. We will, therefore, for‐ mulate 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:

$$\Theta\_r = \mathbb{C}^{\theta \theta} T = \mathbb{C}^{\theta \theta} (\frac{\partial \, \mathbb{C}\_l}{\partial \alpha} \text{ e } \alpha\_r + \mathbb{C}\_m \text{c}) qS \tag{38}$$

Algebraic manipulations of the expressions lead us to the following formulation:

$$\Theta = \frac{\theta\_r}{1 - \frac{\partial \mathcal{C}\_l}{\partial a} \mathcal{C}^{\theta \theta} \mathcal{e} \eta \mathcal{S}} \tag{39}$$

This leads to:

$$\Theta = \frac{\theta\_r}{1 \cdot \frac{q}{q\_D}} = \frac{\theta\_r}{1 \cdot \left(\frac{U}{U\_D}\right)^2} \tag{40}$$

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 supple‐ mentary 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 in‐ termittency 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<sup>α</sup> = 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:

$$\mathbf{K}\_{ar} = \frac{5.8262}{262.5} \mathbf{N} \bullet \frac{m}{rad} = 0.022195 \ \mathbf{N} \bullet \ \mathbf{m} \Big/ \ \text{rad} \tag{41}$$

and the flexibility coefficient is:

**Figure 26.** ANSYS built geometry with meshing

gence velocity, *UD*. The expression was:

which is dependent upon the used fluid and

In [23], the authors have derived the analytical mathematical equation to calculate the diver‐

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, *ρ*,

∂*Cl*

file 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]. How‐ ever, numerical values are finite and cannot model infinite parameters. We will, therefore, for‐ mulate 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:

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

> <sup>=</sup> *<sup>θ</sup><sup>r</sup>* <sup>1</sup> - ( *<sup>U</sup>*

*<sup>θ</sup>* <sup>=</sup> *<sup>θ</sup><sup>r</sup>* <sup>1</sup> - *<sup>q</sup> qD*

*<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>* which depends both on the shape of the pro‐

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

*<sup>U</sup> <sup>D</sup>* )<sup>2</sup> (40)

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

*4.2.2.4. Results*

112 Advances in Wind Power

This leads to:

$$\text{C}^{\prime \Theta \theta} = \frac{1}{K\_{ar}} = \text{45.0552 } \, rad \left/ \text{N} \bullet m \right. \tag{42}$$

The slope of the lift ∂*Cl* <sup>∂</sup> *<sup>α</sup>* , can be calculated for an angle *<sup>α</sup>* =5<sup>0</sup> in the following way:

$$\frac{\partial \, \mathcal{C}\_l}{\partial \alpha} = \frac{\mathcal{C}\_{l,\alpha \succ \mathbb{S}^0} \cdot \mathcal{C}\_{l,\alpha \prec \mathbb{S}^0}}{\alpha \rhd \mathfrak{S}^0 \cdot \alpha \prec \mathfrak{S}^0} \tag{43}$$

We have calculated the lift coefficient at 4.0<sup>0</sup> and 6.0<sup>0</sup> such that:

$$\begin{array}{l} \text{C}\_{l,\alpha=4.0^{\circ}} = 0.475\\ \text{and } \text{C}\_{l,\alpha=6.0^{\circ}} = 0.703 \end{array}$$

(37)

(39)

Hence the gradient can be expressed and calculated as follows:

$$\frac{\partial C\_l}{\partial \alpha} = \frac{0.703 - 0.475}{6.0 - 4.0} = 0.114 \text{ deg}^{-1} = 6.532 \text{ rad}^{-1}$$

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 calcu‐ lated as follows:

*S* =0.2032•0.5334=0.0004129024*m*<sup>2</sup>

Hence the divergence velocity is calculated as:

$$\mathcal{U}\_D = \sqrt{\frac{1}{c^{\frac{\alpha \mathcal{S}\_l}{\alpha \mathcal{S}\_a} \frac{\rho}{2} \epsilon \mathcal{S}}}} = 18.78 \text{ m/s} \tag{44}$$

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 ∂*Cl* <sup>∂</sup> *<sup>α</sup>* taken into considera‐ tion, 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 in‐ teraction. Furthermore, using the model, domain and mesh parameters detailed in the previ‐ ous 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 com‐ pared 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].

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

#### **4.3. 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 mod‐ eling. In the divergence section, the example was used to reinforce the notion of multiphy‐ sics 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.

#### *4.3.1. 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 enor‐ mous 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.

#### *4.3.2. Matlab-Simulink and Ansys-CFX tools*

The theoretical divergence speed given in Table 4 of the NASA experiment [35] is 19.15 m/s.

tion, which in the NASA work was 2π, or 6.283 rad-1, whereas we used a value of 6.532

pared 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

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 mod‐ eling. In the divergence section, the example was used to reinforce the notion of multiphy‐ sics 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

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

 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 in‐ teraction. Furthermore, using the model, domain and mesh parameters detailed in the previ‐ ous 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

. Suddenly, the fixing is removed and the constant flow can be then com‐

∂*Cl*

<sup>∂</sup> *<sup>α</sup>* taken into considera‐

This slight difference is due to the value of slope of the lift profile

Hz which is close to the 7.9 Hz frequency presented in [35].

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

*4.3.1. Computational requirement and Lumped model*

rad-1 .

114 Advances in Wind Power

velocity 15 m *s*-1

**4.3. Aerodynamic flutter**

lumped method.

For flutter modelling, again, ANSYS-CFX model was used to simulate the complex fluidstructure interaction. However, due to excessively important computational time that ren‐ dered 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 em‐ bedded 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 mod‐ els, and aerospace analysis pre-programmed tools as well as aerodynamic coefficient im‐ porting 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:

$$\psi\_{\rm lo} = \frac{\sigma^{\frac{\sigma}{\omega}\frac{2}{\omega}}}{V^{\frac{\sigma}{\omega}}\omega} \cdot \frac{0.8 \left(\frac{\pi L}{4b}\right)^{0.3}}{1 + \left(\frac{4b\omega}{\pi V}\right)^{2}}\tag{45}$$

Here,Lω represents the turbulence scale length and σis the turbulence intensity. The corre‐ sponding transfer function used in Simulink is:

$$\psi\_{\rm lo} = \frac{\sigma\_u \sqrt{\frac{2}{\pi} \frac{L\_v}{V}} \Big|\_{1 \to 0.25 \frac{L\_v}{V} s}}{1 + 1.357 \frac{L\_v}{V} s + 0.1987 \left(\frac{L\_v}{V} s\right)^2 s^2} \tag{46}$$

For the lateral speed, the turbulence spectrum is defined as:

$$\psi\_{\rm la} = \frac{\mp \left(\frac{\omega}{\nabla}\right)^2}{1 + \left(\frac{3\kappa\omega}{\pi V}\right)^2} \cdot \varphi\_{\upsilon}(\omega) \tag{47}$$

and the corresponding transfer function can be expressed as :

$$\psi\_{\rm la} = \frac{\frac{\pi}{4} \left(\frac{s}{V}\right)^{\rm l}}{1 + \left(\frac{3b}{\pi V}\right)^{\rm l}} H\_v \text{(s)}\tag{48}$$

Finally, the vertical turbulence spectrum is expressed as follows:

<sup>ψ</sup><sup>v</sup> <sup>=</sup> <sup>∓</sup>( *<sup>ω</sup> V* )2 <sup>1</sup> <sup>+</sup> ( <sup>4</sup>*b<sup>ω</sup> <sup>π</sup><sup>V</sup>* )<sup>2</sup> .*φω*(*ω*) (49)

and the corresponding transfer function is expressed as follows:

$$\psi\_{\mathbf{v}} = \frac{\overline{\frac{\mathbf{r}}{\mathbf{r}} \begin{pmatrix} \frac{\mathbf{s}}{\mathbf{r} \cdot \mathbf{r}} \end{pmatrix}^{\mathbf{l}}}{\mathbf{1} + \begin{pmatrix} \frac{4\mathbf{b}}{\pi \mathbf{v} \cdot \mathbf{s}} \end{pmatrix}^{\mathbf{l}}}.\tag{50}$$

The Aerodynamic Forces and Moments block computes the aerodynamic forces and mo‐ ments around the center of gravity. The net rotation from body to wind axes is expressed as:

$$\mathbf{C}\_{\omega \gets b} = \begin{bmatrix} \cos \left( \alpha \right) \cos \left( \beta \right) & \sin \left( \beta \right) & \sin \left( \alpha \right) \cos \left( \beta \right) \\ -\cos \left( \alpha \right) \sin \left( \beta \right) & \cos \left( \beta \right) & -\sin \left( \alpha \right) \sin \left( \beta \right) \\ -\sin \left( \alpha \right) & 0 & \cos \left( \alpha \right) \end{bmatrix} \tag{51}$$

On the other hand, the fluid structure interaction to model aerodynamic flutter was made using ANSYS multi domain (MFX). As previously mentioned, the drawback of the AN‐ SYS model is that it is very time and memory consuming. However, it provides a very good option to compare and validate simplified model results and understand the intrin‐ sic theories of flutter modelling. On one hand, the aerodynamics of the application is modelled using the fluid module CFX and on the other side, the dynamic structural part is modelled using ANSYS structural module. An iterative exchange of data between the two modules to simulate the flutter phenomenon is done using the Workbench interface.

#### *4.3.3. Lumped model results*

We will first present the results obtained by modeling AoA for configuration # 2 of reference [35] (also, discussed in the divergence section 4.2) for an initial AoA of 0°. As soon as diver‐ gence is triggered, within 1 second the blade oscillates in a very spectacular and dangerous manner. This happens at a dynamic pressure of 5,59lb/pi2 (268 N/m2 ). Configuration #2 uses, on the airfoil, 20 elements, unity as the normalized element size and unity as the normalized airfoil length. Similarly, the number of elements in the wake is 360 and the corresponding normalized element size is unity and the normalized wake length is equal to 2. The results obtained in [35] are illustrated in Figure 28:

**Figure 28.** Flutter response- an excerpt from [23]

<sup>ψ</sup>la <sup>=</sup> <sup>∓</sup>( *<sup>s</sup>*

<sup>ψ</sup><sup>v</sup> <sup>=</sup> <sup>∓</sup>( *<sup>ω</sup>*

<sup>ψ</sup><sup>v</sup> <sup>=</sup> <sup>∓</sup>( *<sup>s</sup>*

Finally, the vertical turbulence spectrum is expressed as follows:

and the corresponding transfer function is expressed as follows:

Cω←<sup>b</sup> =

116 Advances in Wind Power

*4.3.3. Lumped model results*

manner. This happens at a dynamic pressure of 5,59lb/pi2

obtained in [35] are illustrated in Figure 28:

*V* )1 <sup>1</sup> <sup>+</sup> ( <sup>3</sup>*<sup>b</sup>*

*V* )2 <sup>1</sup> <sup>+</sup> ( <sup>4</sup>*b<sup>ω</sup>*

*V* )1 <sup>1</sup> <sup>+</sup> ( <sup>4</sup>*<sup>b</sup>*

The Aerodynamic Forces and Moments block computes the aerodynamic forces and mo‐ ments around the center of gravity. The net rotation from body to wind axes is expressed as:

On the other hand, the fluid structure interaction to model aerodynamic flutter was made using ANSYS multi domain (MFX). As previously mentioned, the drawback of the AN‐ SYS model is that it is very time and memory consuming. However, it provides a very good option to compare and validate simplified model results and understand the intrin‐ sic theories of flutter modelling. On one hand, the aerodynamics of the application is modelled using the fluid module CFX and on the other side, the dynamic structural part is modelled using ANSYS structural module. An iterative exchange of data between the two modules to simulate the flutter phenomenon is done using the Workbench interface.

We will first present the results obtained by modeling AoA for configuration # 2 of reference [35] (also, discussed in the divergence section 4.2) for an initial AoA of 0°. As soon as diver‐ gence is triggered, within 1 second the blade oscillates in a very spectacular and dangerous

on the airfoil, 20 elements, unity as the normalized element size and unity as the normalized airfoil length. Similarly, the number of elements in the wake is 360 and the corresponding normalized element size is unity and the normalized wake length is equal to 2. The results

(268 N/m2

). Configuration #2 uses,

cos (α)cos (β) sin (β) sin (α)cos (β) -cos (α)sin (β) cos (β) -sin (α)sin (β) -sin (α) 0 cos (α)

*<sup>π</sup><sup>V</sup> <sup>s</sup>*)<sup>1</sup> .*Hv*(*s*) (48)

*<sup>π</sup><sup>V</sup>* )<sup>2</sup> .*φω*(*ω*) (49)

*<sup>π</sup><sup>V</sup> <sup>s</sup>*)<sup>1</sup> .*Hω*(*s*) (50)

(51)

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<sup>0</sup> ± 17<sup>0</sup> . . The same modelling was performed using the Simulink model and the result for the AoA varia‐ tion and the plunge displacement is shown below:

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

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<sup>0</sup> ± 17<sup>0</sup> 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: ζ<sup>=</sup> <sup>c</sup> 2m<sup>ω</sup> , ω is measured as the Laplace integral in Simulink, c is the viscous damping and ω= k <sup>m</sup> . 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 rota‐ tion 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 rea‐ son for which flutter is usually very similar to resonance as it occurs due to a coalescing of dy‐ namic modes close to the natural vibrating mode of the system.

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

**Figure 31.** Flutter simulation with ANSYS-CFX at 1) 1.8449 s, 2) 1.88822 s and 1.93154s

We present here the results obtained for the same case study using ANSYS-CFX. The fre‐ quency 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 maxi‐ mum 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 an‐ gle 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 mo‐ ment 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.
