**2. Fluid--Structure Interaction Model**

#### **2.1. Structural model**

Moreover, while the rest of the wind turbine subsystems are highly developed technological products, the blades are unique. There is no other technological application that uses such a device, so practical experience in blade manufacturing is relatively new. Blades also operate under a complex combination of fluctuating loads, and huge size differences complicate ex‐ trapolation of experimental data from the wind-tunnel to the prototype scale. Hence, com‐ puter models of fluid-structure interaction phenomena are particularly relevant to the design and optimization of wind-turbines. The wind-turbine industry is increasingly using computer models for blade structural design and for the optimization of its aerodynamics. Nevertheless, several features of the complex interaction of physical processes that charac‐ terize the coupled aeroelastic problem still exceed the capacities of existing commercial sim‐ ulation codes. Changes in structural response due to the development of new techniques in blade construction and/or the use of new materials would also represent a major factor to

take into account if the development of a new prototype blade is considered.

present fluid-structure solvers.

the blades.

124 Advances in Wind Power

Hence, a key factor for a breakthrough in wind-turbine technology is to reduce the uncer‐ tainties related to blade dynamics, by the improvement of the quality of numerical simula‐ tions of the fluid-structure interaction process, and by a better understanding of the underlying physics. The current state-of-the-art is to solve the aeroelastic equations in a fully non-linear coupled mode using Bernoulli or Timoshenko beam models (see [11], where a thorough coverage of the topic is presented). The goal is to provide the industry with a tool that helps them to introduce new technological solutions to improve the economics of blade design, manufacturing and transport logistics, without compromising reliability. A funda‐ mental step in that direction is the implementation of structural models capable of capturing the complex features of innovative prototype blades, so they can be tested at realistic fullscale conditions with a reasonable computational cost. To this end, we developed a *general‐ ized Timoshenko* code [27] based on a modified implementation the Variational-Asymptotic Beam Sectional technique (VABS) proposed by Hodges et al. (see [13] and references there‐ in). The ultimate goal is to combine this code with an advanced non-linear adaptive model of the unsteady flow, based on the vorticity-velocity formulation of the Navier-Stokes equa‐ tions, called the KLE model [32,33], which would offer performance advantages over the

In this chapter we present a set of tools for the design and full-scale analysis of the dynam‐ ics of composite laminate wind-turbine blades. The geometric design is carried on by means of a novel interpolation technique and the behavior of the blades is then simulated under normal operational conditions. We obtained results for the displacements and rotations of the blade sections along the span, section stresses, and fundamental vibrational modes of Even though a wind turbine blade is a slender structure that may be studied as a beam, they are usually not simple to model due to the inhomogeneous distribution of material properties and the complexity of their cross section (see Fig. 2. The *ad hoc* kinematic assump‐ tions made in classical theories (like the Bernoulli or the standard Timoshenko ap‐ proaches) may introduce significant errors, especially when the blade is vibrating with a wavelength shorter than its length. Complex blade geometry due to reasons of aerodynamic/ mechanical design, new techniques of blade construction, and the use of new materials combine themselves to give a new dimension to the problem. In order t0 obtain a fluidstructure interaction model capable of dealing with the complex features of new-genera‐ tion blades, we developed a code [27] based on a modified implementation of the Variational-Asymptotic Beam Sectional (VABS) model. Proposed and developed by Prof. Hodges and his collaborators [see 13, 39, and references therein], VABS is a model for curved and twist‐ ed composite beams that uses the same variables as classical Timoshenko beam theory, but the hypothesis of beam sections remaining planar after deformation is abandoned. In‐ stead, the real warping of the deformed section is interpolated by a 2-D finite-element mesh and its contribution to the strain energy is put in terms of the classical 1-D Timoshenko''s variables by means of a pre-resolution. The geometrical complexity of the blade section and/ or its material inhomogeneousness are reduced into a stiffness matrix for the 1-D beam. The reduced 1-D strain energy is equivalent to the actual 3-D strain energy in an asymptotic sense. Elimination of the *ad hoc* kinematic assumptions produces a fully populated *6 × 6* symmetric matrix for the 1-D beam, with as many as 21 stiffnesses, instead of the six fun‐ damental stiffnesses of the original Timoshenko theory. That is why VABS is referred to as a *generalized Timoshenko theory*.

**Figure 2.** Example of blade-section structural architecture representative of current commercial blade designs. The primary structural member is a box-spar, with a substantial build-up of spar cap material between the webs. The exte‐ rior skins and internal shear webs are both sandwich construction with triaxial fiberglass laminate separated by balsa core (from [9]).

Even for the case of large displacements and rotations of the beam sections, our model al‐ lows for accurate modeling of the bending and transverse shear in two directions, extension and torsion of the blade structure as a 1-D finite-element problem. Thus, this way we are able to decouple a general 3-D nonlinear anisotropic elasticity problem into a linear, 2-D, cross-sectional analysis (that may be solved *a priori*), and a nonlinear, 1-D, beam analysis for the global problem, which is what we would solve at each time step of a fluid-structure in‐ teraction analysis. This translates into substantial savings in computational cost as the struc‐ tural problem is solved along many timesteps. The cross-sectional 2-D analysis (that may be performed in parallel for the many cross sections along the blade) calculates the 3-D warp‐ ing functions asymptotically and finds the constitutive model for the 1-D nonlinear beam analysis of the blade. After one obtains the global deformation from the 1-D beam analysis, the original 3-D fields (displacements, stresses, and strains) can be recovered *a posteriori* us‐ ing the already-calculated 3-D warping functions.

**Figure 3.** Generalized Timoshenko theory: Schematic of the reference line, orthogonal triads, and beam sections be‐ fore and after deformation (adapted from [39])

In order to make this chapter self-contained, we shall see a brief outline of the theoretical basis of the dimensional reduction technique. More details can be found in [27, 13, 39] and references therein. Referring to Fig. 3, we have a reference line *R* drawn along the axis of the beam in the undeformed configuration. *R* could be twisted and/or curved according to the initial geometry of the beam. Section planes are normal to *R* at every point along its length. At the point where *R* intersects the section, an associated orthogonal triad , is de‐ fined in such a way that is tangent to *R* and are contained into the section plane; with a correspondent coordinate system *(X1 ,X2 ,X3 )* where *X1* is the coordinate along *R* and *X2 ,X3* are the Cartesian coordinates on the section plane. The position of a generic point on each section may be written as

$$\underline{\mathbf{R}}(X^i) = \underline{\mathbf{R}}(X^1) + X^\kappa \underline{\mathbf{B}}\_a(X^1) \tag{1}$$

Even for the case of large displacements and rotations of the beam sections, our model al‐ lows for accurate modeling of the bending and transverse shear in two directions, extension and torsion of the blade structure as a 1-D finite-element problem. Thus, this way we are able to decouple a general 3-D nonlinear anisotropic elasticity problem into a linear, 2-D, cross-sectional analysis (that may be solved *a priori*), and a nonlinear, 1-D, beam analysis for the global problem, which is what we would solve at each time step of a fluid-structure in‐ teraction analysis. This translates into substantial savings in computational cost as the struc‐ tural problem is solved along many timesteps. The cross-sectional 2-D analysis (that may be performed in parallel for the many cross sections along the blade) calculates the 3-D warp‐ ing functions asymptotically and finds the constitutive model for the 1-D nonlinear beam analysis of the blade. After one obtains the global deformation from the 1-D beam analysis, the original 3-D fields (displacements, stresses, and strains) can be recovered *a posteriori* us‐

**Figure 3.** Generalized Timoshenko theory: Schematic of the reference line, orthogonal triads, and beam sections be‐

In order to make this chapter self-contained, we shall see a brief outline of the theoretical basis of the dimensional reduction technique. More details can be found in [27, 13, 39] and references therein. Referring to Fig. 3, we have a reference line *R* drawn along the axis of the beam in the undeformed configuration. *R* could be twisted and/or curved according to the initial geometry of the beam. Section planes are normal to *R* at every point along its length. At the point where *R* intersects the section, an associated orthogonal triad , is de‐ fined in such a way that is tangent to *R* and are contained into the section plane;

> *,X2 ,X3*

are the Cartesian coordinates on the section plane. The position of a generic point on

*)* where *X1*

