The present work focuses on the unsteady aerodynamics and aeroelastic properties of a smallmedium sized windturbine blade operating under ideal conditions. A tapered/twisted blade representative of commercial blades used in an experiment setup at the National Renewable Energy Laboratory is considered. The aerodynamic loads are computed using Computational Fluid Dynamics (CFD) techniques. For this purpose, FLUENT
^{®}
, a commercial finitevolume code that solves the NavierStokes and the ReynoldsAveraged NavierStokes (RANS) equations, is used. Turbulence effects in the 2D simulations are modeled using the Wilcox
kw
model for validation of the CFD approach. For the 3D aerodynamic simulations, in a first approximation, and considering that the intent is to present a methodology and workflow philosophy more than highly accurate turbulent simulations, the unsteady laminar NavierStokes equations were used to determine the unsteady loads acting on the blades. Five different blade pitch angles were considered and their aerodynamic performance compared. The structural dynamics of the flexible windturbine blade undergoing significant elastic displacements has been described by a nonlinear flaplagtorsion slenderbeam differential model. The aerodynamic quasisteady forcing terms needed for the aeroelastic governing equations have been predicted through a striptheory based on a simple 2D model, and the pertinent aerodynamic coefficients and the distribution over the blade span of the induced velocity derived using CFD. The resulting unsteady hub loads are achieved by a first space integration of the aeroelastic equations by applying the Galerkin’s approach and by a time integration using a harmonic balance scheme. Comparison among two and three dimensional computations for the unsteady aerodynamic load, the flap, lag and torsional deflections, forces and moments are presented in the paper. Results, discussions and pertinent conclusions are outlined.
Nomenclature
A
amplitude of oscillation
a
cylinder diameter
a_{’}
lift curve slope
C_{p}
pressure coefficient
C_{d}
drag coefficient
C_{l}
lift coefficient
C_{x}
force coefficient in the x direction
C_{y}
force coefficient in the y direction
c
chord
dt
time step
F_{x}
X component of the resultant pressure force acting on the vehicle
F_{y}
Y component of the resultant pressure force acting on the vehicle
f, g
generic functions
h
height
N, N_{b}
Number of blades of rotor
r
Radial distance from axis of rotation
R
Radius of rotor
V_{tot}
Resultant total velocity,
V
_{∞}
Freestream velocity
x
Dimensionless local spanwise direction
S_{x} , S_{r}, S_{z}
Drag, radial, and vertical shear forces at the rotor hub, respectively
N_{F}, N_{L}
Root bending and torque moments, respectively
T ; H ; Y
Total forces in the nonrotating frame for a N bladed windmill, respectively
M_{x} ; M_{y} ; Q
Total moments in the nonrotating frame for a N bladed windmill, respectively
 Subscripts
i
induced
r
radial distance from axis of rotation
∞
freestream conditions
 Superscript
2D
twodimensional
3D
threedimensional
 Greek Symbols
