Structural Analysis of Complex Wind Turbine Blades: Flexo- Torsional Vibrational Modes

1. Introduction

Limitations in the current blade technology constitute a technological barrier that needs to be broken in order to continue the improvement in wind-energy cost. Blade manufacturing is mostly based on composite laminates, which is labor-intensive and requires highly-qualified manpower. It constitutes a bottleneck to turbine upscaling that reflects into the increasing share of the cost of the rotor, within the total cost of the turbine, as turbine size increases. Figure 1 shows a compilation of data by NREL-DOE [26] on the proportional cost of each subsystem for different sizes of wind-turbines, where the systematic increase of the rotor cost share is clearly reflected.

Figure 1.

Evolution of the proportional cost for the different wind-turbine subsystems, as size increases (data compilation from [26]).

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 extrapolation of experimental data from the wind-tunnel to the prototype scale. Hence, computer 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 characterize the coupled aeroelastic problem still exceed the capacities of existing commercial simulation 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.

Hence, a key factor for a breakthrough in wind-turbine technology is to reduce the uncertainties related to blade dynamics, by the improvement of the quality of numerical simulations 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 fundamental 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 full-scale conditions with a reasonable computational cost. To this end, we developed a generalized Timoshenko code [27] based on a modified implementation the Variational-Asymptotic Beam Sectional technique (VABS) proposed by Hodges et al. (see [13] and references therein). 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 equations, called the KLE model [32,33], which would offer performance advantages over the present fluid-structure solvers.

In this chapter we present a set of tools for the design and full-scale analysis of the dynamics 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 the blades.

2. Fluid--Structure Interaction Model

2.1. Structural model

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 assumptions made in classical theories (like the Bernoulli or the standard Timoshenko approaches) 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 fluid-structure interaction model capable of dealing with the complex features of new-generation 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 twisted composite beams that uses the same variables as classical Timoshenko beam theory, but the hypothesis of beam sections remaining planar after deformation is abandoned. Instead, 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 fundamental 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 exterior 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 allows 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 interaction analysis. This translates into substantial savings in computational cost as the structural 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 warping 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 using the already-calculated 3-D warping functions.

Figure 3.

Generalized Timoshenko theory: Schematic of the reference line, orthogonal triads, and beam sections before 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 defined in such a way that is tangent to R and are contained into the section plane; with a correspondent coordinate system (X¹,X²,X³) where X¹ is the coordinate along R and X²,X³ are the Cartesian coordinates on the section plane. The position of a generic point on each section may be written as

where denotes the position of the center of the tern along R, and the index α assumes values 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

where w_i are the contribution to the displacement due to warping. Now, we are able to compute 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 ${\mathrm{\Gamma }}_{ij}=\frac{1}{2}\left(F_{ij}+F_{ji}\right)-{\delta }_{ij}$ , which provides a suitable measure of the 3-D strain field in terms of the beam strain measures and arbitrary warping functions. $\mathrm{\Gamma }$ is then used to compute the strain energy density function as

$$2U=⟨\mathbf{⟨}{\Gamma }^T|\mathbf{S}\Gamma ⟩⟩,$$ E3

