Distributed blade aerodynamic properties.
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  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.
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 , 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
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
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
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
where denotes the position of the center of the tern along
When the structure is deformed due to loading, the original reference line
The next step is to find a strain energy expression asymptotically correct up to the second order of
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 are the classical Timoshenko strain measures due to extension, torsion and bending, and the transverse shear strains.
What we need to find is
where is the array of Timoshenko measures of deformation regrouped in a more convenient way, and .
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 .
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 ). 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 ; 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 , and weighted in an energy-consistent way according to Patil et al. , which produces the following variational formulation:
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
This variational formulation was discretized by the spectral-element method (see [20,29]). The magnitudes in (9) where replaced by their interpolated counterparts:
are the interpolation-function arrays, and
is a vector containing the nodal values of both the generalized velocities and the generalized forces. Superscript
corresponds to the discretization of terms 1 and 2 giving the equivalent of a mass matrix.
, corresponding to terms 3 and 8, is the stiffness matrix of the 1-D problem.
, corresponding to terms 4 and 9, is the additional stiffness related with the twist and curvature of the undeformed configuration.
corresponds to the evaluation of term 6, the contribution of the distributed loads; and
is an array containing the nodal values of the generalized distributed loads.
A linearization of around any given configuration gives the matrix
Matrix is a rearrangement of the components of the generalized-forces vector defined above. Matrix is a rearrangement of the components of the generalized-momentum vector , 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 .
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 the discretized expression
and updating the global vector of nodal values of the generalized velocities and forces as .
From the steady-state solution we also obtained the vibrational modes of the blade structure 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 )
2.2. Aerodynamic model
The flow model that interacts with the structural counterpart presented in section 2.1, called
The tendency in the wind-turbine industry to increase the size of the state-of-the-art machine  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
In what follows we will describe the main characteristics of our model, and refer to  and  for details on the classical
For instance, we could write the wind velocity vector
is the incoming wind velocity projected into the
Then, to compute the relative velocity affecting a blade element, we will project
Thus, the conning rotation matrix is a linear operator with a basic rotation taking into account the conning angle for the rotor, and the Pitching rotation matrix , represents a rotation around the pitch of the blade.
where 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
After applying the
, one more projection is needed to get to the instantaneous coordinate system associated with
After all these projections of the
where the addition of corresponds to the blade section structural deformation velocities, coming from the structural model.
Then, the magnitude
and the angle of attack
The aerodynamic loads acting on the blade element is then projected back onto the
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
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 . 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 , .) Another modification needed in the BEM theory is the one developed by Glauert  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.  and Powles  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  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 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
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 ). 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.
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 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 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 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 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
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 . 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 .
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 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
Figure 6 shows the displacement of the blade's reference-line (blade axis)
when it is subjected to the aerodynamic steady load in normal operational conditions. Figure 7, shows the corresponding rotations of the blade sections
From figure 6 we can see that the displacement
In figure 7, angles
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.
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  using FAST  and ADAMS  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.
Figures 8 and 9 show the amplitude of the deformation along the span for the three components of
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
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
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
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.
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.
project [internet] http://www.upwind.eu/ May 2012
Bak C. Madsen H. Aagaard Johansen J. Influence from blade-tower interaction on fatigue loads and dynamic (poster) 2 6 2001
Bathe K. J. Finite element procedures Prentice Hall, Englewood Cliffs New Jersey, USA 1996
Berdichevsky V. L. Variational--asymptotic method of constructing a theory of shells 43 664 687 1979
Burton T. Sharpe D. Jenkins N. Bossanyi E. Wind Energy Handbook Wiley Chichester, UK 2001
Cesnik C. E. S. Sutyrin V. G. Hodges D. H. Refined theory of composite beams: The role of short-wavelength extrapolation 33 1387 1407 1996
Glauert H. A general theory of the autogyro 1111 1926
Glauert H. Airplane propellers 4 169 360 1935
Griffin D. A. Blade system design studies volume I: Composite technologies for large wind turbine blades. Report SAND2002-1879, Sandia National Laboratories 2002
Griffin D. A. Evaluation of design concepts for adaptive wind turbine blades. Report SAND2002-2424, Sandia National Laboratories 2002
Hansen M. O. L. Sorensen J. N. Vousitas S. Sorensen N. Madsen H. Aa. State of the art in wind turbine aerodynamics and aeroelasticity 42 285 330 2006
Hodges D. H. Geometrically exact, intrinsic theory for dynamics of curved and twisted anisotropic beams 41 1131 1137 2003
Hodges D. H. Nonlinear Composite Beam Theory AIAA, Reston Virginia 2006
subtask 2: Research for deeper waters [internet] http://www.ieawind.org/AnnexXXIII/Subtask2.html December 2011
IEC. Wind turbine generator systems -- part 13: Measurement of mechanical loads. Report IEC/TS 61400--13, International Electrotechnical Commission (IEC) 2001
Laino D. J. Jonkman J. NWTC design codes ADAMS2AD [internet] http://wind.nrel.gov/designcodes/simulators/adams2ad/ May 2012
Jonkman J. Butterfield S. Musial W. Scott G. Definition of a 5-MW reference wind turbine for offshore system development. Technical Report NREL/TP-500-38060, National Renewable Energy Laboratory 2009
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 4--6 December 2007 Berlin, Germany 2007
Jonkman J. M. Buhl M. L. Jr. Fast user's guide. Technical Report NREL/EL-500-38230, National Renewable Energy Laboratory (NREL) Golden, Colorado, USA 2005
Karniadakis G. E. Bullister E. T. Patera A. T. A spectral element method for solution of two- and three-dimensional time-dependent incompressible navier-stokes equations 803 New York/Berlin 1985 Springer-Verlag
Kooijman H. J. T. Lindenburg C. Winkelaar D. Hooft E. L. van der Dowec 6 MW pre-design. Technical report, ECN-CX-01-135, Energy Research Center of the Netherlands, Petten 2003
Leishman J. G. Principles of helicopter aerodynamics Cambridge University Press Cambridge, UK 2006
Lindenburg C. Aeroelastic modelling of the LMH64-5 blade 2002
Locke J. Hidalgo I. Contreras The implementation of braided composite materials in the design of a bend-twist coupled blade. Report SAND2002-2425, Sandia National Laboratories 2002
Manwell J. F. McGowan J. G Rogers A. L. Wind energy explained: Theory, design and application Wiley Chichester, UK 2002
Wind power today. Report DOE/GO-102005-2115, U.S. Department of Energy 2005
Otero A. D. Ponta F. L. Structural analysis of wind-turbine blades by a generalized Timoshenko beam model 132 011015 2010
Passon P. Kühn M. Butterfield S. Jonkman J. Camp T. Larsen T.J. OC3--Benchmark exercise of aero-elastic offshore wind turbine codes 75 012071 IOP Publishing 2007
Patera A. T. A spectral element method for fluid dynamics: laminar flow in a channel expansion 54 468 488 1984
Patil M. J. Althoff M. Energy--consistent, Galerkin approach for the nonlinear dynamics of beams using mixed, intrinsic equations 1 9 AIAA 2006
Pitt D. M. Peters D. A. Theoretical prediction of dynamic-inflow derivatives 5 1 21 34 1981
Ponta F. L. The kinematic Laplacian equation method 207 405 426 2005
Ponta F. L. The KLE method: a velocity-vorticity formulation for the Navier-Stokes equations 73 1031 1038 2006
Powles S. R. J. The effects of tower shadow on the dynamics of a horizontal-axis wind turbine 7 26 42 1983
Parametric study for large wind turbine blades. Report SAND2002-2519, Sandia National Laboratories 2002
Viterna L. A. Janetzke D. C. Theoretical and experimental power from large horizontal-axis wind turbines. Technical report, National Aeronautics and Space Administration, Cleveland, OH (USA) Lewis Research Center 1982
Yamaguchi F. Curves and surfaces in computer aided geometric design Springer-Verlag Berlin 1988
Yu W. Hodges D. H. Generalized Timoshenko theory of the variational asymptotic beam sectional analysis 50 46 55 2005
Yu W. Hodges D. H. Volovoi V. Cesnik C. E. S. On Timoshenko-like modeling of initially curved and twisted composite beams 39 5101 5121 2002