**4. Computational fluid dynamics (CFD) methods in aeroelastic modeling**

Aeroelastic modeling of wind blades require complex representation of both fluid flows, in‐ cluding 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 ad‐ vent of high efficiency computers, and improvement in numerical techniques, computation‐ al 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 pop‐ ular 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 anal‐ ysis is based on modeling using ANSYS and CFX software. CFX uses a finite volume meth‐ od 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 aer‐ oelastic model. These calibration considerations will be illustrated in the stall modeling sec‐ tion as an example.

#### **4.1. Dynamic stall**

*L* =2*πρU* <sup>2</sup>

*b*{*d*<sup>1</sup>

*M* =2*πρU* <sup>2</sup>

98 Advances in Wind Power

Simulink as follows:

*L* =2*πρU* <sup>2</sup>

*M* =2*πρU* <sup>2</sup>

*3.3.1. Flutter movement*

*<sup>b</sup>*{ *<sup>i</sup>ωC*(*<sup>k</sup>* )*<sup>h</sup>* <sup>0</sup>

*<sup>i</sup>ωC*(*<sup>k</sup>* )*<sup>h</sup>* <sup>0</sup>

*<sup>b</sup>*{ *<sup>C</sup>*(*<sup>k</sup>* )

*b*{*d*<sup>1</sup> *C*(*k* )*h*˙ .

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

angle of attack of the model that can be written as:

*α*(*x*, *y*, *t*)=*θ<sup>T</sup>* + *θ*(*t*) +

*3.3.2. Flutter equations*

*<sup>U</sup>* <sup>+</sup> *<sup>C</sup>*(*k*)*α*<sup>0</sup> <sup>+</sup> 1 + *<sup>C</sup>*(*k*)(1-2*a*) *<sup>i</sup>ωbα*<sup>0</sup>

<sup>2</sup>*<sup>U</sup>* + *d*<sup>2</sup>

<sup>2</sup>*<sup>U</sup>* ∝˙ *<sup>k</sup>* <sup>+</sup> *<sup>d</sup>*<sup>2</sup>

Theodorsen's equation can be rewritten in a form that can be used and analyzed in Matlab

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

The flutter equations originate in the relation between the generalized coordinates and the

The Lagrangian form equations are constructed for the mechanical system. The first one cor‐

*l*(*x*)*θ* ˙(*t*) *<sup>U</sup>*<sup>0</sup> - *wg* (*x*, *y*, *t*) *U*<sup>0</sup>

*h*˙ (*t*) *<sup>U</sup>*<sup>0</sup> +

responds to the vertical displacement z and the other is for the angle of attack α:

rigidity of the structure whereas the rotation spring represents the rotation rigidity.

*<sup>U</sup>* <sup>+</sup> *<sup>C</sup>*(*k*)*α*<sup>0</sup> <sup>+</sup> 1 + *<sup>C</sup>*(*k*)(1-2*a*) *<sup>i</sup>ωbα*<sup>0</sup>

*<sup>U</sup>* <sup>+</sup> *<sup>C</sup>*(*k*)<sup>∝</sup> + 1 + *<sup>C</sup>*(*k*)(1-2*a*) *<sup>b</sup>*

*<sup>U</sup> <sup>h</sup>*˙ <sup>+</sup> *<sup>C</sup>*(*k*)<sup>∝</sup> <sup>+</sup> 1 + *<sup>C</sup>*(*k*)(1-2*a*) *<sup>b</sup>*

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

<sup>2</sup>*<sup>U</sup>* ∝˙ <sup>+</sup>

*b* <sup>2</sup>*<sup>U</sup>* ∝˙ <sup>+</sup> *ab*<sup>2</sup> <sup>2</sup>*<sup>U</sup>* <sup>2</sup> *<sup>h</sup>* ¨ - ( <sup>1</sup> <sup>8</sup> <sup>+</sup> *<sup>a</sup>* 2) *<sup>b</sup>* <sup>3</sup> ∝¨

*b* <sup>2</sup>*<sup>U</sup>* <sup>2</sup> *h* ¨ - *<sup>b</sup>* <sup>2</sup> *a* <sup>2</sup>*<sup>U</sup>* <sup>2</sup> ∝

*iωbα*<sup>0</sup> <sup>2</sup>*<sup>U</sup>* - *<sup>ω</sup>* <sup>2</sup> *b* 2 *a* <sup>2</sup>*<sup>U</sup>* <sup>2</sup> *<sup>h</sup>*<sup>0</sup> <sup>+</sup> ( <sup>1</sup>

*ω* <sup>2</sup> *bh* <sup>0</sup> <sup>2</sup>*<sup>U</sup>* <sup>2</sup> <sup>+</sup> *<sup>ω</sup>* <sup>2</sup>

*b* 2 *aα*<sup>0</sup>

<sup>8</sup> <sup>+</sup> *<sup>a</sup>* 2) *<sup>ω</sup>* <sup>2</sup> *b* 3 ∝<sup>0</sup>

<sup>2</sup>*<sup>U</sup>* <sup>2</sup> } (25)

<sup>2</sup>*<sup>U</sup>* <sup>2</sup> } (26)

¨ } (27)

<sup>2</sup>*<sup>U</sup>* <sup>2</sup> } (28)