were, S is the matrix of the characteristics of the material expressed in the coordinates, and $$⟨⟨·⟩⟩={\displaystyle {\int }_S·\sqrt{G}}d{X}^{2}d{X}^3$$ , 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/R₀ , where h is the characteristic size of the section, l the characteristic wavelength of deformation along the beam axis, and R₀ the characteristic radii of initial curvatures 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 measures (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:

$$2U={\epsilon }^{T}A\epsilon +{\epsilon }^{T}2B{\epsilon }^{\prime }+{\epsilon }^{\prime T}C{\epsilon }^{\prime }+{\epsilon }^{T}2D{\epsilon }^{\prime \prime },$$ E4

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 X¹ , and $$\epsilon ={\left[{\overline{\gamma }}_{11}{\overline{\kappa }}_1{\overline{\kappa }}_2{\overline{\kappa }}_3\right]}^T$$ , are the strain measures defined in the classical Bernoulli beam theory: ${\overline{\gamma }}_{11}$ the extension of the beam reference line, ${\overline{k}}_1$ its torsion, and ${\overline{k}}_2$ and ${\overline{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 discretization by finite-element techniques of the warping functions w_i defined in expression (2).

During this procedure, a set of four constraints must be applied on w_i . These restrictions, defined as $⟨w_i⟩=0$ and $⟨X^2{w}_3-X^3{w}_2⟩=0$ , where $⟨\bullet ⟩={\int }_s\bullet d{X}^{2}d{X}^3$ , 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 ε. Previous 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 constrained 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 requires 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 generalized Timoshenko model of the form

$$2U=\left[\begin{array}{cc}{\epsilon }^T& {\gamma }_s^T\end{array}\right]\left[\begin{array}{cc}\mathbf{X}& \mathbf{Y}\\ {\mathbf{Y}}^T& \mathbf{G}\end{array}\right]\left[\begin{array}{c}\epsilon \\ {\gamma }_s\end{array}\right]={\epsilon }^T\mathbf{X}\epsilon +2{\epsilon }^T\mathbf{Y}{\gamma }_s+{\gamma }_s^T\mathbf{G}{\gamma }_s,$$ E5

where $ϵ={\left[\begin{array}{cccc}{\gamma }_{11}& {\kappa }_1& {\kappa }_2& {\kappa }_3\end{array}\right]}^T$ are the classical Timoshenko strain measures due to extension, torsion and bending, and $\gamma ={\left[\begin{array}{cc}2{\gamma }_{12}& 2{\gamma }_{13}\end{array}\right]}^T$ 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

$$\epsilon =\epsilon +{\mathbf{Q}}_{\gamma }{{\gamma }^{\prime }}_s+{\mathbf{P}}_{\gamma }{\gamma }_s,$$ E6

where

$$\begin{array}{ccccc}{\mathbf{Q}}_{\gamma }^T=\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& -1& 0\end{array}\right],& & \text{y}& & {\mathbf{P}}_{\gamma }^T=\left[\begin{array}{cccc}0& K_2& -K_1& 0\\ 0& K_3& 0& -K_1\end{array}\right]\end{array},$$ E7

being K₁ the twist, and K₂ and K₃ the curvatures of the original reference line R. Thus, using (6), we may rewrite expression (4) in terms of the generalized Timoshenko strain measures 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 approximations 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 reordering of the matrix $\left[\begin{array}{cc}\mathbf{X}& \mathbf{Y}\\ {\mathbf{Y}}^T& \mathbf{G}\end{array}\right]$ }, in such a way as to get a functional for the strain energy density of expression~(5)

E8

where $$\overline{\gamma }=\left[\begin{array}{l}\gamma \\ \kappa \end{array}\right]$$ is the array of Timoshenko measures of deformation regrouped in a more convenient way, ${\gamma }^T=\left[\begin{array}{ccc}{\gamma }_{11}& 2{\gamma }_{12}& 2{\gamma }_{13}\end{array}\right]$ and ${\kappa }^T=\left[\begin{array}{ccc}{\kappa }_1& {\kappa }_2& {\kappa }_3\end{array}\right]$ .

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 geometric coordinates, provide the additional ability of reducing the so-called skin-error on curvilinear 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 tri-quadrilateral 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 triangle, we can significantly reduce the number of nodes to solve in the final system, subsequently recovering the values for the internal nodes from the solution on the non-condensable 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 approximately 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-element 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-quadrilateral technique may be found in Ponta [33]; including a schematic example of a mesh of tri-quadrilateral 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 following variational formulation:

where

                                                        $$\begin{array}{l}\mathbf{F}=\left[\begin{array}{l}\mathbf{F}\\ \mathbf{M}\end{array}\right],\overline{\mathbf{V}}=\left[\begin{array}{l}\mathbf{V}\\ \mathbf{\Omega }\end{array}\right],\overline{\mathbf{f}}=\left[\begin{array}{l}\mathbf{f}\\ \mathbf{m}\end{array}\right],\\ \widehat{\mathbf{\gamma }}=\left[\begin{array}{l}\tilde{\mathbf{\kappa }}\mathbf{0}\\ \widehat{\gamma }\tilde{\mathbf{\kappa }}\end{array}\right],\widehat{\mathbf{V}}=\left[\begin{array}{l}\tilde{\mathbf{\Omega }}0\\ \tilde{\mathbf{V}}\tilde{\mathbf{\Omega }}\end{array}\right],\widehat{\mathbf{K}}=\left[\begin{array}{l}\tilde{\mathbf{K}}0\\ {\tilde{\mathbf{e}}}_1\tilde{\mathbf{K}}\end{array}\right].\end{array}$$

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 $\tilde{\mathbf{A}}\mathbf{B}$ is equivalent to the cross product A B. Thus, $\tilde{\gamma }$ is associated with $\tilde{\kappa }$ with κ , $\tilde{\mathbf{V}}$ with V, and so forth. Hence, matrix $\tilde{\gamma }$ is a rearrangement of the components of the strain-measures vector $\overline{\gamma }$ defined above, the generalized-velocities vector $\overline{\mathbf{V}}$ and matrix $\tilde{\mathbf{V}}$ represent the components of the linear and angular velocities, and matrix $\tilde{\mathbf{K}}$ represents the initial torsion and curvatures of the beam (matrix $\widetilde{{\mathbf{e}}_1}$ is the skew-symmetric matrix associated to ${\mathbf{e}}_1^T=\left[\begin{array}{ccc}1& 0& 0\end{array}\right]$ , the unit vector along X¹ ). The generalized-forces vector $\overline{\mathbf{F}}$ represents the forces and moments related with the strain measures , and the generalized-distributed-loads vector $\overline{\mathbf{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 $\overline{\mathbf{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 X¹ .

This variational formulation was discretized by the spectral-element method (see [20,29]). The magnitudes in (9) where replaced by their interpolated counterparts: $\overline{\mathbf{V}}={\mathbf{H}}_{\overline{\mathbf{V}}}^e{\mathbf{Q}}^e$ , and $\overline{\mathbf{F}}={\mathbf{H}}_{\overline{\mathbf{F}}}^e{\mathbf{Q}}^e$ , where ${\mathbf{H}}_{\overline{\mathbf{V}}}^e$ and ${\mathbf{H}}_{\overline{\mathbf{F}}}^e$ are the interpolation-function arrays, and ${\mathbf{Q}}^e$ is a vector containing the nodal values of both the generalized velocities and the generalized forces. Superscript 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 derivatives of the magnitudes were interpolated in a similar way: ${\overline{\mathbf{V}}}^{\prime }={\mathbf{B}}_{\overline{\mathbf{V}}}^e{\mathbf{Q}}^e$ , and ${\overline{\mathbf{F}}}^{\prime }={\mathbf{B}}_{\overline{\mathbf{F}}}^e{\mathbf{Q}}^e$ , where ${\mathbf{B}}_{\overline{\mathbf{V}}}^e$ and ${\mathbf{B}}_{\overline{\mathbf{F}}}^e$ are the arrays for the interpolation-function derivatives. Then, we arrived to the discretized version of (9):

$$\delta Q^e{}^T{M}_1^e{\dot{Q}}^e=\delta Q^e{}^T\left(K_1^e+K_2^e\right)Q^e+\delta Q^e{}^T{K}_q^e{q}^{-e}+\delta Q^e{}^T{B}_Q^e\left(Q^e\right),$$ E10

where

                                      

${\mathbf{M}}_{\mathbf{1}}^e$ corresponds to the discretization of terms 1 and 2 giving the equivalent of a mass matrix. ${\mathbf{K}}_{\mathbf{1}}^e$ , corresponding to terms 3 and 8, is the stiffness matrix of the 1-D problem. ${\mathbf{K}}_{\mathbf{2}}^e$ , corresponding to terms 4 and 9, is the additional stiffness related with the twist and curvature of the undeformed configuration. ${\mathbf{K}}_{\mathbf{q}}^e$ corresponds to the evaluation of term 6, the contribution of the distributed loads; and ${\overline{\mathbf{q}}}^e$ 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 mapping from the problem coordinate X¹ 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

            

A linearization of ${\mathbf{B}}_Q^e\left({\mathbf{Q}}^e\right)$ around any given configuration ${\mathbf{Q}}_1^e$ gives the matrix

where

                                                  

Matrix $\overline{\mathbf{F}}$ is a rearrangement of the components of the generalized-forces vector $\overline{\mathbf{F}}$ defined above. Matrix $\overline{\mathbf{P}}$ is a rearrangement of the components of the generalized-momentum vector $\mathbf{P}=\left[\begin{array}{l}{\mathbf{P}}_v\\ {\mathbf{P}}_{\omega }\end{array}\right]$ , which represents the linear and angular momenta related with the generalized-velocities . Tilde operates in the same way defined before, and the subscript 1 indicates the value of the magnitudes at a given state ${\mathbf{Q}}_1^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 $\mathrm{\Delta }\mathbf{Q}$ the discretized expression

E11

and updating the global vector of nodal values of the generalized velocities and forces as ${\mathbf{Q}}^{(i+1)}={\mathbf{Q}}^{(i)}+\mathbf{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}}_{\mathbf{1}}\dot{\mathbf{Q}}+\left[{\mathbf{K}}_{\mathbf{1}}+{\mathbf{K}}_{\mathbf{2}}+{\mathbf{K}}_{\mathbf{N}}\left({\mathbf{Q}}^{(i)}\right)\right]\mathbf{Q}=\mathbf{0}.$$ E12

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])

E13

E14

where u is the vector of displacements of each point along the reference line from its position in the reference configuration to the one in the deformed configuration, and ${\mathbf{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 X¹ ). The strains γ and κ were computed from the generalized forces and the stiffness of the corresponding blade section. Equations (13) and (14) were also linearized, and like the other expressions, discretized by the spectral-element method.

2.2. Aerodynamic model

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 deformation 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 level 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 machine [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 perfectly valid, the actual mathematical formulation implies the assumption of blade sections remaining 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 analyzing BEM theory for this cases, the principle of equating the forces obtained by Blade Element 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-Mechanical train that define the rotor's angular speed, Ω .

For instance, we could write the wind velocity vector W_h 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 International Electrotechnical Commission (IEC) [15].

$${\mathbf{W}}_h=\left[\begin{array}{l}{W}_{w{h}_{(1)}}(1-a)\\ W_{w{h}_{(2)}}+\mathrm{\Omega }{r}_h(1+a^{\prime })\\ W_{w{h}_{(3)}}\end{array}\right],$$ E15

where ${\mathbf{W}}_{wh}$ is the incoming wind velocity projected into the h coordinate system, Ω is the angular velocity of the rotor and r_h is the radial distance of the airfoil section in the h coordinate system.

Then, to compute the relative velocity affecting a blade element, we will project W_h going through the different coordinate systems, from the hub, until reaching the blade's section coordinate 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 ${\mathbf{C}}_{{\theta }_{cn}}$ is a linear operator with a basic rotation taking into account the conning angle for the rotor, and the Pitching rotation matrix ${\mathbf{C}}_{{\theta }_p}$ , represents a rotation around the pitch of the blade.

$${\mathbf{C}}_{{\theta }_p}=\left[\begin{array}{ccc}\cos ({\theta }_p)& -\sin ({\theta }_p)& 0\\ \sin ({\theta }_p)& \cos ({\theta }_p)& 0\\ 0& 0& 1\end{array}\right],$$ E16

where ${\theta }_p$ is the pitch angle.

Two more re-orientations are needed in order to get to the instantaneous coordinate system 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 rectilinear). This curvatures can reflect either an initial twist along the longitudinal axis or a combination 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 position 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 ${\mathbf{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.

After applying the ${\mathbf{C}}_{Rb}$ , one more projection is needed to get to the instantaneous coordinate system associated with r. This last transformation is given by the ${\mathbf{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 W_h 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, ${\mathbf{W}}_{rel}$ :

$${\mathbf{W}}_{rel}=\left({\mathbf{C}}_{rR}{\mathbf{C}}_{Rb}{\mathbf{C}}_{{\theta }_p}{\mathbf{C}}_{{\theta }_{cn}}{\mathbf{W}}_h\right)+{\mathbf{v}}_{str}$$ E17

where the addition of ${\mathbf{v}}_{str}$ corresponds to the blade section structural deformation velocities, coming from the structural model.

Then, the magnitude $|{\mathbf{W}}_{rel}|$ and the angle of attack α are used to compute the forces on the airfoil section through the aerodynamic coefficients C_l , C_d . Another innovation of our model is that the data tables from static wind-tunnel are corrected at each timestep to consider either rotational-augmentation or dynamic-stall effects, or both.

The aerodynamic loads acting on the blade element is then projected back onto the h coordinate system,

$${\mathbf{d}\mathbf{F}}_h={\mathbf{C}}_{{\theta }_{cn}}^T{\mathbf{C}}_{{\theta }_p}^T{\mathbf{C}}_{Rb}^T{\mathbf{C}}_{rR}^T{\mathbf{C}}_{Lthal}{\mathbf{d}\mathbf{F}}_r$$ E18

where ${\mathbf{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 classical BEM theory, ${\mathbf{d}\mathbf{F}}_h$ is equated to the rate of change of momentum in the annular streamtube 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 create 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 velocity field [8]. In the same way, a hub-loss model serves to correct the induced velocity resulting 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 coefficient in the "turbulent-wake" state. This correction plays a key role when the turbine operates 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 aerodynamics 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 airfoils, 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 $\pm {180}^{\circ }$ .

3. Numerical Experimentation

In this section, we report some recent results of the application of our model to the analysis of a set of rotor blades based on the 5-MW Reference Wind Turbine (RWT) proposed by NREL [17]. We will start describing the structural features of the blade, its general aerodynamic properties, the blade internal structure, and the finite element meshes associated to the structural computations.

3.1. NREL Reference Wind Turbine

Based on the REpower 5MW wind turbine, NREL RWT was conceived for both onshore and offshore installations and is well representative of typical utility-scale multi megawatt commercial wind turbines. Although full specific technical data is not available for the REpower 5MW rotor blades, a baseline from a prototype blade was originally released by LM Glassfiber in 2001 for the Dutch Offshore Wind Energy Converter (DOWEC) 6MW wind turbine project~[21,23] and later re-adapted by NREL. In addition, the NREL 5-MW RWT project has been adopted as a reference model by the integrated European Union UpWind research program~[1] and the International Energy Agency (IEA) Wind Annex XXIII Subtask 2 Offshore Code Comparison Collaboration (OC3)~[14,18,28].

As stated in the NREL''s RWT project, the length of our rotor blade is set to be 61.5m. All basic aerodynamic properties as blade section chords, twist angles and basic spanwise stations distribution, correspond to the original data (see [17]). These aerodynamic properties, as well as the denomination of the basic airfoils at the design stations are included in table 1. Complementing the information in this table, figure 4 shows the blade section chords distribution along the span.

Station	Location [m]	Twist angle [°]	Chord length [m]	Airfoil type
1	0	13.3080	3.5420	Cylinder
2	1.3653	13.3080	3.5420	Cylinder
3	4.1020	13.3080	3.8540	Ellipsoid-1
4	6.8327	13.3080	4.1670	Ellipsoid-2
5	10.2520	13.3080	4.5570	DU 00-W-401
6	14.3480	11.4800	4.6520	DU 00-W-350
7	18.4500	10.1620	4.4580	DU 00-W-350
8	22.5521	9.0110	4.2490	DU 97-W-300
9	26.6480	7.7950	4.0070	DU 91-W-250
10	30.7500	6.5440	3.7480	DU 91-W-250
11	34.8520	5.3610	3.5020	DU 93-W-210
12	38.9479	4.1880	3.2560	DU 93-W-210
13	43.0500	3.1250	3.0100	NACA 64-618
14	47.1521	2.3190	2.7640	NACA 64-618
15	51.2480	1.5260	2.5180	NACA 64-618
16	54.6673	0.8630	2.3130	NACA 64-618
17	57.3980	0.3700	2.0860	NACA 64-618
18	60.1347	0.1060	1.4190	NACA 64-618
19	60.5898	0.0903	1.1395	NACA 64-618
20	61.0449	0.0783	0.7787	NACA 64-618

Table 1.

Distributed blade aerodynamic properties.

Figure 4.

Chord distribution along the blade.

The blade structure is a combination of two external aerodynamic shells, mounted on a box-beam spar which provides the main structural component to the aerodynamic forces. Analysing a blade section (see figure 2) we can see the aerodynamic shells plus two spar caps which, together with two shear webs, form the box-beam spar. Constructive characteristics as thickness as well as number and orientation of fiberglass layers for the different structural components of the blade sections are covered in detail in reports published by SANDIA National Labs. [35,9]. According to these reports, the aerodynamic shells are mainly composed by $\pm {45}^{\circ }$ layers, plus a small amount of randomly oriented fibers, gelcoat and filling resin. Shear webs, the the box-beam lateral walls, are made up of $\pm {45}^{\circ }$ layers with a balsa wood core which provides the needed buckling resistance. Shear webs are usually located at the 15% and 45% of the airfoil's chord but, for sections closer to the blade's root, the positions are modified in order to increase the section's stiffness. Focusing now on the spar caps, these are made of $0^{\circ }$ layers and are the most important structural element as they give support to the bending loads on the blade. Finally, the blade sections has a reinforcement at its rear part, i.e. the trailing edge spline, also made up of $0^{\circ }$ oriented fibers which supports the bending loads in the chord-wise direction. Reports [35,9] also provide a comprehensive description of lamination sequences and material properties.

Material properties within the subregions corresponding to each of the blade section components were assumed homogeneous and equal to those of an equivalent material. The properties of this equivalent material, a 6 × 6 symmetric matrix with 21 independent coefficients, were computed by a weighted average of the actual laminates properties. Since the thicknesses of the region layers are very small compared to the actual size of the blade section, this assumption does not introduce significant errors. Besides, if more detail is required, our computational codes allow for independent meshing of every single layer of material separately using the exact properties.

Figure 5.

Finite element meshes for morphed sections.]{Finite element meshes for morphed sections corresponding to: (a) 20% of the blade span, and (b) 60% of the blade span.

After the internal regions and materials were defined, a triquadrilateral mesh was generated for a number of blade sections along the span. The preset master sections in table 1 served as the basis for an innovative 3D-morphing technique based on variational cubic-spline interpolation which allows us to obtain smooth transitions between the known 2D airfoil sections along the span of the blade. This way one can divide the blade into any number of sections larger than the known ones and generate finite element meshes for a more refined study of the structural features. As an example, figure 5 shows the finite element meshes of two morphed airfoil sections located at the 20% and 60% of the blade span.

Using the technique described for the internal blade structure components, we refined 46 blade sections along the span to match the structural properties of the ones reported by NREL [17]. The main targeted properties to refine were edgewise, flapwise and torsional stiffness as well as mass density for every blade section. The pitch axis centering and the location of the aerodynamic coefficients reference points were also computed according to information in reference [17].

The general specifications of the turbine also match the ones on NREL's report. Thus, the rotor has an upwind orientation and is composed of three blades. The hub diameter is 3m and is located at 90m from the ground level. Total rotor diameter is 126m. It has a precone of ${2.5}^{\circ }$ and an overhang distance of 5m from the tower axis. The rated wind speed for this turbine is 11.4m/s.

3.2. Aeroelastic Steady State

After computing stiffness and inertia matrices for a discrete number of cross-sections along the span of the blade as described in section 2.1, the calculation of the aeroelastic steady state of the NREL RWT blades working under nominal conditions was solved by fully-coupling the structural and aerodynamic models presented in sections 2.1 and 2.2. Tip speed ratio for the nominal operational condition is λ=7, so the tangential velocity at the tip of the blade is 80m/s. For this nominal working condition, the power output computed for our rotor is 5.455MW which, taking into account that as in any BEM approach the interference of the tower and the nacelle is computed only approximately, is in very good agreement with the reported power for the NREL-5MW reference turbine rated at 5.296MW according to [17].

Figure 6 shows the displacement of the blade's reference-line (blade axis) ${\mathbf{U}}_h$ when it is subjected to the aerodynamic steady load in normal operational conditions. Figure 7, shows the corresponding rotations of the blade sections θ _h . These geometrical magnitudes were referred to a coordinate system, h from hub, aligned with the rotor's plane, according to standards from the International Electrotechnical Commission (IEC) [15]. Hence, the first unit vector is normal to the rotor's plane (i.e. axial) pointing downwind, the second is in the rotor's tangential direction pointing to the blade's trailing edge, and the third unit vector is in the radial direction pointing to the blade tip.

From figure 6 we can see that the displacement U_h1 of the blade's tip, normal to the rotor's plane, is 5.73m. This is perfectly consistent with results shown in [17]. Added to this, the tangential displacement U_h2 is 0.78m in the negative direction, meaning that aerodynamic forces are bending the blade in the direction towards its rotation as the rotor is producing a positive driving torque.

Figure 6.

Linear displacements of the reference-line U_h when the beam is subjected to a steady load in normal operational conditions (referred to a coordinate system aligned with the rotor's plane).

In figure 7, angles θ _h2 and θ _h1 are directly associated with blade bending in the normal and tangential directions to the rotor plane, that correspond to displacements U_h1 and U_h2 respectively. It is important to note that angle θ _h2 represents the angular displacements which takes the blade's axis out of the rotor's plane.

Figure 7.

Rotations of the beam sections θ _h when the beam is subjected to a steady load in normal operational conditions (referred to a coordinate system aligned with the rotor's plane).

3.3. Natural frequencies & Linear Modes

Vibrational modes around the aeroelastic steady-state are obtained from the solution of an eigenvalue problem as described in section 2.1. The resulting eigenvalues are complex conjugate, their imaginary part represent frequencies while their non-zero real part correspond to aerodynamic damping effects coming from non-conservative force fields in the 1D functional.

Mode	frequency [Hz]	Dominant U	Dominant θ
1	0.7066	Uh1	θ h2
2	1.0188	Uh2	θ h1
3	1.8175	Uh1	θ h2
4	3.3403	Uh2	θ h3
5	3.9493	Uh1	θ h2
6	6.4682	Uh2	θ h3
7	6.6851	Uh1	θ h3
8	8.0129	Uh2	θ h3
9	8.2403	Uh1	θ h2
10	9.7819	Uh1	θ h2

Table 2.

List of frequencies and dominant components of U_h and θ _h for the first ten modes of vibration.

Vibrational mode analysis provides relevant information about both the natural vibrational frequencies of the blade around a steady-state condition, and for the modes of deformation along the blade span. Table 2 summarizes the first 10 modes obtained showing the frequencies together with the corresponding dominant component for the displacements and the rotations of the blade section.

Table 3 shows a comparison of the frequencies obtained for the first 3 modes with the values reported by NREL in [17] using FAST [19] and ADAMS [16] software. FAST and ADAMS are considered today state-of-the-art softwares for structural blades analysis. From this comparison we see that the frequencies computed with our model match previous studies with a difference of 1% for the first mode and a maximum difference of 5% for the second and third modes. This difference is not surprising as the level of detail and richness of information that our computational tools can register is not present in the previous software like FAST or ADAMS.

Mode	frequency [Hz]	FAST	ADAMS
1	0.7066	0.6993	0.7019
2	1.0188	1.0793	1.0740
3	1.8175	1.9223	1.8558

Table 3.

Frequencies comparison for the first three modes according to NREL report [17].

Figures 8 and 9 show the amplitude of the deformation along the span for the three components of U_h and θ _h , normalized by the dominant component, for some of the deformation modes. Every mode shown includes displacements and rotations of the blade sections normalized by the value of the dominant component.

Figure 8.

Amplitude of U _h and θ _h for three vibrational modes around the aeroelastic steady-state configuration (normalized by the dominant component). From top to bottom modes # 1, 2 and 3.

Figure 9.

Amplitude of U _h and θ _h for three vibrational modes around the aeroelastic steady-state configuration (normalized by the dominant component). From top to bottom modes # 4, 7 and 10.

3.4. Recovery of 3-D fields

After computing the global deformation from the 1-D beam analysis, we recovered the corresponding 3-D fields (displacements, strains and stresses) using the 3-D warping functions previously calculated with our model. The knowledge of the stress state of every layer is of utter importance in the analysis of wind turbine blades in order to improve life-time and reliability of the design. Our model can provide the full six tensorial components of the stress tensor. Besides it can provide the 3 components of the displacement and 6 components of the strain.

For the previously solved aeroelastic steady-state, we present in figure 11 and 11 the six components of the Jaumann-Biot-Cauchy stress tensor Z=SΓ for the section located at 40% of the blade span. This region is particularly interesting as it combines energy production and structurally supports significant stress accumulation compared to other regions along the blade span. The dominant stress component, ${\mathbf{Z}}_{11}$ , at the top of figure 10 is the one primarily associated with the out of rotor-plane bending loads. Note here how the lower spar-cap is subjected to tensile stress while the upper one is under compression stress.

Figure 10.

Components of the Jaumann-Biot-Cauchy stress tensor Z = S Γ for the section located at 40% of the blade span (referred to the undeformed coordinate system (X¹,X²,X³) in Pa). From top to bottom: Z₁₁ , Z₁₂ and Z₁₃ .

Fig. 11.

Components of the Jaumann-Biot-Cauchy stress tensor Z= SΓ for the section located at 40% of the blade span (referred to the undeformed coordinate system (X¹,X²,X³) in Pa). From top to bottom: Z₂₂ , Z₂₃ and Z₃₃ .

4. Conclusions

With the method presented in this work we are able to model the structural behavior of wind turbine blades with the simplicity and economy of a one-dimensional model but with a level of description equivalent to a much more costly three-dimensional approach. The one-dimensional model is used to compute a fast, but accurate, solution for the deformed state of the blade when subjected to a steady load in normal operational conditions, and an analysis of the vibrational modes around this steady configuration. This provides a valuable tool to use during the design process. In that sense, the capacity of the Generalized-Timoshenko theory to capture the bending-twisting coupled modes in its fully populated 6 × 6 stiffness matrix for the 1-D beam problem plays a fundamental role.

Due to the geometrical complexity and material inhomogeneousness in the section, all the deformation modes of the blade are combined modes, i.e. there are no pure-flexural or pure-torsional modes. Plots of the vibrational modes may serve to identify eventual unstable states in the dynamic behavior of innovative prototype blades. Figures 8 and 9 show that, for certain modes, in some portions of the span, bending due to lift force occurs simultaneously with twisting in the sense that increases the angle of attack, and then, the lift force. A complete dynamic analysis of the fluid-structure interaction process would be needed to determine if those particular modes would be activated or not during the blade operation. Nevertheless, having the possibility of quickly identifying those modes (and their associated frequencies) at an early stage of the design process seems very helpful.

Regarding the linear vibrational modes depicted in figure 8, the first mode shows mainly out-of-plane curvature corresponding to the fundamental frequency of the blade but as the blade is initially twisted and has complex inhomogeneous sections, it is not a pure bending mode. Therefore, as a result of the non conventional couplings, curvature deformation in the rotor plane and also torsion appears in this mode. The second mode is also commanded by out-of-plane curvature but it has an important torsion component, while mode number three is similar to the first one with a higher wave number.

Regarding the modes presented in figure 9, the fourth and seventh modes are mainly commanded by torsion and hence they could be responsible for fluttering of the blade in case they are excited by the interaction with the surrounding fluid. The tenth mode is also principally a flap-wise curvature mode but with higher wave number than the first and third. It shows a more complex behavior arising from complicated couplings among different deformations. The frequencies of the first three modes presented in table 3 are in good agreement with published results obtained by other models.

The above-mentioned flexo-torsional characteristic also gives this model the ability to simulate the dynamic performance of adaptive blades, at an affordable computational cost. In the Adaptive-Blade concept (see [10,24], among others), tailoring of the flexo-torsional modes of the blade is used to reduce aerodynamic loads by controlling the coupling between bending and twisting. As the blade bends under load, the angle of attack of the airfoil sections changes, reducing the lift force. Limiting extreme loads and improving fatigue performance, this added passive control reduces the intensity of the actuation of the active control system. Plots like figure 6 and figure 7 provide valuable information about the simultaneous deformation of twisting and bending under a given load.

Recovering of the stress tensor components for the different zones of the blade section helps in the prediction of stress concentration in the basic design that may ultimately lead to eventual material failure. More exhaustive fatigue studies can be conducted analyzing the stress both in the steady state or in time-marching solutions of the problem. The capability of computing the whole 6 components of the stress tensor makes it possible to apply sophisticated failure theories.

Our aeroelastic model may also be used to simulate the dynamic response of the wind turbine tower. In that case, the structural model would be applied to the tower to obtain the stiffness matrices of the equivalent beam as it is done with the blades. The aerodynamic loads would be computed from the aerodynamic coefficients of the cylindrical sections of the tower using basically the same subroutines. As in the case of the blades, all the complex flexo-torsional modes of deformation of the tower would be taken into account, and the associated vibrational effects included in the general analysis of the whole turbine.

We plan to continue our work with a dynamic simulation of the fluid-structure problem. In a first stage, we will couple the phenomena by feeding back changes in geometry due to blade deformation in our aerodynamic model and recomputing the forces. At this stage, we also plan to include statistically-generated perturbations to represent fluctuations in wind speed and direction based on anemometry data for wind resource in several representative locations. Besides providing us with a fast model for a quick analysis, this model will serve as an intermediate step before the ultimate goal of coupling the structural model with the velocity-vorticity KLE approach mentioned above.

Acknowledgments

The authors are very grateful for the financial support made available by the National Science Foundation through grants CEBET-0933058 and CEBET-0952218 and University of Buenos Aires through grant 20020100100536 UBACyT 2011/14.