Distributed blade aerodynamic properties.
1. Introduction
Limitations in the current blade technology constitute a technological barrier that needs to be broken in order to continue the improvement in windenergy cost. Blade manufacturing is mostly based on composite laminates, which is laborintensive and requires highlyqualified 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 NRELDOE [26] on the proportional cost of each subsystem for different sizes of windturbines, where the systematic increase of the rotor cost share is clearly reflected.
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 windtunnel to the prototype scale. Hence, computer models of fluidstructure interaction phenomena are particularly relevant to the design and optimization of windturbines. The windturbine 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 windturbine technology is to reduce the uncertainties related to blade dynamics, by the improvement of the quality of numerical simulations of the fluidstructure interaction process, and by a better understanding of the underlying physics. The current stateoftheart is to solve the aeroelastic equations in a fully nonlinear 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 fullscale conditions with a reasonable computational cost. To this end, we developed a
In this chapter we present a set of tools for the design and fullscale analysis of the dynamics of composite laminate windturbine 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. FluidStructure 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
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 1D finiteelement problem. Thus, this way we are able to decouple a general 3D nonlinear anisotropic elasticity problem into a linear, 2D, crosssectional analysis (that may be solved
In order to make this chapter selfcontained, 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
where denotes the position of the center of the tern along
When the structure is deformed due to loading, the original reference line
where
were,
The next step is to find a strain energy expression asymptotically correct up to the second order of
where
The variational asymptotic procedure to get the matrices in equation (4) involves the discretization by finiteelement techniques of the warping functions
During this procedure, a set of four constraints must be applied on
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
where
What we need to find is
where
being
where
For the discretization of the 2D sections, we adopted the triquadrilateral finiteelement technique, which is based on the use of ninenode 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 socalled skinerror on curvilinear boundaries when compared to linear elements. For a detailed description of the isoparametricelement technique and its corresponding interpolation functions see Bathe [3].
In order to combine the advantages of the ninenode quadrilateral isoparametric element with the geometrical ability of a triangular grid to create suitable nonstructured meshes with gradual and smooth changes of mesh density, we implemented what we called triquadrilateral isoparametric elements. The triquadrilateral elements consist of an assembling of three quadrilateral ninenode 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 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 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 finiteelement and spectralelement methods but also because it improves the condition number of the final condensed system. This is related with the properties of the Schurcomplement technique. The condensed system is essentially the Schur complement of the interiornode submatrix in the noncondensed original system. A detailed description of the triquadrilateral 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 triquadrilateral element.
To solve the onedimensional 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 energyconsistent way according to Patil et al. [30], which produces the following variational formulation:
where
Tilde indicates the skewsymmetric matrix associated to a vector magnitude in such a way that, for example, if we have any pair of vectors
This variational formulation was discretized by the spectralelement method (see [20,29]). The magnitudes in (9) where replaced by their interpolated counterparts:
where
A linearization of
where
Matrix
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
and updating the global vector of nodal values of the generalized velocities and forces as
From the steadystate solution we also obtained the vibrational modes of the blade structure and their corresponding frequencies by solving the eigenvalue problem
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])
where
2.2. Aerodynamic model
The flow model that interacts with the structural counterpart presented in section 2.1, called
The tendency in the windturbine industry to increase the size of the stateoftheart 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 preconforming processes, produces high rotations of the blade sections. Moreover, blades could be preconformed with specific curvatures given to any of their axis
In what follows we will describe the main characteristics of our model, and refer to [25] and [5] for details on the classical
For instance, we could write the wind velocity vector
where
Then, to compute the relative velocity affecting a blade element, we will project
Thus, the conning rotation matrix
where
Two more reorientations 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 preconformed curvatures along its longitudinal axis
After applying the
After all these projections of the
where the addition of
Then, the magnitude
The aerodynamic loads acting on the blade element is then projected back onto the
where
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 tiploss model originally developed by Prandtl is implemented as a correction factor to the induced velocity field [8]. In the same way, a hubloss 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 "turbulentwake" 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 crossstream 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 realtime 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
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
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 utilityscale 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
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.