α
angleofattack
λ
inflow ratio
θ_{tw}
blade twist
ω_{R}, Ω
angular velocity of rotor
K
trailingedge (TE) nondimensional angular deflection rate
ψ
azimuthal position of the blade
1. Introduction
Wind turbines operate in an aerodynamic environment which is difficult to analyze and predict. The periodic and aperiodic unsteadiness of the atmospheric boundary layer makes the design of windenergy systems very difficult since their design is typically optimized for steady state conditions. In particular, the performance of the turbine blades is strongly affected by the stochastic nature of the airloads; the blade structure needs to be built for the best efficiency in converting the wind energy into usable mechanical energy while maintaining structural integrity
[1

4]
. The need to account for the random nature of the airloads leads to heavier, stiffer blade designs that are not optimized for the prevailing operating conditions. Also, the heavier blades limit the rotational speed of the rotor as the inertial forces acting at the root of the blades are too high for safe operation. The lessthanoptimal design of the blades leads to higher initial and operating costs since centrifugal and gravity loads are primarily responsible for fatigue failure
[1

4]
. These design weaknesses are contributing factors to the competitive disadvantage of wind energy against other forms of renewable and nonrenewable energy sources.
Options to improve the efficiency of wind turbine blades are varied. On one side of the spectrum there is the option to design lighter, stronger blades using more advanced materials. However, the high costs associated with the advanced structure and materials make this option economically not viable On the opposite side of the spectrum is the implementation of an active control system that, rather than modifying the blade structure to make it stronger and stiffer, allows the blade to morph and to change its shape in order to alleviate the instantaneous loads. By using this control strategy, the blade structural loads can be reduced and the operation of the windturbine can be continued under adverse environmental conditions. An additional advantage of an active control system is the effective stabilization of the blade which protects the wind turbine system against aeroelastic instabilities
[5
,
6]
.
However, the proposed control system concept needs development and verification. In order to achieve these goals, an aeroelastic model that couples the blade structural response to the unsteady aerodynamic loads is being developed. In order to validate the control system as a proofof concept, the aeroelastic model needs to be able to capture the characteristic behavior of both the fluids and structural domains and yet be relatively lowcost for quick turnaround design iterations.
One of the bottlenecks in the development of the aeroelastic model is how to best capture the aerodynamic behavior of the blade since such a prediction is critical to the reliable description of the system. Often aerodynamic modeling has been accomplished by means of engineeringtype empirical models like the ONERA or the Beddoes Leishman models
[7
,
8]
. The airloads models are based on the conventional rotorcraft lifting line approach of treating each 3D rotor blade as a series of 2D airfoil sections and calculating the section airloads based on local section aerodynamic parameters. This is mainly due to the fact that the considered blades have large aspect ratio. In addition, the aerodynamic loads are based on the blade element momentum theory which assumes a 2D flow over each blade airfoil section. The induced velocity is obtained from the simple momentum theory, generalized dynamic wake, prescribed wake, or free wake. Static stall and compressibility effects are implicitly included in tabulated data as functions of angleofattack. Another assumption is that the quasisteady aerodynamics governs the reverse flow region. NavierStokes solvers have been usually restricted to twodimensional analysis and only a few attempts have been made to advance the stateoftheart in aeroelasticity from twodimensional to quasithreedimensional and fully threedimensional NavierStokes solutions
[9
,
10]
.
In this paper, the aerodynamic loads are computed using Computational Fluid Dynamics (CFD) techniques. A commercial finitevolume code, FLUENT
^{®}
, that solves the NavierStokes and ReynoldsAveraged NavierStokes (RANS) equations is used. Turbulence effects, considered only in the 2D simulations, are modeled using Wilcox’s
kw
model
[10]
. For the 3D simulations, the unsteady laminar NavierStokes equations were used to compute the unsteady periodic loads. Moreover, the structural dynamics of the flexible windturbine blade undergoing significant elastic displacements has been described by the nonlinear flaplag torsion slenderbeam differential model introduced by Hodges and Dowell
[11]
. The aerodynamic quasisteady forcing terms needed for the aeroelastic governing equations have been predicted through a striptheory based on a simple 2D model
[11

13]
whose pertinent aerodynamic coefficients could be derived from Computational Fluid Dynamics (CFD) techniques
[14]
. The resulting unsteady hub loads are finally achieved by a first space integration of the aeroelastic equations by applying the Galerkin approach – using as shape functions the bending and torsion free vibration eigenmodes of a uniform cantilever beam – and then in time using the harmonic balance approach
[11
,
15
,
16]
. Referring to the threebladed wind turbine used for the experimental test at the National Renewable Energy Laboratory (NREL)
[17
,
18]
represented in
Fig. 1
, the goals of this paper are twofold: i) to develop an aerodynamic and aeroelastic tool for wind turbine design based on the fullythreedimensional NavierStokes equations, and ii) to highlight the differences and deficiencies of a twodimensional versus a threedimensional wind turbine aerodynamic modeling.
Isometric view of the threebladed rotor S809 airfoil. The insert shows the orientation of the axes of reference.
Desired outputs from this numerical tool are the vibration levels due to the unsteady loads acting on the rotating windturbine blades, and the prediction of the fatiguelife of the structure and the exterior noise of the operating wind turbines. Ultimately this effort is meant to provide the wind turbine industry with a tool for the design of more efficient wind turbines.
2. Computational Fluid Dynamics Model
The CFD model of the blade was developed using a commercial CFD code, Fluent 6.2
[14]
. Fluent was chosen due to its validated performance and availability of physical models and widespread presence in the power generation industry. The blade geometry was reconstructed in the mesh generator, Gambit 2.3, from the data provided in Giguère and Selig
[17]
. The 3D simulations were run using the fully 3D unsteady laminar NavierStokes (NS) equations. Note that while the flow is turbulent, unsteady laminar simulations were used in consideration of the fact that the present work presents a methodology and develops a tool for the aeroelastic analysis of slender wind turbine blades which is independent of the unsteady nature of the flowfield. That is, the same methodology and workflow philosophy that can be applied to a laminar flow can also be used for a turbulent flow, and this flexibility is part of the appeal of the present method. The intent here is to provide the reader with a recipe on the procedure to follow in wind turbine rotor fluidstructure simulations, which has a general value, independently of the qualitative nature of the CFD simulations. In future work the extension to a turbulent fluid model will be presented. The 2D simulations were run using the unsteady Reynolds Averaged NavierStokes (URANS) equations and turbulence effects were computed using the Wilcox
kw
turbulence model
[10
,
14]
. In both cases the flow was considered to be
Blade and Rotor dimensions
Blade and Rotor dimensions
incompressible and a pressurebased solver was employed. Secondorder upwind spatial discretization was used for the transport equation and a secondorder RungeKutta scheme was used for time discretization
[14]
.
 2.1 Boundary conditions