is the coordinate along *R* and

ing the already-calculated 3-D warping functions.

126 Advances in Wind Power

fore and after deformation (adapted from [39])

with a correspondent coordinate system *(X1*

each section may be written as

*X2 ,X3* where denotes the position of the center of the tern along *R*, and the index *α* assumes val‐ ues 2 and 3.

When the structure is deformed due to loading, the original reference line *R* adopts a new geometry *r*, and we have a new triad associated to each point, where is tangent to *r* and are contained into the normal plane. The material point whose original position was given by has now the position vector

$$\underline{\mathbf{r}}(X^l) = \underline{\mathbf{R}} + \underline{\mathbf{u}} + X^\mu \underline{\mathbf{t}}\_n + \underline{\mathbf{v}}\_l \underline{\mathbf{t}}\_{l'} \tag{2}$$

where *wi* are the contribution to the displacement due to warping. Now, we are able to com‐ pute the components of the gradient-of-deformation tensor as , where are respectively the covariant base vectors for the deformed configuration and the contravariant base vectors in the undeformed configuration, obtained from the kinematic description of equations (1) and (2). The Jaumann-Biot-Cauchy strain tensor is <sup>Γ</sup>*ij* <sup>=</sup> <sup>1</sup> <sup>2</sup> (*Fij* + *F ji* ) −*δij* , which provides a suitable measure of the 3-D strain field in terms of the beam strain measures and arbitrary warping functions. Γ is then used to compute the strain energy density function as

$$\text{22I} = \langle \{ \Gamma^T \mid \mathbf{S} \cdot \Gamma \} \rangle,\tag{3}$$

were, *S* is the matrix of the characteristics of the material expressed in the coordinates, and · =*∫ <sup>S</sup>* · *Gd <sup>X</sup>* <sup>2</sup> *d X* <sup>3</sup> , where *s* indicates the 2-D domain of the cross-section.

The next step is to find a strain energy expression asymptotically correct up to the second order of *h/l* and *h/R0*, where *h* is the characteristic size of the section, *l* the characteristic wavelength of deformation along the beam axis, and *R0* the characteristic radii of initial cur‐ vatures and twist of the beam. A complete second-order strain energy is sufficient for the purpose of constructing a generalized Timoshenko model because it is generally accepted that the transverse shear strain measures are one order less than classical beam strain meas‐ ures (extension, torsion and bending in two directions) [38]. A strain energy expression that asymptotically approximates the 3-D energy up to the second order is achieved using the Variational Asymptotic Method proposed in [4]. The complete derivation of this procedure is presented in [13], resulting in the following expression for the asymptotically correct strain energy:

$$\mathbf{2UI = \varepsilon \, ^T A \varepsilon + \varepsilon \, ^T \mathbf{2B} \varepsilon \, ^' + \varepsilon \, ^' \mathbf{C} \, \varepsilon \, ^' + \varepsilon \, ^T \mathbf{2D} \varepsilon \, ^'')}\tag{4}$$

where **A**, **B**, **C**, **D** are matrices that carry information on both the geometry and the material properties of the cross section, *()''* indicates the partial derivative with respect to the axial coordinate *X1* , and *<sup>ε</sup>* <sup>=</sup> *<sup>γ</sup>*¯ <sup>11</sup> *<sup>κ</sup>*¯ <sup>1</sup> *<sup>κ</sup>*¯ <sup>2</sup> *<sup>κ</sup>*¯ <sup>3</sup> *<sup>T</sup>* , are the strain measures defined in the classical Ber‐ noulli beam theory: *γ*¯ <sup>11</sup> the extension of the beam reference line, *<sup>k</sup>* ¯ <sup>1</sup> its torsion, and *k* ¯ <sup>2</sup> and *k* ¯ 3 the bending of the reference line in axes 2 and 3 due to the deformation.

The variational asymptotic procedure to get the matrices in equation (4) involves the discre‐ tization by finite-element techniques of the warping functions *wi* defined in expression (2).

During this procedure, a set of four constraints must be applied on *wi* . These restrictions, de‐ fined as *wi* =0 and *<sup>X</sup>* <sup>2</sup> *<sup>w</sup>*<sup>3</sup> <sup>−</sup> *<sup>X</sup>* <sup>3</sup> *w*<sup>2</sup> =0, where ∙ =*∫ <sup>s</sup>* <sup>∙</sup> *<sup>d</sup> <sup>X</sup>* <sup>2</sup> *<sup>d</sup> <sup>X</sup>* <sup>3</sup> , are intended to eliminate four rigid modes of displacement of the warped section (i.e. the three linear displacements plus the turn around ), which are already included in the Bernoulli strain measures *ε*. Pre‐ vious implementations of VABS (e.g. [39,13]) use the technique described by Cesnik *et al.* [6] to impose these constraints. In Cesnik et al.'s method, the rigid modes of displacement are suppressed explicitly. Then, the eigenvectors associated with the rigid modes in the matrix of the linear system that needs to be solved are computed, and used to get a reduced system. Instead of that, in our implementation of VABS, we use the Lagrangian-multiplier technique in its classical way to impose the constraints, solving the expanded system for the constrain‐ ed variational formulation itself. This simplifies the procedure by basically combining the whole solution in a single step. This simplification produces by itself a certain reduction in the overall computational cost, but most important, it has the advantage of allowing the use of the internal-node condensation technique in the finite-element discretization. As we shall see later, internal-node condensation allows us to substantially improve the efficiency of our solution by the tri-quadrilateral finite-element technique.

Expression (4) for the strain energy is asymptotically correct. Nevertheless, it is difficult to use in practice because it contains derivatives of the classical strain measures, which re‐ quires complicated boundary conditions. But, the well known Timoshenko beam theory is free from such drawbacks. Hence, the next step is to fit the strain energy in (4), into a gener‐ alized Timoshenko model of the form

$$2\mathbf{U} = \begin{bmatrix} \varepsilon^T & \boldsymbol{\gamma}\_s^T \end{bmatrix} \begin{bmatrix} \mathbf{X} & \mathbf{Y} \\ \mathbf{Y}^T & \mathbf{G} \end{bmatrix} \begin{bmatrix} \varepsilon \\ \boldsymbol{\gamma}\_s \end{bmatrix} = \boldsymbol{\varepsilon}^T \mathbf{X} \boldsymbol{\varepsilon} + 2\boldsymbol{\varepsilon}^T \mathbf{Y} \boldsymbol{\gamma}\_s + \boldsymbol{\gamma}\_s^T \mathbf{G} \boldsymbol{\gamma}\_s \tag{5}$$

where = *γ*<sup>11</sup> *κ*<sup>1</sup> *κ*<sup>2</sup> *κ*<sup>3</sup> *<sup>T</sup>* are the classical Timoshenko strain measures due to extension, torsion and bending, and *γ* = 2*γ*<sup>12</sup> 2*γ*<sup>13</sup> *<sup>T</sup>* the transverse shear strains.

What we need to find is **X**, **Y** and **G** in such a way that the strain energy in (4) and (5) would be equivalent up to at least second order. There is an identity that relates both the Bernoulli and the Timoshenko measures of deformation

Structural Analysis of Complex Wind Turbine Blades: Flexo-Torsional Vibrational Modes http://dx.doi.org/10.5772/51142 129

$$\varepsilon = \varepsilon\_\* + \mathbf{Q}\_\gamma \mathbf{y}^\gamma \mathbf{s}^\prime + \mathbf{P}\_\gamma \mathbf{y}\_{s\prime} \tag{6}$$

where

2*U* =*ε <sup>T</sup> Aε* + *ε <sup>T</sup>* 2*Bε* ′ + *ε* ′*TCε* ′ + *ε <sup>T</sup>* 2*Dε* ′′

coordinate *X1*

128 Advances in Wind Power

fined as *wi* =0 and *<sup>X</sup>* <sup>2</sup>

, and *<sup>ε</sup>* <sup>=</sup> *<sup>γ</sup>*¯ <sup>11</sup> *<sup>κ</sup>*¯ <sup>1</sup> *<sup>κ</sup>*¯ <sup>2</sup> *<sup>κ</sup>*¯ <sup>3</sup>

noulli beam theory: *γ*¯ <sup>11</sup> the extension of the beam reference line, *<sup>k</sup>*

tization by finite-element techniques of the warping functions *wi*

*<sup>w</sup>*<sup>3</sup> <sup>−</sup> *<sup>X</sup>* <sup>3</sup>

solution by the tri-quadrilateral finite-element technique.

*<sup>T</sup>* **X Y Y***<sup>T</sup>* **G**

ε γ*s*

alized Timoshenko model of the form

2*U* = ε*<sup>T</sup>* γ*<sup>s</sup>*

torsion and bending, and *γ* = 2*γ*<sup>12</sup> 2*γ*<sup>13</sup>

and the Timoshenko measures of deformation

where = *γ*<sup>11</sup> *κ*<sup>1</sup> *κ*<sup>2</sup> *κ*<sup>3</sup>

During this procedure, a set of four constraints must be applied on *wi*

the bending of the reference line in axes 2 and 3 due to the deformation.

where **A**, **B**, **C**, **D** are matrices that carry information on both the geometry and the material properties of the cross section, *()''* indicates the partial derivative with respect to the axial

The variational asymptotic procedure to get the matrices in equation (4) involves the discre‐

four rigid modes of displacement of the warped section (i.e. the three linear displacements plus the turn around ), which are already included in the Bernoulli strain measures *ε*. Pre‐ vious implementations of VABS (e.g. [39,13]) use the technique described by Cesnik *et al.* [6] to impose these constraints. In Cesnik et al.'s method, the rigid modes of displacement are suppressed explicitly. Then, the eigenvectors associated with the rigid modes in the matrix of the linear system that needs to be solved are computed, and used to get a reduced system. Instead of that, in our implementation of VABS, we use the Lagrangian-multiplier technique in its classical way to impose the constraints, solving the expanded system for the constrain‐ ed variational formulation itself. This simplifies the procedure by basically combining the whole solution in a single step. This simplification produces by itself a certain reduction in the overall computational cost, but most important, it has the advantage of allowing the use of the internal-node condensation technique in the finite-element discretization. As we shall see later, internal-node condensation allows us to substantially improve the efficiency of our

Expression (4) for the strain energy is asymptotically correct. Nevertheless, it is difficult to use in practice because it contains derivatives of the classical strain measures, which re‐ quires complicated boundary conditions. But, the well known Timoshenko beam theory is free from such drawbacks. Hence, the next step is to fit the strain energy in (4), into a gener‐

What we need to find is **X**, **Y** and **G** in such a way that the strain energy in (4) and (5) would be equivalent up to at least second order. There is an identity that relates both the Bernoulli

<sup>=</sup>ε*<sup>T</sup>* **<sup>X</sup>**<sup>ε</sup> <sup>+</sup> <sup>2</sup>ε*<sup>T</sup>* **<sup>Y</sup>**γ*<sup>s</sup>* <sup>+</sup> <sup>γ</sup>*<sup>s</sup>*

*<sup>T</sup>* are the classical Timoshenko strain measures due to extension,

*<sup>T</sup>* the transverse shear strains.

*w*<sup>2</sup> =0, where ∙ =*∫*

, (4)

<sup>1</sup> its torsion, and *k*

defined in expression (2).

. These restrictions, de‐

, are intended to eliminate

*<sup>T</sup>* **<sup>G</sup>**γ*<sup>s</sup>* , (5)

¯ <sup>2</sup> and *k* ¯ 3

*<sup>T</sup>* , are the strain measures defined in the classical Ber‐

*<sup>s</sup>* <sup>∙</sup> *<sup>d</sup> <sup>X</sup>* <sup>2</sup> *<sup>d</sup> <sup>X</sup>* <sup>3</sup>

¯

$$\mathbf{Q}\_{\boldsymbol{\gamma}}^{T} = \begin{bmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & -1 & 0 \end{bmatrix} \qquad \qquad \mathbf{y} \qquad \qquad \mathbf{P}\_{\boldsymbol{\gamma}}^{T} = \begin{bmatrix} 0 & K\_{2} & -K\_{1} & 0 \\ 0 & K\_{3} & 0 & -K\_{1} \end{bmatrix} \tag{7}$$

being *K1* the twist, and *K2* and *K3* the curvatures of the original reference line *R*. Thus, using (6), we may rewrite expression (4) in terms of the generalized Timoshenko strain meas‐ ures using the 1-D equilibrium equations. This provides a way to relate the derivatives of strain measures with the strain measures themselves, to fit the resulting expression into the generalized Timoshenko form (5). Then, an asymptotic method is used to get approxima‐ tions to **X**, **Y** and **G**; using as input the already computed matrices **A**, **B**, **C**, **D** (see [39] for details). Finally, a stiffness matrix for the 1-D beam problem is formed as a simple reor‐ dering of the matrix **X Y <sup>Y</sup>***<sup>T</sup>* **<sup>G</sup>** }, in such a way as to get a functional for the strain energy density of expression~(5)

$$2U = \bar{\gamma}^T \boxtimes \bar{\gamma} \tag{8}$$

where *γ*¯ <sup>=</sup> *<sup>γ</sup> <sup>κ</sup>* is the array of Timoshenko measures of deformation regrouped in a more con‐ venient way, *γ <sup>T</sup>* = *γ*<sup>11</sup> 2*γ*<sup>12</sup> 2*γ*13 and *κ <sup>T</sup>* = *κ*<sup>1</sup> *κ*<sup>2</sup> *κ*<sup>3</sup> .

For the discretization of the 2-D sections, we adopted the tri-quadrilateral finite-element technique, which is based on the use of nine-node biquadratic isoparametric finite elements that possess a high convergence rate and, due their biquadratic interpolation of the geomet‐ ric coordinates, provide the additional ability of reducing the so-called skin-error on curvi‐ linear boundaries when compared to linear elements. For a detailed description of the isoparametric-element technique and its corresponding interpolation functions see Bathe [3].

In order to combine the advantages of the nine-node quadrilateral isoparametric element with the geometrical ability of a triangular grid to create suitable non-structured meshes with gradual and smooth changes of mesh density, we implemented what we called triquadrilateral isoparametric elements. The tri-quadrilateral elements consist of an assembling of three quadrilateral nine-node isoparametric elements in which each triangle of a standard unstructured mesh is divided into. By static condensation of the nodes that lie inside the tri‐ angle, we can significantly reduce the number of nodes to solve in the final system, subse‐ quently recovering the values for the internal nodes from the solution on the noncondensable nodes. The internal nodes may be expressed in terms of nodes which lay on the elemental boundary following the classical procedure for elemental condensation (see [3]). This process of condensation allows us to reduce the size of the new system to solve to ap‐ proximately a 40% of the original system. The use of the static condensation procedure is attractive not only because it reduces the size of the stiffness matrices arising in finite-ele‐ ment and spectral-element methods but also because it improves the condition number of the final condensed system. This is related with the properties of the Schur-complement technique. The condensed system is essentially the Schur complement of the interior-node submatrix in the non-condensed original system. A detailed description of the tri-quadrilat‐ eral technique may be found in Ponta [33]; including a schematic example of a mesh of triquadrilateral finite elements obtained from the original triangular discretization, and a description of the internal topology of the tri-quadrilateral element.

To solve the one-dimensional problem for the equivalent beam, we use a formulation based on the intrinsic equations for the beam obtained from variational principles [12], and weighted in an energy-consistent way according to Patil et al. [30], which produces the fol‐ lowing variational formulation:

$$\begin{aligned} \int\_{0}^{\ell} \underbrace{\left[\boldsymbol{\upxi}\mathbf{V}^{T}\,\mathbf{I}\,\dot{\mathbf{V}} + \underbrace{\boldsymbol{\upxi}\mathbf{F}^{T}\,\mathbf{S}^{-1}\,\dot{\mathbf{F}}\right]}\_{2} d\mathbf{X}^{1} &= \int\_{0}^{\ell} \left[\underbrace{\boldsymbol{\upxi}\mathbf{V}^{T}\,\mathbf{F}^{\prime}}\_{3} + \underbrace{\boldsymbol{\upxi}\mathbf{V}^{T}\,\dot{\mathbf{K}}\,\mathbf{F}}\_{4} + \underbrace{\boldsymbol{\upxi}\mathbf{V}^{T}\,\dot{\mathbf{V}}\,\mathbf{F}}\_{5} + \underbrace{\boldsymbol{\upxi}\mathbf{V}^{T}\,\dot{\mathbf{V}}\,\mathbf{F}}\_{6} + \underbrace{\boldsymbol{\upxi}\mathbf{V}^{T}\,\dot{\mathbf{V}}\,\mathbf{V}}\_{9} \right] d\mathbf{X}^{1}, \end{aligned} \tag{9}$$

where

$$\begin{aligned} \mathbf{F} &= \begin{bmatrix} \mathbf{F} \\ \mathbf{M} \end{bmatrix} & \stackrel{\scriptstyle \mathbf{V}}{\mathbf{V}} &= \begin{bmatrix} \mathbf{V} \\ \boldsymbol{\Omega} \end{bmatrix} & \stackrel{\scriptstyle \mathbf{f}}{\mathbf{f}} &= \begin{bmatrix} \mathbf{f} \\ \mathbf{m} \end{bmatrix} \\ \stackrel{\scriptstyle \mathbf{\tilde{V}}}{\underset{\scriptstyle \mathbf{V}}{\mathbf{\tilde{K}}}} & \stackrel{\scriptstyle \mathbf{\tilde{V}}}{\stackrel{\scriptstyle \mathbf{\tilde{K}}}{\mathbf{\tilde{K}}}} & \stackrel{\scriptstyle \mathbf{\tilde{V}}}{\mathbf{\tilde{V}}} & \stackrel{\scriptstyle \mathbf{\tilde{K}}}{\mathbf{\tilde{K}}} & \stackrel{\scriptstyle \mathbf{\tilde{K}}}{\mathbf{\tilde{K}}} \end{bmatrix} \end{aligned}$$

Tilde indicates the skew-symmetric matrix associated to a vector magnitude in such a way that, for example, if we have any pair of vectors **A** and **B**, the matrix--vector product **A**˜ **B** is equivalent to the cross product **A**× **B**. Thus, *γ*˜ is associated with *κ*˜ with *κ* , **V**˜ with **V**, and so forth. Hence, matrix *γ*˜ is a rearrangement of the components of the strain-measures vector *γ*¯ defined above, the generalized-velocities vector **V**¯ and matrix **V**˜ represent the components of the linear and angular velocities, and matrix **K**˜ represents the initial torsion and curva‐ tures of the beam (matrix **e**<sup>1</sup> ˜ is the skew-symmetric matrix associated to **e**<sup>1</sup> *<sup>T</sup>* = 1 0 0 , the unit vector along *X1* ). The generalized-forces vector **F**¯ represents the forces and moments re‐ lated with the strain measures , and the generalized-distributed-loads vector **f** ¯ represents the forces and moments distributed along the axis of the beam. Here, is the same stiffness matrix for the 1-D model, see equation (8); and **I** ¯ is the inertia matrix of each section. The upper dot indicates a time derivative, and the prime a derivative with respect to the longitudinal coordinate of the beam *X1* .

This variational formulation was discretized by the spectral-element method (see [20,29]). The magnitudes in (9) where replaced by their interpolated counterparts: **V**¯ <sup>=</sup>**HV**¯ *<sup>e</sup>* **Q***<sup>e</sup>* , and **<sup>F</sup>**¯ <sup>=</sup>**HF**¯ *<sup>e</sup>* **Q***<sup>e</sup>* , where **HV**¯ *<sup>e</sup>* and **HF**¯ *e* are the interpolation-function arrays, and **Q***<sup>e</sup>* is a vector contain‐ ing the nodal values of both the generalized velocities and the generalized forces. Super‐ script *e* indicates discretization of the terms at the elemental level, which will disappear after the final assembly of the terms into the global matrix for the whole beam. The axial deriva‐ tives of the magnitudes were interpolated in a similar way: **V**¯′ =**BV**¯ *<sup>e</sup>* **Q***<sup>e</sup>* , and **F**¯′ =**BF**¯ *<sup>e</sup>* **Q***<sup>e</sup>* , where **BV**¯ *<sup>e</sup>* and **BF**¯ *e* are the arrays for the interpolation-function derivatives. Then, we arrived to the discretized version of (9):

$$
\delta \mathbf{Q}^{\epsilon \tau} \mathbf{M}\_1^{\epsilon} \mathbf{Q}^{\epsilon} = \delta \mathbf{Q}^{\epsilon \tau} (\mathbf{K}\_1^{\epsilon} + \mathbf{K}\_2^{\epsilon}) \mathbf{Q}^{\epsilon} + \delta \mathbf{Q}^{\epsilon \tau} \mathbf{K}\_q^{\epsilon} \mathbf{q}^{-\epsilon} + \delta \mathbf{Q}^{\epsilon \tau} \mathbf{B}\_Q^{\epsilon} \mathbf{Q}^{\epsilon}), \tag{10}$$

where

This process of condensation allows us to reduce the size of the new system to solve to ap‐ proximately a 40% of the original system. The use of the static condensation procedure is attractive not only because it reduces the size of the stiffness matrices arising in finite-ele‐ ment and spectral-element methods but also because it improves the condition number of the final condensed system. This is related with the properties of the Schur-complement technique. The condensed system is essentially the Schur complement of the interior-node submatrix in the non-condensed original system. A detailed description of the tri-quadrilat‐ eral technique may be found in Ponta [33]; including a schematic example of a mesh of triquadrilateral finite elements obtained from the original triangular discretization, and a

To solve the one-dimensional problem for the equivalent beam, we use a formulation based on the intrinsic equations for the beam obtained from variational principles [12], and weighted in an energy-consistent way according to Patil et al. [30], which produces the fol‐

description of the internal topology of the tri-quadrilateral element.

**<sup>F</sup>**<sup>=</sup> **<sup>F</sup>**

**γ ^** <sup>=</sup> **<sup>κ</sup>**˜ **<sup>0</sup>** *γ* ^ **<sup>κ</sup>˜** , **<sup>V</sup>**

same stiffness matrix for the 1-D model, see equation (8); and **I**

**<sup>M</sup>** , **<sup>V</sup>**¯ <sup>=</sup> **<sup>V</sup>**

**<sup>Ω</sup>** , **<sup>f</sup>** ¯ = **f m** ,

**<sup>V</sup>**˜ **<sup>Ω</sup>˜** , **<sup>K</sup>**

˜ is the skew-symmetric matrix associated to **e**<sup>1</sup>

). The generalized-forces vector **F**¯ represents the forces and moments re‐

*<sup>T</sup>* = 1 0 0 , the

¯ is the inertia matrix of each

¯

^ <sup>=</sup> **<sup>K</sup>**˜ <sup>0</sup> **<sup>e</sup>**˜ <sup>1</sup> **<sup>K</sup>**˜ .

^ <sup>=</sup> **<sup>Ω</sup>**˜ <sup>0</sup>

Tilde indicates the skew-symmetric matrix associated to a vector magnitude in such a way that, for example, if we have any pair of vectors **A** and **B**, the matrix--vector product **A**˜ **B** is equivalent to the cross product **A**× **B**. Thus, *γ*˜ is associated with *κ*˜ with *κ* , **V**˜ with **V**, and so forth. Hence, matrix *γ*˜ is a rearrangement of the components of the strain-measures vector *γ*¯ defined above, the generalized-velocities vector **V**¯ and matrix **V**˜ represent the components of the linear and angular velocities, and matrix **K**˜ represents the initial torsion and curva‐

lated with the strain measures , and the generalized-distributed-loads vector **f**

represents the forces and moments distributed along the axis of the beam. Here, is the

section. The upper dot indicates a time derivative, and the prime a derivative with respect to

.

lowing variational formulation:

130 Advances in Wind Power

tures of the beam (matrix **e**<sup>1</sup>

the longitudinal coordinate of the beam *X1*

unit vector along *X1*

where

$$\begin{aligned} \mathbf{M}\_{\mathbf{1}}^{e} &= \int\_{-1}^{1} \left[ \mathbf{H}\_{\widehat{\mathbf{V}}}^{e} \,^{T} \, \mathbf{I} \, \mathbf{H}\_{\widehat{\mathbf{V}}}^{e} + \mathbf{H}\_{\widehat{\mathbf{F}}}^{e} \, \mathbf{S}^{-1} \, \mathbf{H}\_{\widehat{\mathbf{F}}}^{e} \right] J \, dt \\\\ \mathbf{K}\_{1}^{e} &= \int\_{-1}^{1} \left[ \mathbf{H}\_{\widehat{\mathbf{V}}}^{e} \, ^{T} \, \mathbf{B}\_{\widehat{\mathbf{F}}}^{e} + \mathbf{H}\_{\widehat{\mathbf{F}}}^{e} \, ^{T} \, \mathbf{B}\_{\widehat{\mathbf{V}}}^{e} \right] J \, dt \, \mathbf{Q}^{e}, \\\\ \mathbf{K}\_{2}^{e} &= \int\_{-1}^{1} \left[ \mathbf{H}\_{\widehat{\mathbf{V}}}^{e} \, ^{T} \hat{\mathbf{K}} \, \mathbf{H}\_{\widehat{\mathbf{F}}}^{e} - \mathbf{H}\_{\widehat{\mathbf{F}}}^{e} \, \hat{\mathbf{K}}^{T} \, \mathbf{H}\_{\widehat{\mathbf{V}}}^{e} \right] J \, dt \\\\ \mathbf{K}\_{2}^{e} &= \int\_{-1}^{1} \, \mathbf{H}\_{\widehat{\mathbf{V}}}^{e} \, ^{T} \, \mathbf{H}\_{\widehat{\mathbf{D}}}^{e} \, J \, dt, \end{aligned}$$

**M1** *e* corresponds to the discretization of terms 1 and 2 giving the equivalent of a mass matrix. **K1** *e* , corresponding to terms 3 and 8, is the stiffness matrix of the 1-D problem. **K2** *e* , corre‐ sponding to terms 4 and 9, is the additional stiffness related with the twist and curvature of the undeformed configuration. **Kq** *<sup>e</sup>* corresponds to the evaluation of term 6, the contribution of the distributed loads; and **q**¯*<sup>e</sup>* is an array containing the nodal values of the generalized distributed loads. *t* is the natural coordinate in the elements and *J* is the Jacobian of the map‐ ping from the problem coordinate *X1* to *t* (see [3]). The discretized version of the terms in (9) related to non-linear interactions, i.e. terms 5, 7 and 10, gives

$$\mathbf{B}\_Q^\varepsilon \left( \mathbf{Q}^\varepsilon \right) = \int\_{-1}^1 \left[ \mathbf{H}\_{\widehat{\mathbf{V}}}^{\varepsilon^T} \, \widehat{\boldsymbol{\gamma}} \, \mathbf{H}\_{\widehat{\mathbf{F}}}^{\varepsilon} - \mathbf{H}\_{\widehat{\mathbf{V}}}^{\varepsilon^T} \, \widehat{\mathbf{V}} \, \mathbb{I} \, \mathbf{H}\_{\widehat{\mathbf{V}}}^{\varepsilon} - \mathbf{H}\_{\widehat{\mathbf{F}}}^{\varepsilon^T} \, \widehat{\boldsymbol{\gamma}}^T \, \mathbf{H}\_{\widehat{\mathbf{V}}}^{\varepsilon} \right] \mathbf{Q}^{\varepsilon^\*} \, \mathbb{J} \, dt$$

A linearization of **B***<sup>Q</sup> <sup>e</sup>* (**Q***<sup>e</sup>*) around any given configuration **Q**<sup>1</sup> *e* gives the matrix

$$\begin{aligned} \mathbf{K}\_{\mathbf{N}}^{\varepsilon} \left( \mathbf{Q}\_{1}^{\varepsilon} \right) &= \int\_{-1}^{1} \left\{ \mathbf{H}\_{\mathbf{\hat{V}}}^{\varepsilon^{T}} \left[ \hat{\boldsymbol{\gamma}}\_{1} \ \mathbf{H}\_{\mathbf{\hat{F}}}^{\varepsilon} - \hat{\mathbf{V}}\_{1} \ \mathbf{I} \ \mathbf{H}\_{\mathbf{\hat{V}}}^{\varepsilon} - \hat{\mathbf{F}}\_{1} \ \mathbf{S}^{-1} \ \mathbf{H}\_{\mathbf{\hat{F}}}^{\varepsilon} + \hat{\mathbf{P}}\_{1} \ \mathbf{H}\_{\mathbf{\hat{V}}}^{\varepsilon} \right] + \cdots \\ & \quad \mathbf{H}\_{\mathbf{\hat{F}}}^{\varepsilon^{T}} \left[ \hat{\mathbf{V}}\_{1}^{T} \ \mathbf{S}^{-1} \ \mathbf{H}\_{\mathbf{\hat{F}}}^{\varepsilon} - \hat{\mathbf{P}}\_{1}^{T} \ \mathbf{H}\_{\mathbf{\hat{V}}}^{\varepsilon} \right] \end{aligned}$$

where

$$
\hat{\mathbf{F}} = \begin{bmatrix} \mathbf{0} & \mathbf{\tilde{F}} \\ \mathbf{\tilde{F}} & \mathbf{\tilde{M}} \end{bmatrix}, \quad \hat{\mathbf{P}} = \begin{bmatrix} \mathbf{0} & \mathbf{\tilde{P}}\_v \\ \mathbf{\tilde{P}}\_v & \mathbf{\tilde{P}}\_\omega \end{bmatrix}.
$$

Matrix **F**¯ is a rearrangement of the components of the generalized-forces vector **F**¯ defined above. Matrix **P**¯ is a rearrangement of the components of the generalized-momentum vector

**P**= **P***v* **P***ω* , which represents the linear and angular momenta related with the generalized-ve‐

locities . Tilde operates in the same way defined before, and the subscript 1 indi‐ cates the value of the magnitudes at a given state **Q**<sup>1</sup> *e* .

Finally, after the assembly of the elemental terms into the global system, the solution for the nonlinear problem (9) in its steady state was obtained by solving iteratively for Δ**Q** the dis‐ cretized expression

$$
\left[\mathbf{K}\_1 + \mathbf{K}\_2 + \mathbf{K}\_N \left(\mathbf{Q}^{(l)}\right)\right] \Delta \mathbf{Q} = -\mathbf{K}\_\mathbf{q} \mathbf{q} - \left(\mathbf{K}\_1 + \mathbf{K}\_2\right) \mathbf{Q}^{(l)} - \mathbf{B}\_\zeta \left(\mathbf{Q}^{(l)}\right), \tag{11}
$$

and updating the global vector of nodal values of the generalized velocities and forces as **Q**(*i*+1) =**Q**(*i*) + Δ**Q**.

From the steady-state solution we also obtained the vibrational modes of the blade structure and their corresponding frequencies by solving the eigenvalue problem

$$\mathbf{M}\_1 \dot{\mathbf{Q}} + \left[ \mathbf{K}\_1 + \mathbf{K}\_2 + \mathbf{K}\_N \{ \mathbf{Q}^{(i)} \} \right] \mathbf{Q} = \mathbf{0}.\tag{12}$$

From these results for the intrinsic equations we recovered the displacements and rotations of the blade sections by solving the kinematic equations for the beam (see [13])

$$\mathbf{u}' - \mathbf{C}\_{rR}^T(\gamma + \mathbf{e}\_1) + \mathbf{e}\_1 + \tilde{\mathbf{K}}\ \mathbf{u} = \mathbf{0},\tag{13}$$

Structural Analysis of Complex Wind Turbine Blades: Flexo-Torsional Vibrational Modes http://dx.doi.org/10.5772/51142 133

$$
\tilde{\mathbf{K}} + \tilde{\kappa} + \mathbf{C}'\_{rR} \mathbf{C}^T\_{rR} - \mathbf{C}\_{rR} \tilde{\mathbf{K}} \, \mathbf{C}^T\_{rR} = \mathbf{0},\tag{14}
$$

where **u** is the vector of displacements of each point along the reference line from its posi‐ tion in the reference configuration to the one in the deformed configuration, and **C***rR* is the orthogonal matrix that rotates the local triad from its original orientation in the reference configuration to the one in the deformed configuration (both are defined in function of the longitudinal coordinate *X1* ). The strains *γ* and *κ* were computed from the generalized forces and the stiffness of the corresponding blade section. Equations (13) and (14) were also line‐ arized, and like the other expressions, discretized by the spectral-element method.

#### **2.2. Aerodynamic model**

A linearization of **B***<sup>Q</sup>*

132 Advances in Wind Power

where

**P**= **P***v* **P***ω*

**Q**(*i*+1)

cretized expression

=**Q**(*i*) + Δ**Q**.

*<sup>e</sup>* (**Q***<sup>e</sup>*) around any given configuration **Q**<sup>1</sup>

Matrix **F**¯ is a rearrangement of the components of the generalized-forces vector **F**¯ defined above. Matrix **P**¯ is a rearrangement of the components of the generalized-momentum vector

locities . Tilde operates in the same way defined before, and the subscript 1 indi‐

Finally, after the assembly of the elemental terms into the global system, the solution for the nonlinear problem (9) in its steady state was obtained by solving iteratively for Δ**Q** the dis‐

and updating the global vector of nodal values of the generalized velocities and forces as

From the steady-state solution we also obtained the vibrational modes of the blade structure

From these results for the intrinsic equations we recovered the displacements and rotations

and their corresponding frequencies by solving the eigenvalue problem

**M1Q˙** <sup>+</sup> **K1** <sup>+</sup> **K2** <sup>+</sup> **KN**(**Q**(*i*)

of the blade sections by solving the kinematic equations for the beam (see [13])

cates the value of the magnitudes at a given state **Q**<sup>1</sup>

, which represents the linear and angular momenta related with the generalized-ve‐

*e* . *e*

gives the matrix

(11)

(13)

) **Q**=**0**. (12)

The flow model that interacts with the structural counterpart presented in section 2.1, called *Large Sectional Rotation BEM (LSR-BEM)*, is responsible to provide the aerodynamic loads along the rotor blades, and is sensitive enough to take into account all the complex deforma‐ tion modes that the structural model is able to solve. The basis for our aerodynamic model is the well known Blade Element Momentum theory (BEM). Nevertheless, due to the high lev‐ el of detail that our structural model can provide, a complete reformulation was needed in the aerodynamic model to get a compatible level of description.

The tendency in the wind-turbine industry to increase the size of the state-of-the-art ma‐ chine [17] drives not only to bigger, but also to more flexible blades which are relatively lighter. It is observed for this type of wind turbine blades that big deformations, either due to blade flexibility or to pre-conforming processes, produces high rotations of the blade sections. Moreover, blades could be pre-conformed with specific curvatures given to any of their axis *(i.e. conning/sweeping)*. This tendency puts in evidence one of the most important limitations of the current BEM theory. While the basics of this theory keeps being perfect‐ ly valid, the actual mathematical formulation implies the assumption of blade sections re‐ maining perpendicular to an outwards radial line contained in the plane of the actuator disk coincident with the rotor's plane. That is, even though the basics of the BEM theory (i.e. the equation of the aerodynamic loads and the change of momentum in the streamtubes) keeps being valid, the mathematical formulation cannot represent large rotations of the blade sections. This basically leads to a misrepresentation of the effects of the large deformation associated to flexible blades on the computation of the aerodynamic loads. Hence, a new mathematical formulation is required to project the velocities obtained from momentum theory onto the blade element's plane and then re-project backwards the resulting forces from Blade Element theory onto the plane of the stream tube actuator disk. When analyz‐ ing BEM theory for this cases, the principle of equating the forces obtained by Blade Ele‐ ment theory with the ones coming from the the changing of momentum in the stream tube is still valid.

In what follows we will describe the main characteristics of our model, and refer to [25] and [5] for details on the classical *BEM theory*.We start by defining a set of orthogonal matrices that perform the rotation of the physical magnitudes involved (velocities, forces, etc.) The interaction with external control modules will require a constant update of this projection matrices. For example, the rotor azimuth matrix, besides the instantaneous position of the blade along its rotation, can reflect control actions on the dynamics of the Electro-Mechani‐ cal train that define the rotor's angular speed, *Ω* .

For instance, we could write the wind velocity vector **Wh** facing the differential annulus of our actuator disk, affecting its components, according to BEM theory, by the axial induction factor *a* and the rotational induction factor *a'*. The *h* subscript here indicates that the wind velocity vector is described in the *hub* coordinate system according to standards from the In‐ ternational Electrotechnical Commission (IEC) [15].

$$\mathbf{W}\_h = \begin{bmatrix} \mathbf{W}\_{wh\_{\langle 1\rangle}} (1 - a) \\ \mathbf{W}\_{wh\_{\langle 2\rangle}} + \boldsymbol{\Omega} \, r\_h \, (1 + a^{'}) \\ \mathbf{W}\_{wh\_{\langle 3\rangle}} \end{bmatrix} \tag{15}$$

where **W***wh* is the incoming wind velocity projected into the *h* coordinate system, *Ω* is the angular velocity of the rotor and *rh* is the radial distance of the airfoil section in the *h* coordi‐ nate system.

Then, to compute the relative velocity affecting a blade element, we will project **Wh** going through the different coordinate systems, from the hub, until reaching the blade's section co‐ ordinate system. Let's see first which rotations we shall go through, and which matrices will transform our velocity vector from one coordinate system to the other.

Thus, the conning rotation matrix **C***θcn* is a linear operator with a basic rotation taking into account the conning angle for the rotor, and the Pitching rotation matrix **C***θ<sup>p</sup>* , represents a rotation around the pitch of the blade.

$$\mathbf{C}\_{\boldsymbol{\theta}\_{p}} = \begin{bmatrix} \cos \left( \boldsymbol{\theta}\_{p} \right) & -\sin \left( \boldsymbol{\theta}\_{p} \right) & 0 \\ \sin \left( \boldsymbol{\theta}\_{p} \right) & \cos \left( \boldsymbol{\theta}\_{p} \right) & 0 \\ 0 & 0 & 1 \end{bmatrix} \tag{16}$$

where *θp* is the pitch angle.

Two more re-orientations are needed in order to get to the instantaneous coordinate sys‐ tem of the blade sections, associated with the deformed reference line r of section 2.1. The first of this matrices contains information on blade section's geometry at the time the blade was designed and manufactured. As it was mentioned previously, the blade could have pre-conformed curvatures along its longitudinal axis *(i.e. the blade axis is no longer rectilin‐ ear)*. This curvatures can reflect either an initial twist along the longitudinal axis or a com‐ bination of twist plus pre-bending on the other two axes *(i.e. conning/sweeping)*. To this end, we compute during the blade design stage a set of transformation matrices which contain the information of the three dimensional orientation of the blade''s sections for each posi‐ tion on the longitudinal axis as we move along the span. To this end, we compute the Frenet-Serret formulas that define the curvature of the (now curvilinear) longitudinal axis. These differential formulas provide the means to describe the *tangent*, *normal* and *binormal* unit vectors on a given curve. Due to this unit vectors, the Frenet-Serret coordinate system is also known as the *TNB frame*. More information about the calculation of the TNB unit vectors, their properties and other applications can be found in [37]. Around the tangential axis of the TNB, there is a further rotation of each blade section to orient it accordingly to the particular twist specified on the blade's aerodynamic design. Combining these rotations we then create a transformation matrix for every blade section at different span positions. We call this matrix the **C***Rb*, as it relates the global coordinate system of the blade *b*, with the system of coordinates of the blade sections in the undeformed configuration defined by line *R*, as in section 2.1.

interaction with external control modules will require a constant update of this projection matrices. For example, the rotor azimuth matrix, besides the instantaneous position of the blade along its rotation, can reflect control actions on the dynamics of the Electro-Mechani‐

For instance, we could write the wind velocity vector **Wh** facing the differential annulus of our actuator disk, affecting its components, according to BEM theory, by the axial induction factor *a* and the rotational induction factor *a'*. The *h* subscript here indicates that the wind velocity vector is described in the *hub* coordinate system according to standards from the In‐

cal train that define the rotor's angular speed, *Ω* .

ternational Electrotechnical Commission (IEC) [15].

Thus, the conning rotation matrix **C***θcn*

rotation around the pitch of the blade.

where *θp* is the pitch angle.

nate system.

134 Advances in Wind Power

**W***<sup>h</sup>* =

*Wwh* (1)

*Wwh* (2)

*Wwh* (3)

transform our velocity vector from one coordinate system to the other.

**C***θ<sup>p</sup>* =

account the conning angle for the rotor, and the Pitching rotation matrix **C***θ<sup>p</sup>*

cos (*θp*) −sin (*θp*) 0 sin (*θp*) cos (*θp*) 0 0 0 1

Two more re-orientations are needed in order to get to the instantaneous coordinate sys‐ tem of the blade sections, associated with the deformed reference line r of section 2.1. The first of this matrices contains information on blade section's geometry at the time the blade was designed and manufactured. As it was mentioned previously, the blade could have pre-conformed curvatures along its longitudinal axis *(i.e. the blade axis is no longer rectilin‐ ear)*. This curvatures can reflect either an initial twist along the longitudinal axis or a com‐ bination of twist plus pre-bending on the other two axes *(i.e. conning/sweeping)*. To this end, we compute during the blade design stage a set of transformation matrices which contain the information of the three dimensional orientation of the blade''s sections for each posi‐

(1−*a*)

<sup>+</sup> <sup>Ω</sup> *rh* (1 <sup>+</sup> *<sup>a</sup>* ′

where **W***wh* is the incoming wind velocity projected into the *h* coordinate system, *Ω* is the angular velocity of the rotor and *rh* is the radial distance of the airfoil section in the *h* coordi‐

Then, to compute the relative velocity affecting a blade element, we will project **Wh** going through the different coordinate systems, from the hub, until reaching the blade's section co‐ ordinate system. Let's see first which rotations we shall go through, and which matrices will

)

, (15)

, (16)

, represents a

is a linear operator with a basic rotation taking into

After applying the **C***Rb*, one more projection is needed to get to the instantaneous coordinate system associated with *r*. This last transformation is given by the **C***rR* matrix, computed by the 1D structural model, see equation 14. It contains information to transform vectors from the *R* to the *r* systems after structural deformations had occurred. Note that this matrix is updated at every timestep of the 1D model during dynamic simulations, being one of the key variables transporting information between the structural and aerodynamic models.

After all these projections of the **Wh** vector, we have the relative wind velocity expressed in the blade''s section coordinate system. The expression for the flow velocity relative to the blade section, **W***rel* :

$$\mathbf{W}\_{rel} = \left(\mathbf{C}\_{rR}\mathbf{C}\_{Rb}\mathbf{C}\_{\boldsymbol{\theta}\_{p}}\mathbf{C}\_{\boldsymbol{\theta}\_{cn}}\mathbf{W}\_{h}\right) + \mathbf{v}\_{str} \tag{17}$$

where the addition of **v***str* corresponds to the blade section structural deformation velocities, coming from the structural model.

Then, the magnitude |**W***rel* | and the angle of attack *α* are used to compute the forces on the airfoil section through the aerodynamic coefficients *Cl* , *Cd*. Another innovation of our model is that the data tables from static wind-tunnel are corrected at each timestep to consider ei‐ ther rotational-augmentation or dynamic-stall effects, or both.

The aerodynamic loads acting on the blade element is then projected back onto the *h* coordi‐ nate system,

$$\mathbf{dF}\_h = \mathbf{C}\_{\partial\_{cn}}^T \mathbf{C}\_{\partial\_p}^T \mathbf{C}\_{Rb}^T \mathbf{C}\_{rR}^T \mathbf{C}\_{Lthal} \mathbf{dF}\_r \tag{18}$$

where **C***Lthal* is the matrix which projects the lift and drag forces onto the chord-wise and chord-normal directions, which are aligned with the coordinates of *r*. Finally, as in the clas‐ sical BEM theory, **dF***h* is equated to the rate of change of momentum in the annular stream‐ tube corresponding to the blade element. The component normal to the rotor's disk, is equated to the change in axial momentum, while the tangential component, is equated to the change of angular momentum.

In order to apply this theory to HAWT rotors, we must introduce some corrective factors into the calculation process. BEM theory does not account for the influence of vortices being shed from the blade tips into the wake on the induced velocity field. These tip vortices cre‐ ate multiple helical structures in the wake which play a major role in the induced velocity distribution at the rotor. To compensate for this deficiency in BEM theory, a tip-loss model originally developed by Prandtl is implemented as a correction factor to the induced veloci‐ ty field [8]. In the same way, a hub-loss model serves to correct the induced velocity result‐ ing from a vortex being shed near the hub of the rotor (see [25], [5].) Another modification needed in the BEM theory is the one developed by Glauert [7] to correct the rotor thrust co‐ efficient in the "turbulent-wake" state. This correction plays a key role when the turbine op‐ erates at high tip speed ratios and the induction factor is greater than about 0.45.

BEM theory was originally conceived for axisymmetric flow. Often, however, wind turbines operate at yaw angles relative to the incoming wind, which produces a skewed wake behind the rotor. For this reason, the BEM model needs also to be corrected to account for this skewed wake effect [31,22]. The influence of the wind turbine tower on the blade aerody‐ namics must also be modeled. We implemented the models developed by Bak et al. [2] and Powles [34] which provide the influence of the tower on the local velocity field at all points around the tower. This model contemplate increases in wind speed around the sides of the tower and the cross-stream velocity component in the tower near flow field.

Our model also incorporates the possibility to add multiple data tables for the different air‐ foils, and use them in real-time according to the instantaneous aerodynamic situations on the rotor. It also uses the Viterna's extrapolation method [36] to ensure the data availability for a range of angles of attack ±180 .