1  0  13.3080  3.5420  Cylinder 
2  1.3653  13.3080  3.5420  Cylinder 
3  4.1020  13.3080  3.8540  Ellipsoid1 
4  6.8327  13.3080  4.1670  Ellipsoid2 
5  10.2520  13.3080  4.5570  DU 00W401 
6  14.3480  11.4800  4.6520  DU 00W350 
7  18.4500  10.1620  4.4580  DU 00W350 
8  22.5521  9.0110  4.2490  DU 97W300 
9  26.6480  7.7950  4.0070  DU 91W250 
10  30.7500  6.5440  3.7480  DU 91W250 
11  34.8520  5.3610  3.5020  DU 93W210 
12  38.9479  4.1880  3.2560  DU 93W210 
13  43.0500  3.1250  3.0100  NACA 64618 
14  47.1521  2.3190  2.7640  NACA 64618 
15  51.2480  1.5260  2.5180  NACA 64618 
16  54.6673  0.8630  2.3130  NACA 64618 
17  57.3980  0.3700  2.0860  NACA 64618 
18  60.1347  0.1060  1.4190  NACA 64618 
19  60.5898  0.0903  1.1395  NACA 64618 
20  61.0449  0.0783  0.7787  NACA 64618 

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. 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 boxbeam 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
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
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 3Dmorphing technique based on variational cubicspline 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
3.2. Aeroelastic Steady State
After computing stiffness and inertia matrices for a discrete number of crosssections 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
Figure 6 shows the displacement of the blade's referenceline (blade axis)
From figure 6 we can see that the displacement
In figure 7, angles
3.3. Natural frequencies & Linear Modes
Vibrational modes around the aeroelastic steadystate 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 nonzero real part correspond to aerodynamic damping effects coming from nonconservative force fields in the 1D functional.





1  0.7066 



2  1.0188 



3  1.8175 



4  3.3403 



5  3.9493 



6  6.4682 



7  6.6851 



8  8.0129 



9  8.2403 



10  9.7819 




Vibrational mode analysis provides relevant information about both the natural vibrational frequencies of the blade around a steadystate 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 stateoftheart 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  

Figures 8 and 9 show the amplitude of the deformation along the span for the three components of
3.4. Recovery of 3D fields
After computing the global deformation from the 1D beam analysis, we recovered the corresponding 3D fields (displacements, strains and stresses) using the 3D 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 lifetime 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 steadystate, we present in figure 11 and 11 the six components of the JaumannBiotCauchy stress tensor
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 onedimensional model but with a level of description equivalent to a much more costly threedimensional approach. The onedimensional 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 GeneralizedTimoshenko theory to capture the bendingtwisting coupled modes in its fully populated
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 pureflexural 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 simultaneously with twisting in the sense that increases the angle of attack, and then, the lift force. A complete dynamic analysis of the fluidstructure 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 outofplane 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 outofplane 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 flapwise 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 abovementioned flexotorsional characteristic also gives this model the ability to simulate the dynamic performance of
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 timemarching 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 flexotorsional 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 fluidstructure 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 statisticallygenerated 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 velocityvorticity KLE approach mentioned above.
Acknowledgments
The authors are very grateful for the financial support made available by the National Science Foundation through grants CEBET0933058 and CEBET0952218 and University of Buenos Aires through grant 20020100100536 UBACyT 2011/14.
References
 1.