The geometry of the wind turbine and the definition of the boundary conditions with a uniform velocity inlet profile, allows the use of rotationally periodic boundary conditions so that only one single blade has been simulated. The blade geometry is shown schematically in
Fig. 1
and the blade dimensions are summarized in
Table 1
. The blade twist is plotted together with the nominal angleofattack for a rotational speed of 72 rpm in
Fig. 2
.
Additional geometrical parameters are provided in the numerical results, Section 4. The computational domain, shown in
Fig. 3
, comprises a section of a cylinder 120° in arc. The radial dimension is 30 m and the domain extends 15 m in front and aft of the rotor so as to minimize any influence from the boundary conditions. The operating conditions for the simulations are summarized in
Table 2
. A uniform steady velocity profile of 6.0 m/s was used for the upstream boundary condition while a constantpressure boundary condition was used for the downstream boundary condition
Summary of the flow conditions, ISA sealevel
Summary of the flow conditions, ISA sealevel
Twist and nominal angle of attack distribution of the blade for a pitch angle of 0 degree. The nominal blade pitch angle is given by the twist angle at r/R=0.75.
and for the side at the plane of maximum radius. An axis boundary condition was used for the axis of rotation of the rotor and rotationallyperiodic boundary conditions were applied at the two side planes.
 2.2 Computational grid
The computational mesh is a hybrid mesh that contains both hexahedral and tetrahedral elements for a total of 1.14M cells. The hexahedral elements are primarily used to construct a finer mesh on the solid surfaces of the turbine so as to capture the viscous effects of the boundary layer. A canonical grid convergence study was not performed for the present work although solutionsteered grid adaption was performed to study the sensitivity of the force and momentum results to mesh refinement which was found to be minor for the present conditions.
 2.3 2D/3D Aerodynamic coupling
In this work 2D and a 3D aerodynamic models have been compared. The 2D model enables one to calculate the inflow ratio, denoted as
through a combined bladeelement and momentum theory. It is possible to write
[15
,
23]
where
(see
Table 1
);
In the previous equation,
is a function of the local spanwise distance
x
. The evaluation of the inflow ratio is one of the critical issues when dealing with the aeroelastic response of the rotating blade. In this paper, the 3D aerodynamic computations have been used to estimate the actual inflow ratio. Indeed, the results obtained from the CFD computations have been used in the same way the experimental aerodynamicist would use if the equivalent 2D coefficients have to be estimated from inflight or wind tunnel tests. This approach has the great advantage in reducing the time required for the aeroelastic computations since the time consuming CFD simulations are run only once. Following
[15]
and
[23]
, the inflow ratio, as function of the spanwise location
could be estimated from the knowledge of the effective angle of attack gained from the 3D CFD. The sequence to evaluate the inflow ratio for the 3D case is as follows:
First, the effective angle of attack,
α_{eff}
, is estimated from Eq. (2a) considering the 2D inflow ratio as the starting point, then
is evaluated from Eq. (2b). At this point, the values of
and
are different, therefore, Eqs. (2a) and (2b) are iteratively solved until Eq. (2c) is satisfied. Note that ε is a small arbitrary quantity and for the present simulations it was set to 1x10
^{5}
. The coefficients
c
_{1}
and
c
_{2}
are extracted via a linear interpolation from the FLUENT threedimensional results.
3. Structural Dynamics Model
Assuming the structural behavior of the windmill blade equivalent to a rotating beam, the corresponding mathematical model has been derived from the nonlinear flaplagtorsion equations of motion presented in Hodges et al.
[11]
. In this representation, the blade is a cantilever, straight, slender, homogeneous, isotropic with nonuniform properties, with twist and both mass and tensile offsets. The resulting equilibrium equations are further simplified neglecting those terms considered to be of the third order with respect to the bending slope, that is an arbitrary small parameter, and not contributing to damping. Moreover, in this set of equations the radial displacement is not explicitly considered since the corresponding local tension is used. Therefore, the bending deflections are not responsible for the windmill blade axial extensions that are simply geometrical consequences of the transverse bending deflections
[19]
. Taking into account the previous assumptions, the dynamic system is represented by a set of coupled nonlinear integropartial differential equations suitable for describing significant deflections
[19
,
20]
of the windmill blades where the unknown state is represented by the inplane displacement of the elastic axis,
v
(
x
,
t
), the outofplane displacement of the elastic axis,
w
(
x
,
t
), and the crosssection elastic torsion deflection,
ϕ
(
x
,
t
).
The aerodynamic load, used to couple the previous equations of motion, is based on the lowfrequency approximation of twodimensional Greenberg theory
[12]
, leading to a simple quasisteady striptheory capable to take into account a pulsating freestream for an unsteady aerodynamic loads acting on a thin airfoil. In this formulation, a first coupling with the 3D aerodynamic calculations is introduced. The local section aerodynamic parameters are evaluated from the results of the CFD calculations described in the previous section. Moreover, the aeroelastic model is completed by including the effects of the flow induced by the wake on the blade angleofattack as resulting from Eqs. (2a) (2c). Again, considering the aerodynamic 3D calculations as an experimental test bed, the changes of the angleofattack due to the induced velocity could be estimated from the velocity changes of the flow through the actuator disk, as reported in the previous section and in
[12]
and
[13]
.
The resulting final aeroelastic system could be integrated in the space domain using the Galerkin approach adopting the nonrotating modes of the blade as shape functions. If
x
is the vector of the modal amplitudes (Lagrangean coordinates),
are the linear timeinvariant mass, damping and stiffness structural and aerodynamic contributions, and
represents the set of all nonlinear and/or timedependent coefficient terms, then the set of the nonlinear ordinary differential equation is written as
where
ψ = Ωt
is the azimuthal position of the blade, with
Ω
as the angular velocity of the rotor, and
θ =θ(ψ)
is the pitch angle. Since the periodic solution is sought, then the harmonic balance approach has been used to solve these equations in the time domain.
The dynamic hub loads (root shears and moments) corresponding to the aeroelastic deformation, achieved from the previous solution, could be evaluated in the rotating frame. Then, the hub loads in the fixed frame of reference could be achieved by adding the contributions from all the blades expressed in the inertial frame of reference. If the drag, radial, and vertical shear forces at the rotor hub are denoted as
S_{x}
,
S_{r}
and
S_{z}
, whereas the root bending and torque moments are
N_{F}
and
N_{L}
respectively, then the total forces and moments in the nonrotating frame for a N bladed windmill are given in
[6]
,
[23]