(29)

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 aerody‐ namic 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 fol‐ lowing 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 solici‐ tations. The MFX module allows fluid-structure modelling, i.e., the results of the aerody‐ namic 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.

#### *4.1.1. 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 mod‐ el, an S809 profile, designed by NREL, was used.

**Figure 13.** S809 airfoil

This airfoil has been chosen as experimental results and results from other sources are avail‐ able for comparison. The experimental results have been obtained at the Low Speed Labora‐ tory 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 character‐ istics of the S809 profile for angles of incidence from -200 to 200 . The second study [32], real‐ ized by Ramsey, gives the characteristics of the airfoil for angles of incidence ranging from -200 to 400 . 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 (80 , 140 and 200 ) oscillating at (±5.50 and ± 100 ) at different frequencies for Reynolds numbers between 0.75 and 1.4 million.

### *4.1.1.2. Convergence studies*

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 aerody‐ namic 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 fol‐ lowing 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 solici‐ tations. The MFX module allows fluid-structure modelling, i.e., the results of the aerody‐ namic 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

In an attempt to calibrate the domain size, mesh size, turbulence model and transition mod‐

This airfoil has been chosen as experimental results and results from other sources are avail‐ able for comparison. The experimental results have been obtained at the Low Speed Labora‐ tory 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 character‐

ized by Ramsey, gives the characteristics of the airfoil for angles of incidence ranging from

numbers of 0.75 to 1.5 million. Moreover, this study provides experimental results for the

) at different frequencies for Reynolds numbers between 0.75 and 1.4 million.

. The experiments were conducted on a 0.457 meters chord lenght for Reynolds

to 200

, 140 and 200

. The second study [32], real‐

) oscillating at (±5.50

and ±

flow.

100 Advances in Wind Power

**Figure 13.** S809 airfoil


100

to 400

*4.1.1. Model and convergence studies*

*4.1.1.1. Model and experimental results*

el, an S809 profile, designed by NREL, was used.

istics of the S809 profile for angles of incidence from -200

study of the dynamic stall for incidence angles of (80

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.

#### *4.1.1.3. Computational domain*

The computational domain is defined by a semi-disc of radius I1×c around the airfoil and two rectangles in the wake, of length I2×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 be‐ tween the domain boundary and the airfoil affects the results, only I1 and I2 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 prelimina‐ ry arbitrary one by a homothetic transformation with respect to the centre G and a factor b. Fig‐ ure 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.

**Figure 14.** Shape of the calculation domain

Figure 15 below respectively illustrates the drag and lift coefficients as a function of the ho‐ mothetic factor b for different angles of attack.


**Table 1.** Description of the trials through homothetic transformation

**Figure 15.** Drag and lift coefficients vs. homothetic factor for different angles of attack

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 ho‐ mothetic 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.20 . The evolution of the coefficients towards stabilization illustrate an important physical phe‐ nomenon: 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.

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

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

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 com‐ puted according to different AoA for a given Reynolds number and the results were com‐ pared 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.

## *4.1.1.5. Turbulence model calibration*

**Figure 15.** Drag and lift coefficients vs. homothetic factor for different angles of attack

9,597 m and width 4.799 m was used.

*4.1.1.4. Meshing*

102 Advances in Wind Power

**Table 2.** Mesh parameters

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 ho‐ mothetic 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.20

The evolution of the coefficients towards stabilization illustrate an important physical phe‐ nomenon: 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

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.

Figure 16 gives us an appreciation of the mesh we have used of in our simulations:

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

.

CFX proposes several turbulence models for resolution of flow over airfoils. Scientific litera‐ ture makes it clear that different turbulence models perform differently in different applica‐ tions. 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 suc‐ cessfully 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 Wil‐ cox 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 ob‐ tained with the different models were compared with the experimental OSU and DUT re‐ sults. D'Hamonville et al. [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.

#### *4.1.1.6. Transition model*

ANSYS-CFX proposes in the advanced turbulence control options several transitional mod‐ els 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 de‐ fine the onset of turbulence, referring to [33], we have only assessed the relative perform‐ ance 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 intermitten‐ cy 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 da‐ ta. 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 dif‐ fer from 8.20 . 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. Anal‐ ysis of the two figures brings us to the conclusion that the model with γ=0.94 provides re‐ sults very close to the DUT results. Therefore, we will compare the intermittency model with γ=0.94 with the other transitional models.

**Figure 17.** Drag and lift coefficients for different AoA using different intermittency values

Figure 18 illustrates the drag and lift coefficients for different AoA using different transition‐ al models.

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

Figure 18 shows that the drag coefficients for the three models are very close until 180 after which the results become clearly distinguishable. As from 20<sup>0</sup> , the γ-θ model over predicts the experimental values whereas such phenomena appear only after 22.10 for the two other models. For the lift coefficients, Figure 18 shows that the k-ω SST intermittency models pro‐ vide results closest to the experimental values for angles smaller than 14<sup>0</sup> . The k-ω SST mod‐ el under predicts the lift coefficients for angles ranging from 6<sup>0</sup> to 14<sup>0</sup> . For angles exceeding 200 , 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<sup>0</sup> . However, for AoA greater than 20<sup>0</sup> , a purely turbulent model needs to be used.