project [internet] http://www.upwind.eu/ May 2012  2.
Bak C. Madsen H. Aagaard Johansen J. Influence from bladetower interaction on fatigue loads and dynamic (poster)  3.
Bathe K. J. Finite element procedures Prentice Hall, Englewood Cliffs New Jersey, USA 1996  4.
Berdichevsky V. L. Variationalasymptotic method of constructing a theory of shells  5.
Burton T. Sharpe D. Jenkins N. Bossanyi E. Wind Energy Handbook Wiley Chichester, UK 2001  6.
Cesnik C. E. S. Sutyrin V. G. Hodges D. H. Refined theory of composite beams: The role of shortwavelength extrapolation  7.
Glauert H. A general theory of the autogyro  8.
Glauert H. Airplane propellers  9.
Griffin D. A. Blade system design studies volume I: Composite technologies for large wind turbine blades. Report SAND20021879, Sandia National Laboratories 2002  10.
Griffin D. A. Evaluation of design concepts for adaptive wind turbine blades. Report SAND20022424, Sandia National Laboratories 2002  11.
Hansen M. O. L. Sorensen J. N. Vousitas S. Sorensen N. Madsen H. Aa. State of the art in wind turbine aerodynamics and aeroelasticity  12.
Hodges D. H. Geometrically exact, intrinsic theory for dynamics of curved and twisted anisotropic beams  13.
Hodges D. H. Nonlinear Composite Beam Theory AIAA, Reston Virginia 2006  14.
subtask 2: Research for deeper waters [internet] http://www.ieawind.org/AnnexXXIII/Subtask2.html December 2011  15.
IEC. Wind turbine generator systems  part 13: Measurement of mechanical loads. Report IEC/TS 6140013, International Electrotechnical Commission (IEC) 2001  16.
Laino D. J. Jonkman J. NWTC design codes ADAMS2AD [internet] http://wind.nrel.gov/designcodes/simulators/adams2ad/ May 2012  17.
Jonkman J. Butterfield S. Musial W. Scott G. Definition of a 5MW reference wind turbine for offshore system development. Technical Report NREL/TP50038060, National Renewable Energy Laboratory 2009  18.
Jonkman J. Butterfield S. Passon P. Larsen T. Camp T. Nichols J. Azcona J. Martinez A. Offshore code comparison collaboration within IEA Wind Annex XXIII: Phase II results regarding monopile foundation modeling  19.
Jonkman J. M. Buhl M. L. Jr. Fast user's guide. Technical Report NREL/EL50038230, National Renewable Energy Laboratory (NREL) Golden, Colorado, USA 2005  20.
Karniadakis G. E. Bullister E. T. Patera A. T. A spectral element method for solution of two and threedimensional timedependent incompressible navierstokes equations  21.
Kooijman H. J. T. Lindenburg C. Winkelaar D. Hooft E. L. van der Dowec 6 MW predesign. Technical report, ECNCX01135, Energy Research Center of the Netherlands, Petten 2003  22.
Leishman J. G. Principles of helicopter aerodynamics Cambridge University Press Cambridge, UK 2006  23.
Lindenburg C. Aeroelastic modelling of the LMH645 blade  24.
Locke J. Hidalgo I. Contreras The implementation of braided composite materials in the design of a bendtwist coupled blade. Report SAND20022425, Sandia National Laboratories 2002  25.
Manwell J. F. McGowan J. G Rogers A. L. Wind energy explained: Theory, design and application Wiley Chichester, UK 2002  26.
Wind power today. Report DOE/GO1020052115, U.S. Department of Energy 2005  27.
Otero A. D. Ponta F. L. Structural analysis of windturbine blades by a generalized Timoshenko beam model  28.
Passon P. Kühn M. Butterfield S. Jonkman J. Camp T. Larsen T.J. OC3Benchmark exercise of aeroelastic offshore wind turbine codes  29.
Patera A. T. A spectral element method for fluid dynamics: laminar flow in a channel expansion  30.
Patil M. J. Althoff M. Energyconsistent, Galerkin approach for the nonlinear dynamics of beams using mixed, intrinsic equations  31.
Pitt D. M. Peters D. A. Theoretical prediction of dynamicinflow derivatives  32.
Ponta F. L. The kinematic Laplacian equation method  33.
Ponta F. L. The KLE method: a velocityvorticity formulation for the NavierStokes equations  34.
Powles S. R. J. The effects of tower shadow on the dynamics of a horizontalaxis wind turbine  35.
Parametric study for large wind turbine blades. Report SAND20022519, Sandia National Laboratories 2002  36.
Viterna L. A. Janetzke D. C. Theoretical and experimental power from large horizontalaxis wind turbines. Technical report, National Aeronautics and Space Administration, Cleveland, OH (USA) Lewis Research Center 1982  37.
Yamaguchi F. Curves and surfaces in computer aided geometric design SpringerVerlag Berlin 1988  38.
Yu W. Hodges D. H. Generalized Timoshenko theory of the variational asymptotic beam sectional analysis  39.
Yu W. Hodges D. H. Volovoi V. Cesnik C. E. S. On Timoshenkolike modeling of initially curved and twisted composite beams