[26]
:
Computational domain and boundary conditions.
4. Numerical Results
In this section, the results from both the aerodynamic and aeroelastic simulations are reported. In the following subsection, after a general description of the wind turbine model used in the numerical analysis, the results from the aerodynamic simulation of the rotating S809 airfoil are first outlined, then the results from to the aeroelastic simulations are presented.
 4.1 Main characteristics of the NREL S809 wind turbine
The work focuses on a small windturbine blade designed to operate on a stallregulated downwind wind turbine having a rated power of 20 kilowatt.
The blade geometry and dimensions were extracted from the National Renewable Energy Laboratory (NREL) experiment setup
[17
,
18
,
21
,
22]
, and some of these parameters are reported for completeness in
Table 1
. The rotor diameter is 11.06 m with a root extension from the center of rotation to the airfoil transition of 0.88m, see
Fig. 4
. The blade is twisted (root twist of 20°, tip twist of 2.5° about the 30% chord location) and tapered (root chord of 0.737 m, tip chord of 0.305 m) and it features a NREL S809 airfoil with a thickness ratio of 21% (see
Figs. 4
and
5
)
[21]
. The blade pitch ratio is referenced to the 75% span location, see
Fig. 4
, and it is defined with respect to the x and y axes, as shown in
Fig. 5
(c). The threebladed rotor has a solidity of 0.069 and it rotates at 72 rpm with a uniform inflow speed of 6 m/s, giving a tipspeedratio of 7.0. An ellipsoid of revolution is located between the blade root extension and the axis of revolution of the rotor to maintain a smooth and realistic flow field
Blade planform (not to scale).
Schematics of the blade used in the present study
around the root region.
 4.2 Aerodynamic simulation
