**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 com‐ mercial 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 Glass‐ fiber 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 pro‐ gram~[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 sta‐ tions 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 distri‐ bution along the span.


**Table 1.** Distributed blade aerodynamic properties.

equated to the change in axial momentum, while the tangential component, is equated to the

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.

tower and the cross-stream velocity component in the tower near flow field.

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

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

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

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 com‐ mercial 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 Glass‐ fiber in 2001 for the *Dutch Offshore Wind Energy Converter (DOWEC) 6MW* wind turbine

change of angular momentum.

136 Advances in Wind Power

for a range of angles of attack ±180 .

**3. Numerical Experimentation**

**3.1. NREL Reference Wind Turbine**

structural computations.

**Figure 4.** Chord distribution along the blade.

The blade structure is a combination of two external aerodynamic shells, mounted on a boxbeam spar which provides the main structural component to the aerodynamic forces. Analy‐ sing 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 Na‐ tional Labs. [35,9]. According to these reports, the aerodynamic shells are mainly composed by ±45 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 ±45 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 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 oriented fibers which supports the bending loads in the chord-wise direction. Reports [35,9] also provide a comprehensive de‐ scription of lamination sequences and material properties.

Material properties within the subregions corresponding to each of the blade section compo‐ nents were assumed homogeneous and equal to those of an equivalent material. The proper‐ ties 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 thick‐ nesses 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 sepa‐ rately 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 inter‐ polation 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 lo‐ cation of the aerodynamic coefficients reference points were also computed according to in‐ formation in reference [17].

The general specifications of the turbine also match the ones on NREL's report. Thus, the ro‐ tor 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 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**

The blade structure is a combination of two external aerodynamic shells, mounted on a boxbeam spar which provides the main structural component to the aerodynamic forces. Analy‐ sing 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 Na‐ tional Labs. [35,9]. According to these reports, the aerodynamic shells are mainly composed by ±45 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 ±45 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 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 oriented fibers which supports the bending loads in the chord-wise direction. Reports [35,9] also provide a comprehensive de‐

Material properties within the subregions corresponding to each of the blade section compo‐ nents were assumed homogeneous and equal to those of an equivalent material. The proper‐ ties 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 thick‐ nesses 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 sepa‐

**Figure 5.** Finite element meshes for morphed sections.]{Finite element meshes for morphed sections corresponding

scription of lamination sequences and material properties.

rately using the exact properties.

138 Advances in Wind Power

to: (a) 20% of the blade span, and (b) 60% of the blade span.

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 fullycoupling 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 accord‐ ing to [17].

Figure 6 shows the displacement of the blade's reference-line (blade axis) **U***<sup>h</sup>* when it is sub‐ jected to the aerodynamic steady load in normal operational conditions. Figure 7, shows the corresponding rotations of the blade sections *θ h*. These geometrical magnitudes were refer‐ red to a coordinate system, *h* from *hub*, aligned with the rotor's plane, according to stand‐ ards 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 ro‐ tor'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 *Uh1* 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 tan‐ gential displacement *Uh2* 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 posi‐ tive driving torque.

**Figure 6.** Linear displacements of the reference-line **Uh** when the beam is subjected to a steady load in normal opera‐ tional 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 *Uh1* and *Uh2* respec‐ tively. 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 θ *<sup>h</sup>* when the beam is subjected to a steady load in normal operational condi‐ tions (referred to a coordinate system aligned with the rotor's plane).

### **3.3. Natural frequencies & Linear Modes**

gential displacement *Uh2* 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 posi‐

**Figure 6.** Linear displacements of the reference-line **Uh** when the beam is subjected to a steady load in normal opera‐

In figure 7, angles **θ h2** and **θ h1** are directly associated with blade bending in the normal and

tively. It is important to note that angle **θ h2** represents the angular displacements which

**Figure 7.** Rotations of the beam sections θ *<sup>h</sup>* when the beam is subjected to a steady load in normal operational condi‐

and *Uh2*

respec‐

tangential directions to the rotor plane, that correspond to displacements *Uh1*

tional conditions (referred to a coordinate system aligned with the rotor's plane).

takes the blade's axis out of the rotor's plane.

tions (referred to a coordinate system aligned with the rotor's plane).

tive driving torque.

140 Advances in Wind Power

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 conju‐ gate, 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.


**Table 2.** List of frequencies and dominant components of **Uh** 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 frequen‐ cies 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 compari‐ son 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 computa‐ tional tools can register is not present in the previous software like FAST or ADAMS.


**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 compo‐ nents of **Uh** and *θ h*, normalized by the dominant component, for some of the deformation modes. Every mode shown includes displacements and rotations of the blade sections nor‐ malized 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 (nor‐ malized 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 (nor‐ malized by the dominant component). From top to bottom modes # 4, 7 and 10.

#### **3.4. Recovery of 3-D fields**

Figures 8 and 9 show the amplitude of the deformation along the span for the three compo‐ nents of **Uh** and *θ h*, normalized by the dominant component, for some of the deformation modes. Every mode shown includes displacements and rotations of the blade sections nor‐

**Figure 8.** Amplitude of **U**h and θ *h* for three vibrational modes around the aeroelastic steady-state configuration (nor‐

**Figure 9.** Amplitude of **U**h and θ *h* for three vibrational modes around the aeroelastic steady-state configuration (nor‐

malized by the dominant component). From top to bottom modes # 1, 2 and 3.

malized by the dominant component). From top to bottom modes # 4, 7 and 10.

malized by the value of the dominant component.

142 Advances in Wind Power

After computing the global deformation from the 1-D beam analysis, we recovered the cor‐ responding 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 re‐ liability 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, **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 *(X1,X2,X3)* in Pa). From top to bottom: *Z11*, *Z12* and *Z13*.

**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 *(X1,X2,X3)* in Pa). From top to bottom: *Z22*, *Z23* and *Z33*.

## **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-Timo‐ shenko 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 puretorsional 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 simultane‐ ously 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 de‐ termine 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 com‐ manded 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 princi‐ pally 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 defor‐ mations. 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 simu‐ late 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 defor‐ mation of twisting and bending under a given load.

**Fig. 11.** Components of the Jaumann-Biot-Cauchy stress tensor **Z=** *S***Γ** for the section located at 40% of the blade span

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-Timo‐ shenko theory to capture the bending-twisting coupled modes in its fully populated *6 × 6*

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 puretorsional 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 simultane‐ ously 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 de‐

(referred to the undeformed coordinate system *(X1,X2,X3)* in Pa). From top to bottom: *Z22*, *Z23* and *Z33*.

stiffness matrix for the 1-D beam problem plays a fundamental role.

**4. Conclusions**

144 Advances in Wind Power

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 even‐ tual 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 com‐ puting 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 tur‐ bine 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 as‐ sociated 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.