This section presents the results of the aerodynamic modeling of the rotor. First, a validation study of an S809 airfoil was performed to demonstrate the validity of the CFD approach to calculate the loads on the blade structure. In
Fig. 6
the pressure coefficient distribution along the chord of the wind turbine blade is presented achieved with the reference flow conditions reported in
Table 2
.
The predicted C
_{p}
distribution is in overall excellent agreement with the experimental data with the exception of the location at 55%c where a mild flow separation is present and the CFD solution fails to predict its strength. The numerical solution shows noticeable oscillations in the C
_{p}
distribution in the region ahead of 55%c. The oscillations are indicative of an incipient flow separation on the upper surface of the airfoil. Comparisons for the drag, lift and
Comparison of moments and forces for the 2D calculations.
Comparison of moments and forces for the 2D calculations.
moment coefficient of the present work to experimental results
[21]
and the Wolfe and Ochs CFD simulations
[22]
, where the
kε
turbulence model was used, are presented in
Table 3
.
The comparison shows some discrepancy between FLUENT
^{®}
results and the experiments reported in Somers
[21]
, especially for the drag coefficient which is 37% larger than the measured value. One source of error is the lack of detailed information on the inlet boundary conditions used in the experiment. In the present CFD simulations a uniform inlet velocity profile with a freestream velocity of 48.5 m/s was used to match the experimental Reynolds number of 2.0x10
^{6}
. Additional source of error is the use of a viscous laminar flow model even at high Reynolds number where turbulence flow might significant effect, with potential effect also from turbulence boundary layer and trailing edge separation. These are all issues of great importance when looking into accuracy of the fluidflow numerical simulations. Inlet turbulence was set at a turbulence intensity, TI=k/(u
_{∞}
)
^{2}
, of 5% and at a turbulent viscosity ratio,
μ_{T}/μ
, of 10. Note that the 2D structured mesh consists of rectangular elements with a y
^{+}
of less than 1.0 near the solid wall of the airfoil. The 3D simulations of the blade were run at five different blade pitch angles, 0°, 2°, 6°, 12° and 18°. The pitch angle is defined with respect to the airfoil at 75% of the blade radius, as shown in
Fig. 7
(b). For completeness and to ensure the overall similarity in order of magnitude between the turbulent and laminar results, a comparison of a laminar vs. a turbulent simulation with the Wilcox
kw
model was run for the full
Comparison of the experimental and computational pressure coefficient distribution on the airfoil for the conditions as in Wolfe and Ochs [20].
threedimensional blade for a nominal angle of attack of 0 deg. The results are presented in
Fig 6
a.
While the absolute forces produced by the laminar and turbulent simulations are different, as expected, there is an overall similarity in the results, especially for the Cd. Also, it is important to notice that the methodology here proposed is valid for both flow regimes. Therefore, the 3D simulations were run without turbulence modeling. Also, laminar conditions might exist in some specific situations but in general it will be up to the designer to check the flow regime and run the appropriate flow model.
It is worth mentioning that other realistic computational techniques, such as the striptheory, could be adopted in conjunction with 2D turbulence computations. As an example, the performance of several 2D turbulence computations at a number of axial positions along the blade could be performed and the formulation of an approximate 3D loading for the full length of the blade would then be a feasible option. This was however well beyond the scope of the present paper and it is left to the interested reader to pursue this or other viable approaches.
Comparison of the laminar and turbulent solutions for the full threedimensional blade at an angle of attack of 0deg
(a) Tip vortices convected downstream. (b) Vorticity contours on a longitudinal plane showing the core of the wingtip vortices being con vected downstream and dissipated. Side view along blade centerline. Blade pitch angle = 6°.
Figure 8
shows the momentum deficit of the incoming flow once it has passed the rotor. In this isometric view the rotor is shown in full and the hub is included for clarity. The relative highmomentum in the undisturbed region away from the rotor and in the small region of undisturbed flow near the axis of rotation is clearly noticeable. The latter highlights the action of the turbine blades that are extracting linear momentum from the incoming air flow to turn it into rotor rotational momentum. The use of rotationallyperiodic boundary conditions in the 3D simulations allowed the study of the bladewake interaction. However, from monitoring the force time history, it appears that for a rotor angular speed of 72 rpm and an undisturbed airspeed of 6 m/s (tip speed ratio of 7.0) the interaction is almost negligible.
The time history of the lift, drag and moment coefficient on the blade,
Fig. 9
, show a Strouhal number based on a blade average chord of 0.521 m and an average blade velocity of 26.7 m/s of 0.22, corresponding to a dimensional shedding period of 0.09 s, also confirmed by a Fast Fourier Transform analysis of the force data. This period is shorter than that of the wake shed by the preceding blades, 0.278 s
Axial velocity contours.
and it is attributable to the natural shedding frequency of the blade. Notice also that the amplitude of this oscillation is relatively small, 0.27% of the mean value for the drag and moment coefficient and 0.16% of the mean value for the lift coefficient. The bladewake interaction is also visualized in
Fig. 7
(a) where the streamlines originating from the blade tips are plotted as they are convected downstream. The picture shows that the blades are well clear of the wake core produced by the preceding blade.
Fig. 7
(b) maps the vorticity contours on a flat plane normal to the blade axis at r/R = 0.75. The mapping shows the vortex chores from the main blade whose cross section is visible at the center of the map. The vortex chore coming from this blade, blade 1, is visible traveling downstream in the upper center of the figure. The vortex chore produced by the preceding blade is also visible before and after it approaches the azimuthal location of blade 1. Again, the mapping highlights the weak interaction of the blades with the wake.
The lift, and drag coefficients for the blade at different pitch angles are plotted in
Fig. 10
as a function of the blade span. The coefficients are computed using the dynamic pressure at the specific blade section and the angle formed by the local velocity vector. The plot of the lift coefficient shows a consistent trend at all pitch angles of decreasing
Periodicity in the blade lift, drag, and moment coefficient as a function of nondimensional time. C_{Avg} and U_{Avg} correspond to the arithmetic mean of the blade chord (0.521 m) and of the total relative velocity (26.7 m/s). In this case the shedding period correspond to a Strouhal number of 0.22 (t= 0.09 s).
Lift and drag coefficient along the blade span at different pitch angles as a function of blade radial section.
lift with span location. Also, the lift tends to decrease as the pitch angle is increased. The drag coefficient curves indicate two distinct trends depending on the pitch angle. At lower pitch angles, the drag increases steadily up to 50% of the span where it abruptly decreases, possibly indicating separation. This discontinuity decreases in intensity with increasing pitch angle and it disappears completely for a pitch angle of 12° and 18°.
Note that, at this location for the low pitch angles there is a minor drop in lift coefficient. According to the orientation of the axes or reference, a negative drag coefficient corresponds to a “forward” thrust for the rotor, the highest power is obtained for a blade pitch angle of 6 . Lower pitch angle settings are not producing as much thrust since most of the aerodynamic force is transformed in lift which is not contributing in a major way to the production of usable power. Also note that the overall thrust is given by a combination of lift and drag.
Figure 11
shows the liftcurve of the different blade sections as a function of angleofattack. The lift curve slope is 3.00 per radiant (0.95
π
/Rad) and all the blade sections show the same behavior and lift values, as expected. The liftcurve slope is much lower than the potential flow value of 2
π
/Rad and also of the 1.91
π
/Rad value calculated by Wolfe and Ochs
[21]
.
The lower value of the 3D solution is due to the threedimensional effects of the blade. These effects were calculated following Bramwell
[23]
where the induced velocity,
U_{i}
is calculated from
Where
λ_{i} = U_{i}/ωR;λ_{c} = U_{∞}/ωR;σ = Nc/πr; a'
is the liftcurve slope (3.00/Rad in this case) and
N
is the number
Lift coefficient along the blade span at selected pitch angles as a function of blade radial section.
Plot of induced velocity and induced angle of attack as a function of the blade span.
Comparison between the different induced velocity models.
of blades of the rotor. The calculated induced velocity is shown in
Fig. 12
together with the induced angleofattack. As expected the induced velocity is higher closer to the hub where the blade loading is generally higher. Also, the induced velocity is higher for the case with higher pitch angle and tends to decrease as the pitch angle is lowered.
The coefficients
c
_{1}
and
c
_{2}
describe the relationship between
λ_{i}
and
α_{eff}
listed in
Table 4
. The results obtained
Fluent data for 3D model
from the previous analysis are shown in
Fig. 13
; from these figures the different trend that occur for both
λ_{i}
and
α_{eff}
, particularly for low values of the local spanwise position
r
, can be underlined. At the first step, the value of the effective angleofattack, starting from the 2D inflow ratio value, is calculated; The second step consist of calculation the 3D inflow ratio through the existing relationship between each couple of
λ_{i}
and
α_{eff}
. This is done for every station x along the span with variable
α_{eff}
; in the last step, the difference between the 3D value calculated in the previous step and the first tentative 2D value in step one is computed: if the difference is greater than a fix tolerance, ε = 10
^{4}
, the steps are repeated in a loop, in the first step with a new value of
λ_{i}
is used corresponding to the 3D value previously estimated.
 4.3 Structural blade characterization
As stated in the previous section, the nonrotating modes of the blade are required in order to solve the aeroelastic system
Comparison between the reference and the actual modal model
Comparison between the reference and the actual modal model
Inertial and Stiffness distribution of the windmill blade.
Comparison between the static displacements (ψ =0).
of equations. The flapwise and lag stiffness together with the mass distributions, required for setting up the modal model, have been evaluated using both the data reported in Giguère and Selig
[17]
, and estimated from similar blades by using an approach based on the dynamic similitude. In
Fig.14
the mass and the flapwise stiffness distribution as function of the blade span are reported, while, since no torsional stiffness data were available from
[17]
, a typical value for this kind of blade, considered constant over the blade span, has been assumed. These structural characterizations lead to natural frequencies reported in
Table 5
. As one can see, the considered modal model is almost the same as the one identified in
[17]
.
In
Fig. 14
the mass and the flapwise stiffness distribution as function of the blade span are reported, while, since no torsional stiffness data were available from
[17]
, a typical value for this kind of blade, considered constant over the blade span, has been assumed. These structural characterizations lead to natural frequencies reported in
Table 5
. As one can see, the considered modal model is almost the same as the one identified in
[17]
.
 4.4 Static Aeroelastic Simulation
The aeroelastic simulation has been performed considering the windmill rotating at fixed operational conditions in terms of rotational speed and blade aerodynamic incidence. This condition is equivalent to consider the windmill operating at constant, ideal, wind flow regime. The solution of the system of aeroelastic equations is carried out through harmonic balance leading to the evaluation of the elastic deformation of the hingeless blade. Herein, the aerodynamic coefficients and the inflow distributions along the blade span due to induced velocity from the 3D calculations are taken into account. The static component of the flap and torsional deflections, for the blade azimuthal position ψ=0, are presented in
Figs. 15
and
16
.
It is noteworthy to observe that the flap deflection occurs in the wind direction (
Fig. 15
a), whereas the lag occurs in the direction of the rotational velocity (
Fig. 15
b). By comparing these results with the ones of a rotorcraft in hover with ascension velocity equal and opposite to the one of the wind
Comparison between the torsional deflections (ψ=0).
Comparison of the different estimates of the rotating frame reference shears (ψ=0).
Comparison of the different estimates of the rotating reference frame moments (ψ=0).
mill it can be noticed that the flap and lag are the opposite in directions; these results are consistent to the fact that a torque is produced at the hub of the wind mill by the interaction windblade system, while in the helicopter a torque at the shaft is provided by the propulsion system to generate the thrust. It is also clearly evident from these figures that the 2D model does predict lower flapwise deflections that the 3D model counterparts which is associated with the slightly increased effective angle of attack that is present in the 3D case, see
Fig. 13
a. On the other hand, the corresponding lagwise deflection is slightly lower, which might be associated with an overestimate on the 2D aerodynamic loading. The different is quite small also considering also
Sketch of the Rotating and Fixed Reference Frames
that in this direction the loads are one order of magnitude smaller and in the flapwise direction. The combined effect of flap and lag deflection contributes to a large discrepancy in the torsional direction, as it is displayed in
Fig. 16
. Clearly the inherent flexibility of the rotor combined with the complex aerodynamic characteristics in the spanwise direction creates some complexities in the deformation field and it is quite challenging to fully interpret. It is also possible to estimate the shear and moments along the span of the blade in the undeformed rotating reference frame attached to the blade, as reported in
Fig. 17
.
The comparison between the dynamic loads as function of the blade span, considered at ψ=0, are reported in
Figs. 18
and
19
where the different evaluations of the dynamic
Comparison of the different static hub load estimates in the rotating reference frame.
Comparison of the different harmonic hub load estimates in the rotating reference frame.
Comparison of the different harmonic hub load estimates in the rotating reference frame.
rotating shears and moments, are depicted. Moreover, the different hub load estimates, always in the rotating reference frame, are reported in
Figs. 20
through
22
. From
Fig. 20
, only a remarkable difference in the static outofplane shear is reported, whereas a completely different behavior, but with much smaller amplitudes, of the dynamic elastic reactions (acting at 3/rev and 5/rev harmonics), is evident in
Figs. 21
and
22
, respectively.
When considering the resulting loads acting in the fixed reference frame, the convention reported in
Fig. 17
are used. For the positive values of the hub dynamic loads: the thrust, T, is developed by the rotor acting in the orientation of the wind; the torque, Q, (responsible for the power generation of the wind turbine) has the orientation of
Ω
; H and Y are the shears acting in the x and y reference directions, respectively; and, finally, M
_{x}
and M
_{y}
are the bending moments in the x and y axis, respectively. The aeroelastic loads evaluated in the fixed reference frame excite the windmill structure at discrete frequencies that are the number of blades of the rotor times the rotating angular frequency (see e.g.
[20]
through
[22]
,
[27]
) plus the static component.
The different evaluations of the static loads are reported in
Fig. 23
, whereas
Fig. 24
refers to the 3/rev harmonics loads.
For the lateral equilibrium, the rotor forces H, Y and the moments M
_{x}
, M
_{y}
should be approximately zero; indeed this is the case as represented in
Fig. 23
. It appears that the 2D and 3D aerodynamic models predict dissimilar results, in particular between the force T and moment Q, being the value of the 2D model greater than the one of the in 3D model. The different evaluation of the dynamic loads leads to an over estimate of the 2D static thrust with respect to the one predicted by a 3D analysis. This fact could be of great interest in designing the windmill tower since the resulting bending moments, evaluated at the basement of the tower, are significant different, i.e., the 3D simulation provides lesser stressed tower structure than the 2D analysis. Furthermore, the 3D analysis predicts a value of moment Q that is greater than the correspondent 2D value, resulting in the capability of the windmill turbine to extract higher power.
5. Conclusions
Unsteady aerodynamics and aeroelastic characteristics of a mediumsized windturbine blade operating under ideal conditions have been considered. The tapered/twisted blade was designed to operate on a stallregulated downwind wind turbine originally studied in an experiment setup at the NREL. The structural dynamics of the flexible windturbine blade has been described by a nonlinear flaplagtorsion slenderbeam differential model. The aerodynamic quasisteady forcing terms needed for the aeroelastic governing equations
Comparison of the different static hub load estimates in the fixed reference frame.
Comparison of the different harmonic hub load estimates in the fixed reference frame.
have been predicted through a striptheory based on a simple 2D model, whose pertinent aerodynamic coefficients have been obtained using CFD. Aerodynamic loads are computed by CFD and analytical techniques. Turbulence effects in the 2D simulations run for validation purposes were accounted by the Wilcox
kw
model. In this first approximation of the aerodynamic loading, the 3D CFD simulations were run using the unsteady laminar NavierStokes equations solved via a commercial software, FLUENT
^{®}
. Five blade pitch angles were considered and their aerodynamic characteristics compared. The 6° pitch angle was found to be producing the highest power setting due to the relatively high (propulsive) drag force. For pitch angles higher than 6° the blade showed an increased drag coefficient indicating lessthanoptimal performance. The interaction of the rotor blades with their own wake was found to be negligible for the operating conditions under consideration due to the convective effect of uniform airflow. The resulting unsteady hub loads are achieved by a first space integration of the aeroelastic equations by applying the Galerkin’s approach and by a time integration using a harmonic balance scheme. Two and three dimensional computations for the unsteady aerodynamic loads, the flap, lag and torsional deflections, and the forces and moments applied at the hub of rotor along with aeroelastic simulations are illustrated. Structural results have shown that the responses to the 2D and 3D aerodynamic models produce dissimilar results, therefore stressing the importance to keep the dimensionality effects into consideration if quantitative results are expected from the simulations. To this point the authors acknowledge that turbulence models are very problematic for these flows, as the flow typically transitions at some point on the airfoil, and none of the usual turbulence models currently available handle very well the adverse pressure gradient, that occur on the aft of the airfoil, very well. Therefore experimental data will always be needed to back up computational predictions. Clearly neglecting flow separation along with turbulence, especially in 3D simulations can lead to inaccurate results. In particular, 3D analysis predicts a value of a torque, and consequently the power output, that is a conservative estimation of the energy that can be extracted from the wind. However, since the main goal of this paper is not to provide an highly accurate aerodynamic model of the wind turbine during operational condition, but rather to demonstrate the proofofconcept of an innovative way to perform reliable and fast flowstructure simulations for aeroelastic design purposes, the assumptions made are reasonable and can be easily extended to more realistic aerodynamic models. Also the design of the wind turbine tower could benefit from the results of the 3D simulations performed since a lesser stressed tower is predicted. Future work will focus on the turbulence effects on the 3D simulations using unsteady RANS equations and DES techniques, and on the inclusion of an atmospheric boundary layer and asymmetric hub loading.
Acknowledgements
The authors would like to thank the Metropolitan Development Association of Syracuse and Central New York, Inc and The National Science Foundation, grant NSFCMMI1031036, for providing partial support for this research. Also, partial support for the present work was provided by the University of Rome, “La Sapienza” through the project “Progetto Giovani Ricercatori: Dinamica ed Aeroelasticita dei Rotori di Elicottero”.
Burton T.
,
Sharpe D.
,
Jenkins N.
,
Bossanyi E.
2001
Wind Energy Handbook
1st Edition
John Wiley & Sons, Ltd.
Chichester
ISBN: 0471489972
Manwell J. F.
,
McGowan J.G.
,
Rogers A.L.
2002
“Wind Energy Explained”
1st edition
John Wiley & Sons, Ltd.
Chichester
ISBN 0471499722
Eggleston D.M.
,
Stoddard F.S.
1987
“Wind Turbine Engineering Design”
Van Nostrand Reinhold
New York
307 
Germanischer Lloyd
1993
Rules and regulations, VI  Nonmarine Technology, Part 1  Wind Energy
Hamburg
Librescu L.
,
Marzocca P.
,
Silva W. A.
2005
“Aeroelasticity of 2D lifting surfaces with timedelayed feedback control”
Journal of Fluids and Structures
20
(2)
197 
215
DOI : 10.1016/j.jfluidstructs.2004.10.005
Behal A.
,
Rao V. M.
,
Marzocca P.
,
Kamaludeen M.
2006
“Adaptive Control for a Nonlinear Wing Section with Multiple Flaps”
Journal of Guidance, Control, and Dynamics
29
(3)
744 
748
DOI : 10.2514/1.18182
Leishman J.G.
2000
“Principles of Helicopter Aerodynamics”
1st Edition
Cambridge University Press
ISBN 0521523966
Leishman J. G.
2002
“Challenges in modeling the Unsteady Aerodynamics of Wind Turbines”, AIAA 20020037
21st ASME Wind Energy Symposium and 40th AIAA Aerospace Sciences Meeting
Reno, NY
Politis E.S.
,
Nikolaou I.G.
,
Chaviaropoulos P.K.
,
Bertagnolio F.
,
Sørensen N.N.
,
Johansen J.
2005
“KNOWBLADE Task4 report; NavierStokes Aeroelasticity”, RisøR1492 (EN)
Risø National Laboratory
Roskilde, Denmark
Wilcox D.C.
1998
“Turbulence Modeling for CFD”
DCW Industries, Inc.
La Canada, California
Hodges D. H.
,
Dowell E. H.
1974
“Nonlinear Equation for the Elastic Bending and Torsion of Twisted Non Uniform Rotor Blades”
NASA TN D7818
Greenberg J.M.
1947
“Airfoil in Sinusoidal Motion in a Pulsating Stream”
NACA TN1326
Theodorsen T.
1935
“General Theory of Aerodynamic Instability and the Mechanism of Flutter”
NACA Report No.496
Fluent, Inc.
2006
FLUENT User Guide
September 29
Stepniewski W.Z.
,
Keys C.N.
1984
“Rotary – Wing Aerodynamics”
Dover Publications, Inc.
New York
Johnson W.
1980
“Helicopter Theory”
Dover Publications
Giguère P.
,
Selig M. S.
1999
“Design of a Tapered and Twisted Blade for the NREL Combined Experiment Rotor”
National Renewable Energy Laboratory
Golden, CO
NREL/SR50026173
Butterfield C.P.
,
Musia W.P.
,
Simms D.A.
1992
“Combined Experiment Phase I Final Report”
National Renewable Energy Laboratory
Golden, CO
NREL/TP2574655
Hodges D.H.
,
Ormiston R.A.
1976
“Stability of Elastic Bending and Torsion of Uniform Cantilever Rotor Blades in Hover with Variable Structural Coupling”
NASA TN D8192
Friedmann P.P.
1990
“Helicopter Rotor Dynamics and Aeroelasticity: Some Key Ideas and Insights”
Vertica
14
(1)
101 
121
Somers D. M.
1989
“Design and Experimental Results for the S809 Airfoil”, Technical Report
Airfoils, Inc.
State College, PA
Wolfe W.P.
,
Ochs S.S.
1997
“CFD Calculations of S809 Aerodynamic Characteristics”
AIAA Paper 970973
Bramwell A.R.S.
1976
“Helicopter Dynamics”
Butterworth Heinemann Ltd
June 1
Bir G.S.
2005
“Structural Dynamics Verification of Rotorcraft Comprehensive Analysis System (RCAS)”
National Renewable Energy Laboratory
Golden, CO
NREL/TP50035328
Bielawa R.L.
1992
“Rotary Wing Structural Dynamics and Aeroelasticity” AIAA Education Series
AIAA
Coppotelli G.
,
Cellini F.
“A Rotorcraft Trim Procedure for the Prediction of Fuselage Vibrations” Advances in Engineering Research, Vol. 4
Nova Science Publisher
NY, USA
ISBN: 9781621006954
235 
257
Balis Crema L.
,
Coppotelli G.
,
Grappasonni C.
2011
“On the Effects of Structural Modifications in the Large Wind Turbine Dynamics”
Proceedings of the 52nd AIAA/ASME/ASCE/ASC Structures, Structural Dynamics and Materials Conference Denver (CO)  USA
Denver (CO)  USA