Post on 12-Mar-2018
Technische Universität Berlin
Fakultät V
Institut für Strömungsmechanik und Technische Akustik
Fachgebiet für Experimentelle Strömungsmechanik
Diploma Thesis
A BEM Based Simulation-Tool for Wind
Turbine Blades with Active Flow Control
Elements
from : Guido Weinzierl
Date : 19.04.2011
Course : Strömungslehre II
1st Supervisor : Prof. Dr.-Ing. C. O. Paschereit
2nd Supervisor : Dipl.-Ing. G. Pechlivanoglou
Preface
This diploma thesis was written in order to fulll the requirements of obtaining the
degree Dipl.-Ing. at Berlin University of Technology. The work was carried out in
collaboration with Smart Blade GmbH and the Wind Energy Group at the Institute of
Fluid Dynamics and Technical Acoustics (ISTA) at the Berlin University of Technology.
I especially want to thank Georgios Pechlivanoglou for his great support and supervi-
sion of the project. Many thanks as well to Oliver Eisele from Smart Blade GmbH, for
the numerous brainstorming sessions and Smart Blade GmbH in general, for generously
funding this thesis.
Abstract
This thesis describes the development of a software tool which provides a method to
investigate the use of dierent active ow control (AFC) concepts for load reduction
and power regulation of wind turbines. The software features an aeroelastic model
to calculate the dynamic response of the wind turbine structure. The program is an
extension of QBlade, an open-source GUI application for wind turbine calculations.
The user can easily dene a wind turbine blade on which various active elements
can be positioned. The dierent aerodynamic characteristics of the AFC-elements
are considered by their individual lift and drag polars. These polars can either be
calculated using an implemented two-dimensional panel method code (XFoil) or im-
ported to provide an interface to wind tunnel measurement data or CFD calculations.
To model the aerodynamic and structural behavior of the turbine, a binding to the
aerodynamic analysis routines AeroDyn and the structural analysis code YawDyn is
implemented. These simulation codes are provided by the National Wind Technology
Center (NWTC) of the National Renewable Energy Laboratory (NREL).
In order to control the active elements on the blade, two control approaches are
provided: a simple optimization loop, to nd an optimal actuator position for each time
step and a PID controller. Both approaches can be used, for example, to minimize the
root bending moment of the wind turbine blades, either by keeping local blade element
forces constant or by minimizing blade deections or blade deection rates.
Zusammenfassung
Die vorliegende Diplomarbeit beschreibt die Entwicklung einer Software, die es er-
möglicht den Einsatz von Elementen zur aktiven Strömungskontrolle (AFC) zur Las-
treduktion und Leistungsregelung an Windkraftanlagen zu untersuchen. Die Software
beinhaltet ein aeroelastisches Simulationsmodul, um den dynamischen Einuss der
AFC-Elemente auf die Struktur der Windturbine zu berechnen. Das Programm ist
eine Erweiterung von QBlade, einer open-source Anwendung zur Rotorblattentwick-
lung und -berechnung nach der Blatt-Element-Methode.
Der Benutzer kann auf einfache Art und Weise ein Rotorblatt entwerfen, auf dem
mehrere aktive Elemente platziert werden können. Der unterschiedliche aerodynamis-
che Einuss der verschiedenen Elemente, ist durch ihre individuellen Auftriebs- und
Widerstandspolare gekennzeichnet. Die Polaren können dabei entweder direkt über
eine eingebaute zweidimensionale Panelmethode (Xfoil) berechnet werden, oder sie
werden importiert, um die Verbindung zu Windkanalversuchen oder CFD Simula-
tionen herzustellen. Um das aerodynamische und strukturelle Veralten der Anlage
zu modellieren, werden die aerodynamischen Berechnungsroutinen AeroDyn und der
strukturdynamische Berechnungscode YawDyn benutzt. Die Programme werden vom
National Wind Technology Center (NWTC) des National Renewable Energy Labora-
tory (NREL) bereitgestellt.
Um die aktiven Elemente auf dem Rotorblatt zu regeln, werden zwei Herange-
hensweisen verfolgt: Zum einen eine einfache Optimierungsschleife, um für jeden
Berechnungszeitschritt die optimale Aktuatorposition zu nden, zum anderen ein PID
Regler. Beide Regelstrategien können beispielsweise dazu genutzt werden, das Blat-
twurzelbiegemoment zu reduzieren. Dazu werden entweder die lokalen Kräfte am Blat-
telement konstant gehalten, oder die Rotorblattbiegung oder deren Änderungsrate
minimiert.
Contents
1 Introduction 1
I Model Theory 4
2 Aeroelastic model 5
3 Aerodynamics 8
3.1 Wake modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.1.1 Blade element momentum theory (BEM) . . . . . . . . . . . . . 10
3.1.2 Generalized dynamic wake model (GDW) . . . . . . . . . . . . 14
3.2 Airfoil aerodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2.1 2D static airfoil characteristics . . . . . . . . . . . . . . . . . . . 15
3.2.2 Polar extrapolation . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2.3 Unsteady blade element aerodynamics . . . . . . . . . . . . . . 22
3.2.4 Stall delay and 3D eects . . . . . . . . . . . . . . . . . . . . . 26
3.3 Tower shadow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 Structural dynamics 28
4.1 YawDyn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5 Turbulent wind simulation 34
6 Active Flow Control 36
Contents VI
II Software 39
7 QBladeAE 40
7.1 Active Flow Control simulation . . . . . . . . . . . . . . . . . . . . . . 41
7.1.1 Optimization loop . . . . . . . . . . . . . . . . . . . . . . . . . . 43
7.1.2 PID controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
7.2 Blade related simulation parameters . . . . . . . . . . . . . . . . . . . . 49
7.2.1 QBlade and NREL blade format . . . . . . . . . . . . . . . . . . 49
7.2.2 Blade mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
7.2.3 Blade center of gravity . . . . . . . . . . . . . . . . . . . . . . . 50
7.2.4 Blade mass moment of inertia . . . . . . . . . . . . . . . . . . . 52
7.2.5 Torsional root spring constant . . . . . . . . . . . . . . . . . . . 52
7.2.6 Dynamic stall parameters . . . . . . . . . . . . . . . . . . . . . 52
7.3 Program modules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
7.3.1 Blade design with active elements . . . . . . . . . . . . . . . . . 54
7.3.2 Aerodynamic representation of active elements . . . . . . . . . . 56
7.3.3 Wind eld simulation . . . . . . . . . . . . . . . . . . . . . . . . 57
7.3.4 Aeroelastic simulation . . . . . . . . . . . . . . . . . . . . . . . 57
III Simulation 60
8 Standard simulation 61
8.1 Turbine and blade model . . . . . . . . . . . . . . . . . . . . . . . . . . 61
8.2 Blade validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.3 Dynamic stall eects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
8.4 Yawed turbine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
8.5 Wind eld . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
8.6 Baseline simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
9 AFC simulation 71
9.1 Flap parameter study . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
9.1.1 Flap positions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
9.1.2 Flap size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
9.1.3 Actuator speed and range . . . . . . . . . . . . . . . . . . . . . 79
Contents VII
9.1.4 Sensor delay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
9.1.5 Multiple aps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
9.2 Optimization loop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
10 Suggestions for future research 85
11 Conclusion 87
Bibliography 88
A Appendix 94
A.1 QBladeAE input les for YawDynAE . . . . . . . . . . . . . . . . . . . 94
A.2 Geometric blade design . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
List of Figures
1.1 Concept of a segmented wind turbine rotor blade with active elements
in form of trailing edge aps [52]. . . . . . . . . . . . . . . . . . . . . . 3
2.1 Local blade element velocities and inow angles. . . . . . . . . . . . . . 6
2.2 Local blade element forces. . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.1 Summary of the various aerodynamic sources that contribute to the
airloads on a wind turbine [32]. . . . . . . . . . . . . . . . . . . . . . . 8
3.2 Rotor of a three-bladed wind turbine with the rotor radius R [18]. . . . 11
3.3 2D arfoil characteristics of a blade element. . . . . . . . . . . . . . . . . 16
3.4 Lift coecient cl(α) at dierent Reynolds numbers and xed/free tran-
sition for the DU 91-W2-250 airfoil. . . . . . . . . . . . . . . . . . . . . 17
3.5 Drag coecient cd(α) at dierent Reynolds numbers and xed/free tran-
sition for the DU 91-W2-250 airfoil. . . . . . . . . . . . . . . . . . . . . 18
3.6 Moment coecient cm(α) at dierent Reynolds numbers and xed/free
transition for the DU 91-W2-250 airfoil. . . . . . . . . . . . . . . . . . 18
3.7 Wind triangular for dierent radial positions [16]. . . . . . . . . . . . . 20
3.8 Exemplary time series of α for two radial positions at rin = 6.8m and
rout = 36.8m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.9 Extrapolated cl for the DU 91-W2-250 airfoil and cl according to the
at plate theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.10 Extrapolated cd for the DU 91-W2-250 airfoil and cd according to the
at plate theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.11 Dynamic stall events on a NACA 0012 airfoil (reprinted from [9] and [32]). 24
4.1 Components of a HWAT structural model . . . . . . . . . . . . . . . . 29
4.2 The equivalent hinge-spring model for the blade ap degree of freedom
[28]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
List of Figures IX
4.3 View of the HAWT dening selected terms and coordinate systems. All
angles are shown in their positive sense. The bold X,Y,Z axes are xed in
space and are the coordinates in which the wind components are dened
(VX, VY, VZ). Note that blade azimuth is zero when the blade is at the
6 o'clock position [28]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.1 Wind speed at hub height and inow velocity at blade tip including
rotation, wind shear and tower eect. . . . . . . . . . . . . . . . . . . . 35
5.2 3D wind eld from QBladeAE with 20x20 points. . . . . . . . . . . . . 35
6.1 Feedback ow control triad (after [25]). . . . . . . . . . . . . . . . . . . 37
7.1 QBladeAE embedded in QBlade and XFLR5. . . . . . . . . . . . . . . 41
7.2 Working principle of QBladeAE with input and output control to the
modied NREL codes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
7.3 Blade with two active elements, which are represented by using several
airfoil polars. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
7.4 Schematic control circuit and control terminology. . . . . . . . . . . . . 44
7.5 Implementation of the optimization loop for nding the optimal polar
for each active section (element dependent). . . . . . . . . . . . . . . . 45
7.6 Implementation of the PID controller: one for each active element (blade
dependent). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
7.7 Exemplary control circuit for the PID controller using a trailing edge
ap as actuator and the blade ap rate as control variable. . . . . . . . 47
7.8 Dierent blade denition in QBlade and NREL format. . . . . . . . . . 49
7.9 Dierent blade masses over blade length and the used exponential ap-
proximation function [51]. . . . . . . . . . . . . . . . . . . . . . . . . . 50
7.10 Simplied geometric representation (rectangular cone) of a homogeneous
blade section for the calculation of the blade center of gravity. . . . . . 51
7.11 Dynamic stall related parameter using clcd-curve for automatically detect-
ing the critical static stall angle αstall. . . . . . . . . . . . . . . . . . . . 53
7.12 QBladeAE active blade design module. . . . . . . . . . . . . . . . . . . 54
7.13 QBladeAE multiple aerodynamic polar module. . . . . . . . . . . . . . 57
7.14 QBladeAE wind eld generator module (beta). . . . . . . . . . . . . . . 58
7.15 QBladeAE wind eld generator module (beta). . . . . . . . . . . . . . . 58
List of Figures X
8.1 3D view of blade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
8.2 Blade tip deection with three dierent root spring stinesses. The wind
inow is steady and constant over the whole rotor disk. . . . . . . . . . 63
8.3 Rotor power over wind speed calculated with QBlade and QBladeAE. . 64
8.4 Blade pitch over wind speed for power regulation. . . . . . . . . . . . . 65
8.5 Inuence of the dynamic stall model on the blade tip deection over
time and the angle of attack over radial position for the blade with a
pitch angle of θb = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
8.6 Inuence of the dynamic stall model on the blade tip deection over
time and the angle of attack over radial position for the blade with a
pitch angle of θb = 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
8.7 Blade tip deection for dierent yaw angles. . . . . . . . . . . . . . . . 68
8.8 Turbulent wind speed time series in x-direction at a hub height of 89m
and a mean wind speed of 13ms. . . . . . . . . . . . . . . . . . . . . . 69
8.9 Results of baseline simulation. . . . . . . . . . . . . . . . . . . . . . . . 70
9.1 Overlapping airfoil contours for positive deection. Red: The original
DU-96-W-180 airfoil; Green: slightly deected exible ap; Red: fully
deected exible ap [41]. . . . . . . . . . . . . . . . . . . . . . . . . . 71
9.2 Overlapping airfoil contours for negative deection. Red: The original
DU-96-W-180 airfoil; Green: slightly deected exible ap; Red: fully
deected exible ap [41]. . . . . . . . . . . . . . . . . . . . . . . . . . 72
9.3 Lift and drag coecient over angle of attack cl(α) for exible ap at
four ap angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
9.4 Extrapolated 360 cl-polar. . . . . . . . . . . . . . . . . . . . . . . . . . 73
9.5 Blade with 9 equidistant outer sections. . . . . . . . . . . . . . . . . . . 74
9.6 Load reduction of a single ap with a length of 1.5m at dierent radial
positions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
9.7 Load reduction of single aps with dierent lengths and dierent radial
positions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
9.8 Out-of-plane bending moment for baseline simulation and the single
13.5m aps with maximum load reduction. . . . . . . . . . . . . . . . . 78
9.9 Inuence of dierent actuator speeds on the load reduction for two ap
congurations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
List of Figures XI
9.10 Inuence of dierent actuator ranges and speeds on the load reduction
for a 9m ap conguration. . . . . . . . . . . . . . . . . . . . . . . . . 80
9.11 Flap angle for 9m ap over simulation time. . . . . . . . . . . . . . . . 81
9.12 Load reduction over controller time delay. . . . . . . . . . . . . . . . . 82
9.13 Local blade element force DFN at AE# 3 for the baseline, the single
PID controlled 13.5m ap and the multiple individually optimization
loop controlled aps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
List of Tables
8.1 Turbine parameters used for the simulation. . . . . . . . . . . . . . . . 62
8.2 Blade structural parameters used for the simulation. . . . . . . . . . . . 63
9.1 Radial position of the 9 possible active elements. . . . . . . . . . . . . . 75
9.2 Load reduction for 1.5m ap at dierent radial positions. . . . . . . . . 75
9.3 Load reduction for dierent ap lengths and dierent radial positions. . 77
9.4 Load reduction for individual sections using the optimization loop. . . . 83
A.1 Blade geometric parameters in NREL format used for the simulation. . 98
Nomenclature
Greek symbols
α . . . . . . . . . . . . . . . angle of attack
β . . . . . . . . . . . . . . . local blade twist
χ . . . . . . . . . . . . . . . wake skew angle
γ . . . . . . . . . . . . . . . yaw angle
Ω . . . . . . . . . . . . . . blade angular velocity
φ . . . . . . . . . . . . . . . inow angle
ψ . . . . . . . . . . . . . . . blade azimuthal angle
ρ . . . . . . . . . . . . . . . air density
θ . . . . . . . . . . . . . . . local pitch angle
θp . . . . . . . . . . . . . . blade pitch angle
Roman symbols
a . . . . . . . . . . . . . . . axial induction factor
a′ . . . . . . . . . . . . . . tangential induction factor
B . . . . . . . . . . . . . . number of blades
c . . . . . . . . . . . . . . . chord length
cd . . . . . . . . . . . . . . drag coecient
cl . . . . . . . . . . . . . . lift coecient
cm . . . . . . . . . . . . . . moment coecient
Ct . . . . . . . . . . . . . . rotor thrust coecient
D . . . . . . . . . . . . . . blade element drag force
DFN . . . . . . . . . . blade element normal force
DFT . . . . . . . . . . . blade element tangential force
F . . . . . . . . . . . . . . force, Prandtl's tip-loss factor
L . . . . . . . . . . . . . . . blade element lift force
pn . . . . . . . . . . . . . . load normal to rotor plane
pt . . . . . . . . . . . . . . load tangential to rotor plane
Nomenclature XIV
PMA . . . . . . . . . . blade element pitching moment
Q . . . . . . . . . . . . . . rotor torque
r . . . . . . . . . . . . . . . local blade element radius
T . . . . . . . . . . . . . . rotor thrust
U∞ . . . . . . . . . . . . . inow velocity
ve,ip . . . . . . . . . . . . out-of-plane velocity due to structural deection
ve,op . . . . . . . . . . . . in-plane velocity due to structural deection
W . . . . . . . . . . . . . . incident velocity
AE . . . . . . . . . . . . . Active Element
AFC . . . . . . . . . . . Active Flow Control
AFMB . . . . . . . . . axial ap bending moment (out-of-plane)
BEM . . . . . . . . . . . blade element momentum theory
CFD . . . . . . . . . . . Computational Fluid Dynamics
GDW . . . . . . . . . . generalized dynamic wake model
GUI . . . . . . . . . . . . Graphical User Interface
HAWT . . . . . . . . . Horizontal Axis Wind Turbine
ISTA . . . . . . . . . . . Institute of Fluid Dynamics and Technical Acoustics
NREL . . . . . . . . . . National Renewable Energy Laboratory
NWTC . . . . . . . . . National Wind Technology Center
1 Introduction
The use of wind energy is continuously growing and in order to increase it's market
competitiveness, wind turbines must become even more cost eective. The key param-
eter therefore is e/kWh as in any other energy system. The cost of electricity generated
by modern wind turbines ranges from approximately 0.05 - 0.07 e/kWh at sites with very
good wind speeds to 0.09 - 0.11 e/kWh at sites with low wind speeds [14]. To cut costs,
the expenses for production, operation and maintenance have to be reduced while at
the same time the performance of wind turbines has to increase. One possibility to
achieve this is to build larger wind turbines which can extract more wind energy from
a bigger swept area. However, the continuous increase of the rotor diameter leads to
structural problems due to the enormous size of the wind turbine blades. As the wind
uctuations over the swept area get higher, the loads on the blades due to wind shear
layer, tower shadow, yaw misalignment and atmospheric turbulence increase. In ad-
dition to that, the higher blade mass introduces higher cyclic loads, which also have
negative inuence on the lifetime of the blades.
One attempt to overcome the inherent limitations of upscaling and to reduce the
structural load on the wind turbine, is the introduction of an individual pitch control.
Such a system has two major drawbacks. Firstly, the pitch actuator is relatively slow
and cannot cope with high-frequency uctuations; secondly the actuation takes place
at the innermost part of the blade, whereas the highest load contribution comes from
the outermost part. Consequently, the elasticity of the blade waters down the control
circuit.
Another approach is the use of active ow control devices (AFC) on wind turbine
blades. Actuators for active ow control can be integrated directly in the blade, with
sensors, not only behind the rotor as current anemometers and ow vanes but lo-
cally near the actuator. Measuring and controlling the loads directly where they occur,
1 Introduction 2
plus having multiple, smaller and faster actuators makes it possible to further reduce
load uctuations.
Recently a signicant amount of research has been carried out on this topic. Al-
though the main focus seems to lie on the investigation of trailing edge devices, a wide
range of other possibilities might be of interest too, such as:
• Active mini-aps
• Flexible leading-edge aps
• Inatable stall ribs
• Inclined and vertical spoilers
In order to investigate the behavior of dierent AFC-solutions a software is developed
to determine the potential of these concepts. The software allows easy denition of
a wind turbine rotor blade with several so-called active elements (AE). Figure 1.1
shows an example of how such a rotor blade, equipped with several active elements,
may look like. On the outer region, there are four trailing edge ap devices integrated.
The active elements are characterized by a variable aerodynamic performance which
is expressed by multiple lift, drag and moment polars. The response of the blade and
the wind turbine can be investigated using an unsteady aeroelastic simulation. A bind-
ing to the aerodynamic analysis routine AeroDyn [38] and the structural dynamic code
YawDyn [28] is implemented. These simulation codes are provided by the National
Wind Technology Center (NWTC). The aerodynamic package features unsteady simu-
lation with dynamic stall eect modeling as well as and full turbulent wind eld input.
The structural dynamic model is simple, yet useful for analyzing preliminary designs
and assessing aerodynamic responses, as the simulation time is very low and a lot of
parameter investigations can be carried out.
Next to the core routines, a major driving design parameter of the software is user
friendliness. The user is able to perform simulations without manual script handling or
code related operations. To achieve an easy-to-use interface, a GUI was implemented
which automatically handles all the communication between the modules and provides
dynamic data visualization. The software itself is based on the open-source application
1 Introduction 3
Figure 1.1: Concept of a segmented wind turbine rotor blade with active elementsin form of trailing edge aps [52].
QBlade [35] and as it deals with the investigation of the aforementioned user-dened
active elements, the name speaks for itself: QBladeAE.
After describing the theory of the models used by the software in Part I, the second
Part II of this work presents the working principle of the software. In Part III a param-
eter study is performed on a normal blade conguration, to investigate the behavior
of the physical models. Finally, QBladeAE is being used to simulate an exemplary
blade with exible aps in dierent congurations. The results are compared with
the baseline conguration to point out the advantages of the AFC-solution for load
reduction on wind turbines.
The focus of this project lies on the development of the software and not on the
investigation and comparison of dierent AFC-solutions.
2 Aeroelastic model
As mentioned in the introduction QBladeAE works in conjunction with AeroDyn. This
set of FORTRAN subroutines contains a pure aerodynamic wind turbine model de-
scription and is provided by NREL. AeroDyn is no stand alone application and needs
to be coupled with a structural program, which provides information about the dy-
namic structural deections of the wind turbine and the elastic blades, as well as the
operating conditions (e.g. rotational speed and blade pitch angles). As the structural
deections induce changes in the aerodynamic forces, the computation gets fully aeroe-
lastic. Currently AeroDyn works together with three structural programs, which dier
in their level of complexity. These are YawDyn (which is used by QBladeAE), FAST
and ADAMS (4).
The structural program controls the whole turbine simulation and calls the AeroDyn
subroutines during runtime (once for each time step, blade and blade element), in
order to obtain the aerodynamic forces on the blades. The blade is split into several
blade elements and for each element the lift and drag forces as well as the pitching
moment is determined. The computation is broken down into a two-dimensional local
blade element formulation. All element velocities are accumulated and expressed in the
local blade element coordinate system. Finally, a resulting incident velocity W with
a resulting inow angle of attack α is determined. Using airfoil polar tables for the
lift coecient cl(α), drag coecient cd(α) and moment coecient cm(α), the resulting
element forces can be calculated. After integrating the forces over each blade, the
structural model updates the dynamic deections of the blade and the wind turbine
structure in the next time step. Consequently, the local element velocities change,
which in turn aects the blade element aerodynamics again. Figure 2.1 shows the
dierent portions of the velocities and inow angles seen by a single blade element.
2 Aeroelastic model 6
rotation plane
chord line
Ωr(1+a')Ve,ip
Ve,op
U∞(1-a)
ϕ
α
θ = θp + β
W
Figure 2.1: Local blade element velocities and inow angles.
The incident velocity W is given as
W = U∞(1− a) + Ωr(1 + a′) + ve,op + ve,ip (2.1)
where U∞ is the (unsteady) inow velocity, Ωr the local element circumferential speed,
a and a′ the axial and tangential induction factors and ve,op and ve,ip the in-plane and
out-of-plane velocities due to the structural deection. The total inow angle φ is a
combination of the angle of attack α and the local element pitch angle θ, which is again
consists of the blade pitch angle θp and the local blade twist β.
As described in 3.2, the ow around the airfoil produces the aerodynamic lift force
L perpendicular to the incident velocity W and the aerodynamic drag force D in ow
direction (Figure 2.2). The resulting force can then be split into a portion perpendicular
to the rotor plane, the normal load pn = L cos(φ) +D sin(φ) and a smaller tangential
portion pt = L sin(φ) +D cos(φ) which contributes to the rotation of the rotor.
2 Aeroelastic model 7
rotation plane
chord line
W
D
Lϕ
pn
pt
Figure 2.2: Local blade element forces.
To calculate the aerodynamic and structural forces on the wind turbine, is the task
of the aeroelastic model of the software. The methods, working principles and assump-
tions made within the model are discussed in the following chapters. The aeroelastic
problem is split into the aerodynamic model and the structural dynamics model.
3 Aerodynamics
The aerodynamic model used by QBladeAE contains representations of several dierent
aerodynamic eects on a wind turbine (HAWT) which will be described below. This
includes wake modeling, airfoil aerodynamics (static and dynamic), yawed inow, tower
shadow, shear layer and atmospheric turbulence eects. Figure 3.1 gives a general
overview of the periodic and aperiodic aerodynamic sources on a wind turbine.
Flowfield Structure
Mostly periodic Mostly Aperiodic
Wind
Speed
Inflow Yaw Tower
Shadow
Wind
turbulence
Wake
dynamics
Blade/
wake
interactions
Figure 3.1: Summary of the various aerodynamic sources that contribute to theairloads on a wind turbine [32].
As can be seen, wind turbines operate under extreme unsteady aerodynamic condi-
tions, which are hard to dene, to measure and to predict with mathematical models
[32]. The approaches to describe these eects with the models of AeroDyn are described
below.
3.1 Wake modeling
Wind turbines extract kinetic energy from the wind. The air approaches the turbine
and is slowed down. As it passes the rotor, there is a step drop in pressure and the
velocity further decreases, as the pressure has to reach the atmospheric level again.
The dierence in the air ow velocity before and after the turbine accounts for the
extracted energy. To calculate the aerodynamic ow across the rotor, dierent models
3 Aerodynamics 9
are available. They reach from simple blade element/momentum theory over more
advanced engineering models to complicated full scale CFD simulations. The former is
the oldest and most common approach to model the aerodynamics of wind turbines and
to calculate the velocity decit in the rotor plane. The latter have high computational
costs, but provide a more realistic physical model, as they solve the Navier-Stokes
equations. In between these two ends, there are various other models which are mostly
adapted from the helicopter industry. These engineering models usually combine the
blade element theory with either a dynamic inow or a vortex wake model [32]. A major
dierence between the engineering models and the BEM theory is the the modeling of
unsteady wake eects. The time dependent changes in the inow and the blade loading
can be treated in dierent ways, by making one of the following assumptions:
Frozen wake One assumption could be that small changes in the inow have no inu-
ence on the induced velocities. The wake is only dependent on the average wind
speed over a (short) period of time. This means the unsteady wind component
passes the rotor unattenuated [8].
Equilibrium wake On the other side it could be assumed, that the wake instanta-
neously reacts on changes in the aerodynamic loading. The induced velocities
change as the inow changes and therefore the wake is always in equilibrium. For
the simulations, this means that the induction factors have to be re-calculated
for every blade element and time step. Most blade element momentum theories
make use of this assumption.
Dynamic wake In reality neither of these assumptions is correct and the truth lies
somewhere in between. Changes in the inow change the vorticity that is trailed
into the rotor wake and the full eect of these changes takes a nite time to
change the induced ow eld [5]. As indicated above, the most common method
to model dynamic inow eects is a combination of the blade element theory and
a dynamic inow model.
How important the consideration of dynamic inow eects is especially for fast
pitching transients and yawed conditions was investigated within the the European
Union JOULE 1 and JOULE 2 programs [46], where dierent models were compared
to each other, with the model of Pitt and Peters [44] probably being the most common
one. For the sake of completeness it shall be mentioned, that there are as well models,
3 Aerodynamics 10
which implement a dynamic wake formulation in the blade element momentum theory,
like outlined in [18].
AeroDyn provides two ways of calculating the induced velocities in the rotor plane:
the classic blade element momentum (BEM) theory with the equilibrium wake as-
sumption and the general dynamic wake model (GDW). Both are be described in the
following paragraphs.
3.1.1 Blade element momentum theory (BEM)
As mentioned above, the blade element momentum theory is one of the oldest methods
for wind turbine wake modeling. It is a combination of the blade element theory, in
which a blade is split in several independent sections, and a momentum theory, which
attributes the loss of momentum in the rotor plane to the aerodynamic eects of the
ow passing the blades.
Blade element theory
In the blade element theory, a blade is regarded as a number of independent blade
sections or elements. According to the airfoil theory (3.2), the aerodynamic forces
which act on the blade element can be calculated by means of two-dimensional airfoil
characteristics. With the information of the absolute value and the direction of the
incident velocityW at the blade element, the lift and drag forces as well as the pitching
moment can be determined by using airfoil tables. It is assumed that each element
cuts out an annular ring section of the rotor disc (Figure 3.2), and that the overall
rotor performance is the integration over the single annular rotor sections.
The lift and drag forces on the blade element with the chord length c (Figure 2.2)
are given as
dL =ρ
2clcdrW
2 (3.1)
dD =ρ
2cdcdrW
2 (3.2)
3 Aerodynamics 11
dr
R
r
Figure 3.2: Rotor of a three-bladed wind turbine with the rotor radius R [18].
with the trigonometric relations (Figure 2.1) for the inow angle φ = θ + α
sin(φ) =U∞(1− a)
W(3.3)
cos(φ) =Ωr(1 + a′)
W(3.4)
Blade element and momentum theory (BEM)
The combination of the blade element theory with a momentum theory allows the cal-
culation of the induction factors a and a′, which are necessary to determine the incident
velocity W in the rotor plane. It is assumed that each blade element is responsible for
the change of momentum of the air, which passes through the annulus swept by the
element [8]. A detailed derivation can be found in many textbooks, such as [18] and
only a short introduction is given here.
To determine the induction factors, the aerodynamic forces which contribute to the
thrust T and the torque Q are set equal with the change of momentum in the annulus
3 Aerodynamics 12
section. The aerodynamic forces of B blades in axial and tangential direction result in
dT = Bρ
2W 2cdr(cl cosφ+ cd sinφ) (3.5)
dQ = Bρ
2W 2cdr(cl sinφ+ cd cosφ)r (3.6)
The momentum theory gives
dT = 4πρU2∞a(1− a)rdr (3.7)
dQ = 4πρU∞(Ωr)a′(1− a)r2dr (3.8)
These equations can now be solved iteratively using two-dimensional airfoil data. Note
that AeroDyn takes the additional velocities ve,op and ve,ip into account for determining
the absolute value and inow angle φ of the incident velocity W but does not consider
them in the momentum theory, which might not be the appropriate physical model for
the element-wake coupling [38].
Despite it's simplicity, the BEM theory provides relatively accurate results. There
are other aerodynamic eects on a real turbine, which can not be modeled with the
BEM method directly, because of the assumptions made in the theory. These are eects
due to heavy loaded rotors with high induction factors, blade tip and hub losses due to
a limited number of blades and skewed inow which is not perpendicular to the rotor
plane. AeroDyn includes several corrections to account for these eects:
Tip loss model The fact that vortices are being shed from the blade tip causes high
axial induction factors, which leads to lower inow velocities at the rotor. This
causes to smaller inow angles φ and most of the aerodynamic lift contributes to
the thrust. Less torque means less power and therefore the losses near the blade
tips are higher.
AeroDyn features two models to calculate the blade tip losses. First of all, the
classic model developed by Prandtl. An additional correction term F is added to
the momentum Equations 3.7 and 3.8
F =2
πcos−1 e−f (3.9)
(3.10)
3 Aerodynamics 13
where
f =B
2
(R− r
r sinφ
)(3.11)
The other model used in AeroDyn slightly modies the Prandtl correction factor
F , using an empirical relationship for the tip losses on base of the Navier-Stokes
solutions [38]:
Fnew =F 0.85Prantl + 0.5
2, for 0.7 ≤ r/R ≤ 1 (3.12)
Fnew = 1−( rR
) 1− FPrantl(r/R=0.7)
0.7, for 0.7 ≤ r/R ≤ 1 (3.13)
Hub loss model The eect of a vortex in the hub region is described using a nearly
identical implementation of the tip-loss model. Equation 3.11 is replaced with
f =B
2
(r −Rhub
Rhub sinφ
)(3.14)
Turbulent wake state The standard momentum equation used in the BEM theory
gives negative thrust values for induction factors greater than 0.5. However,
this is not what happens in reality. As the wake becomes turbulent for heavily
loaded rotors, air (thus momentum) is transported from the outer ow region
into the wake. To account for this eect, the empirical correction of Glauert is
implemented in AeroDyn and is slightly modied, to avoid numerical instability:
CT =8
9+
(4F − 40
9
)a+
(50
9− 4F
)a2 (3.15)
Skewed wake To be able to describe the eects of yaw misalignment, AeroDyn pro-
vides a skewed wake correction. The model is based on the work of Glauert
(1926) and was extend by Pitt and Peters (1981). For steady inow conditions,
the local element induction factor askew is given with
askew = a
[1 +
15π
32
r
Rtan
χ
2cosψ
](3.16)
where ψ is the azimuthal angle that is zero at the most downwind position of the
3 Aerodynamics 14
rotor and χ being the wake skew angle, which can be approximated using the
yaw angle γ as follows [38]:
χ ≈ (0.6a+ 1)γ (3.17)
Despite the original assumption made by Glauert, AeroDyn applies the induction
factor askew to all local elements [38].
3.1.2 Generalized dynamic wake model (GDW)
The generalized dynamic wake model of AeroDyn is also known as acceleration potential
method and is based on the work of Peters and He (1989), which again is based on
the aforementioned model of Pit and Peters [44]. The main advantage over the BEM
method is the inherent inclusion of dynamic wake eects, tip losses and skewed wake
aerodynamics [38]. The equations describe the distribution of inow and are written
in the form of dierential equations, which can be solved non-iteratively. The GDW
model has several drawbacks as well. Theses are:
• Instabilities at low wind speeds when the turbulent wake state is approached.
AeroDyn uses the BEM method for wind speeds below 8m/s.
• No accounting for wake rotation. AeroDyn uses the BEM method as well to
calculate the tangential induction factor.
• Flat disk assumption makes the eect of large aeroelastic deections inaccurate.
The method itself is based on the unsteady and inviscid Euler equations. Assum-
ing the induced velocities are small against the wind velocity U∞ the conservation of
momentum can be written as
∂u
∂t+ U∞
∂u
∂x= −1
ρ
∂p
∂x(3.18)
∂v
∂t+ U∞
∂v
∂x= −1
ρ
∂p
∂y(3.19)
∂w
∂t+ U∞
∂w
∂x= −1
ρ
∂p
∂z(3.20)
3 Aerodynamics 15
and for continuity of the ow
∂u
∂x+∂v
∂y+∂z
∂w= 0 (3.21)
and the Laplace equation for the pressure distribution
∇2p = 0 (3.22)
The boundary conditions are the aerodynamic forces on the loaded blade, the pres-
sure returns to ambient pressure far behind the rotor and the equality of discontinuous
pressure and rotor thrust. The pressure eld is then split into a term for the spatial
variation and a term for the unsteadiness to split the unsteady Euler equations accord-
ingly.
A pressure distribution, which gives a discontinuous pressure drop across the rotor
and satises the Laplace equation was developed by Kinner (1937). A more detailed
description of the method is found in [38] and [8].
According to [10] the GDW method for calculating yawed and dynamic inow is
surprisingly good for its computational simplicity. However, it contains many simpli-
fying assumptions and it is proposed to implement a free vortex wake method for more
accurate results instead. On the other hand, the disadvantage of a free vortex model
is the long computation time. As noted in [47], a 10 minute time simulation with the
advanced Alcyone free wake model lasted 5 days.
3.2 Airfoil aerodynamics
3.2.1 2D static airfoil characteristics
Most wind turbine models including the ones described above make use of two-
dimensional static airfoil tables. The assumption that the ow around the blade at a
given radial position is two-dimensional, as indicated in Figure 3.3, is not always valid
especially in the blade root and tip region. On the other hand, The advantage of
having static airfoil look-up tables for the aerodynamic forces as a function of the angle
of attack α is very useful for the aerodynamic simulation. This approach in contrary
to an on-the-y calculation allows to import airfoil characteristics from wind tunnel
3 Aerodynamics 16
measurements or complex numerical computations. It is obvious that the key to an
accurate simulation lies in the careful provision of valid airfoil properties. Unfortunately
this is a hard task, as wind tunnel measurements for very high Reynolds numbers
are very costly and valid CFD simulations very time consuming and computationally
intensive.
Lift
U∞
Figure 3.3: 2D arfoil characteristics of a blade element.
Once the aerodynamic lift- drag and moment coecients cl, cd and cm are known,
the resulting forces for lift L, drag D and pitching moment M can then be calculated
using the denition:
cl(α) =L
ρ2U2∞c
(3.23)
cd(α) =D
ρ2U2∞c
(3.24)
cm(α) =M
ρ2U2∞c
2(3.25)
As stated in [38] and [49] the largest source of error in load and performance simu-
lations are errors in the airfoil data tables. Figure 3.4 - 3.6 show exemplary how the
characteristics for the same airfoil can varies under dierent conditions. The results are
computed with XFOIL, a 2D panel method which includes an estimation for viscous
ow [12]. The graphs show a calculation for two Reynolds numbers and an additional
calculation for xed transition near the leading edge.
3 Aerodynamics 17
-1
-0.5
0
0.5
1
1.5
2
-10 -5 0 5 10 15 20
c l
α
Re = 5e6, Ma = 0.1, xtrf = 1.0Re = 2e6, Ma = 0.1, xtrf = 1.0Re = 2e6, Ma = 0.1, xtrf = 0.1
Figure 3.4: Lift coecient cl(α) at dierent Reynolds numbers and xed/free tran-sition for the DU 91-W2-250 airfoil.
As the operational conditions of the wind turbine are changing during the simula-
tion, the airfoil performance changes as well. With increasing wind speed or changing
rotational speed (in spanwise direction), the Reynolds number varies. That makes it
hard to cover the whole range of operation in a simulation with only one set of polars.
AeroDyn provides the possibility to dene multiple tables for one airfoil. The user can
specify dierent tables for dierent Reynolds numbers. These tables are dynamically
accessible during simulation. As will be seen later, this functionality can be used to
dene the dierent aerodynamic characteristics of the aforementioned active elements
too.
3 Aerodynamics 18
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
-10 -5 0 5 10 15 20
c d
α
Re = 5e6, Ma = 0.1, xtrf = 1.0Re = 2e6, Ma = 0.1, xtrf = 1.0Re = 2e6, Ma = 0.1, xtrf = 0.1
Figure 3.5: Drag coecient cd(α) at dierent Reynolds numbers and xed/free tran-sition for the DU 91-W2-250 airfoil.
-0.15
-0.14
-0.13
-0.12
-0.11
-0.1
-0.09
-0.08
-0.07
-0.06
-0.05
-10 -5 0 5 10 15 20
c m
α
Re = 5e6, Ma = 0.1, xtrf = 1.0Re = 2e6, Ma = 0.1, xtrf = 1.0Re = 2e6, Ma = 0.1, xtrf = 0.1
Figure 3.6: Moment coecient cm(α) at dierent Reynolds numbers and xed/freetransition for the DU 91-W2-250 airfoil.
3 Aerodynamics 19
3.2.2 Polar extrapolation
Unlike airplane wings, wind turbine blades experience stalled operation. The rotational
speed of the blade gets higher towards the blade tip, but the average wind inow
velocity u∞,mean remains constant. This results in higher ow angles φ in the root
region, as can be seen in Figure 3.7. To compensate for this and to keep the angle of
attack α constant over the span width, the blades are structurally twisted inwards
more than outwards. As mentioned above the total inow angle φ is given as (assuming
a blade pitch angle of θp = 0):
φ = θp + β + α = β + α (3.26)
To keep the angle of attack α constant, the twist β must increase when the total inow
angle φ increases.
The conventional manufacturing process however, allows only a limited blade twist.
The root region is likely to operate under stalled conditions. Figure 3.8 shows the
angle of attack α for an inner and an outer section of a 40m blade with limited twist
of βroot = 13 in the root region. The turbine is operating at a rotational speed of
n = 16rpm at u∞,mean = 13m/s. The average value of the angle of attack is αmean = 20.
It's obvious that the airfoil tables need to be extended to a wider range of angles of
attack. But stall phenomena are viscous eects and it is anything but trivial to nd
valid numerical or experimental ways to determine the behavior of airfoil characteristics
beyond stall. Methods based on the potential ow theory (like used by Xfoil) are only
able to include viscous eects by semi-empirical models. Wind tunnels measurements
are also complicated, due to the high blockage in the measurement section for high
angles of attacks.
One way to overcome this problem is to use airfoil characteristics for normal operation
and extrapolate them by using the at plate theory. This approach points out the
similarity between a at plate and an airfoil at high angles of attack. This method is
refereed to as the Viterna method [55]. As wind turbine airfoils used for the root region
are relatively thick, the model can be further adapted. For the 360-extrapolation in
QBladeAE, the empirical method described in [37] is used. A similar approach is
described in [50]. Figure 3.9 and 3.10 show the extrapolation of the initial values for cl
3 Aerodynamics 20
Ωrout
Ωrmid
Ωrin
vtot
Ωrout
Ωrmid
Ωrin
u
u
u
ϕ
Figure 3.7: Wind triangular for dierent radial positions [16].
5
10
15
20
25
30
35
0 10 20 30 40 50 60
α [d
eg]
time [s]
ri = 6.8mro = 36.8m
Figure 3.8: Exemplary time series of α for two radial positions at rin = 6.8m androut = 36.8m.
3 Aerodynamics 21
-1.5
-1
-0.5
0
0.5
1
1.5
2
-180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180
c l
α
2*sin(x)*cos(x)Re = 5e6
Figure 3.9: Extrapolated cl for the DU 91-W2-250 airfoil and cl according to theat plate theory.
and cd seen in 3.4 and 3.5. In addition to that, the lift and drag coecients according
to the at plate theory are shown as well:
cl,fp(α) = cd,90 sinα cosα (3.27)
cd,fp(α) = cd,90 sin2 α (3.28)
As described in [35] the user can inuence the extrapolation via several control vari-
ables, to adapt the method to dierent airfoil characteristics. It shall be noted, that
the 360-extrapolation has no eect on the original two-dimensional airfoil data and
the original polar needs to cover an angles of attack range right up to the stall point.
Errors made in the extrapolation process mainly eect the root region of the blade,
whose contribution to the energy yield and the structural load is naturally smaller.
Nevertheless, the accuracy of valid airfoil data has to be pointed out again.
Next to the physical operation under high angles of attacks, the 360-extrapolation
of airfoil data is as well needed for the classical BEM-method (3.1.1). This has numer-
ical reasons and is necessary for a successful convergence of the iterations in the BEM
method.
3 Aerodynamics 22
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
-180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180
c d
α
1.98*sin(x)*sin(x)Re = 5e6
Figure 3.10: Extrapolated cd for the DU 91-W2-250 airfoil and cd according to theat plate theory.
3.2.3 Unsteady blade element aerodynamics
Until now, only the steady nature of the airfoil performance was described. The airfoil
polars are so far not more than simple look up tables. In reality, the blades operate
under highly unsteady conditions. Firstly, there is an angle of attack change α over
time, due to blade pitching or yaw misalignment. Secondly, there are in-plan inow
velocity changes U∞(t) from wind gusts as well as out-of-plane velocity changes due to
wake interactions. According to [32] it is important to distinguish between these two
inuences on the airloads and to treat them separately. However, AeroDyn does treat
every change in angle of attack equally, whether they arise from blade pitch motions,
changes in the relative wind velocity or blade ap or lag motions.
Attached ow conditions
The unsteadiness at blade element level is often put on a level with dynamic stall. It
has to be noted, that there are as well unsteady eects in the attached ow region with
low angles of attack. These eects become noticeable in moderate amplitude and phase
3 Aerodynamics 23
variations compared to the steady eects [32] provided that the reduced frequencies
are small. One way to describe the unsteady eects in the linear lift region was given
by Theodorsen and is known as the Theodorsen's theory [53]. The theory is based
on the incompressible and inviscid ow theory for thin airfoils. There are extensions
which have been developed, including compressible ow and wing rotation eects. As
AeroDyn does not include the computation of unsteady airfoil aerodynamics under
attached ow conditions, only a reference is made to [32], where a general insight in
the unsteady aerodynamics of HAWTs is given.
Dynamic stall
Leaving the linear lift region of an airfoil, stall eects occur. Dynamic stall is describing
the phenomenon of delayed stall occurrence on airfoils under unsteady conditions
either a time varying in the inow or the angle of attack. An airfoil which oscillates or
pitches through the static stall region, experiences delayed stall onset, but considerably
stronger and longer-lasting stall eects compared to the static stall development. A
cycle of a pitching airfoil, which experiences a dynamic stall hysterisis is shown in
Figure 3.11 and can be described as followed: If the static stall angle is exceeded, there
is no immediate change in the viscous or inviscid ow around the airfoil, due to a time
lag until stall takes place. After the rst appearance of ow reversal in the boundary
layer near the trailing edge, the reversal ow moves further upwards to the leading
edge. This is when rst large eddies appear and a vortex is formed at the leading edge.
This vortex rolls up, gains in strength and is shed downstream, producing an increasing
lift slope. At the same time the center of pressure moves along with the vortex, causing
a large nose down pitch moment. After the vortex has passed the trailing edge, the lift
drops rapidly and full stall takes place. When the angle of attack decreases and falls
below the static stall value again, the ow re-attaches from the leading to the trailing
edge again [43].
It is obvious that the unsteady stall eects on airfoils cannot generally be included
in two-dimensional static airfoil tables. To cover theses eects, an dynamic stall model
is implemented in AeroDyn which shall be briey described according to [43]. The
model is based on the work of Beddoes and Leishman [31]. It is a semi-empirical model
which adapts the attached ow indicial1response of an airfoil to the position where ow
1By denition, an indicial function is the response to a disturbance that is applied instantaneously
3 Aerodynamics 24
(a) STATIC STALL ANGLE EXCEEDED
(B) FIRST APPEARANCE OF FLOW
REVERSAL ON SURFACE
(c) LARGE EDDIES APPEAR IN
BOUNDARY LAYER
(d) FLOW REVERSAL SPREADS OVER
MUCH OF AIRFOIL CHORD
(e) VORTEX FORMS NEAR
LEADING EDGE
(f) LFIT SLOPE INCREASES
(g) MOMENT STALL OCCURS
(h) LIFT STALL BEGINS
(i) MAXIMUM NEGATIVE MOMENT
(j) FULL STALL
(k) BOUNDRAY LAYER REATTACHES
FRONT TO REAR
(l) RETURN TO UNSTALLED VAUES
Figure 3.11: Dynamic stall events on a NACA 0012 airfoil (reprinted from [9] and[32]).
3 Aerodynamics 25
separation actually takes place. With an indicial response function φα, the change in
the normal lift coecient ∆cn for a angle of attack change ∆α for the attached ow
region is given as
∆cCn = cnαφCα∆α (3.29)
∆cIn =4
MaφIα∆α (3.30)
with cnα being the slope of the normal coecient in the linear region, Ma the Mach
number. The lift coecient is additionally separated in one component for the circu-
latory part cCn and one for the non-circulatory part cIn. The response of the tangential
force coecient cc in chord wise direction is derived from the circulatory part of cn. The
attached ow indicial response is then adapted to the separation point f of the suction
side of the airfoil. With the static airfoil data, the separation point f is determined by
the relation
cn = cnα(α− α0)(1 +
√f
2)2 (3.31)
cc = cnα(α− α0) tan(α)√f (3.32)
where α0 denotes the angle of attack for zero lift. As these equations are derived from
an inviscid formulation, f is referred to as the static eective separation point, and
might not be the exact point of reversal ow appearance. A empirical time lag is fur-
ther applied to the movement of the eective separation point to account for the time
lag of the real separation point under unsteady conditions. In a last step, the vortex
shedding across the upper surface of the airfoil is modeled, as soon as a critical leading
edge pressure parameter indicates leading edge separation. This results in the typical
lift increase until the airloads return to their static values. The relation between cn
and cl can be found in Figure 2.2.
The dynamic stall model described above, is modied slightly in AeroDyn. The main
dierences are:
• Extension to very high/low angles of attack
• The eective separation point f(α) is not curve tted by an exponential function
at time zero and held constant thereafter; that is a disturbance given by a step function [32].
3 Aerodynamics 26
but treated in a look up table with linear interpolation in between
• Two separate point tables are used, one for cn and cc
The advantage of the model is, that it uses very few empirical coecients, mostly
derived from the static airfoil tables. The airfoil input les for AeroDyn contain the
necessary values for the dynamic stall models. QBladeAE automatically calculates and
exports the values, (7.2.6) which are:
• Angle of attack for zero lift α0
• cn slope for zero lift
• cn at stall value for positive angle of attack
• cn at stall value for negative angle of attack
• Minimum cd value
• Angle of attack for minimum cd value αcd,min
The model provides a fairly accurate way to predict the unsteady eects of dynamic
stall but it has to be noted, that none of the available models are developed to full
extend and future investigations have to be made on this topic [32].
Furthermore the dynamic model used for simulating active elements in QBladeAE is
not changed according to the specic active elements used and has to be applied with
caution. Unfortunately it is not possible to derive a generally adapted model for each
dierent kind of active ow control actuator for example leading edge or trailing edge
ap. A modied dynamic stall model for trailing edge aps can be found in [1].
3.2.4 Stall delay and 3D eects
Measurements show, that conventional HAWT simulations can under-predict the power
output compared to measured data. Beside the under estimation of delayed stall due
to dynamic stall eects, another reason for the under-prediction are three-dimensional
ow eects which are caused by centrifugal and Coriolis forces. These forces can have
a positive eect on the pressure gradient on the suction side, so that stall is delayed
[8]. Numerical investigations support these results [3] and point out their importance.
3 Aerodynamics 27
Another eect responsible for stall delay are the incident ow velocities which result
in a realtive wing sweep [32].
The eects mentioned above, are still undergoing research and not form of three-
dimensional airfoil correction is included in AeroDyn. It has to be noted as well, that
the BEM theory (3.1.1) can not handle three-dimensional eects, as the blade sections
are considered to be independent of each other by denition. If there is any desire
for implementation of three-dimensional eects, the static airfoil tables have to be
corrected and modied manually.
3.3 Tower shadow
The inuence of the tower shadow can be modeled as a velocity change seen by the
blade. The inuence is manifested in a velocity decit normal to the rotor plane.
AeroDyn uses two models to simulate the tower eects. They both use a potential ow
around a cylinder as basis and superimpose either a dam model for upwind turbines
and an additional wake model for downwind turbines. Based on the tower diameter a
drag value cd,tower is used to calculate the dimensionless velocity eld according to
u = 1− (x+ 0.1)2 − y2
((x+ 0.1)2 + y2)2+cd,tower
2π
x+ 0.1
(x+ 0.1)2 + y2(3.33)
v = 2(x+ 0.1)y
((x+ 0.1)2 + y2)2+cd,tower
2π
y
(x+ 0.1)2 + y2(3.34)
(3.35)
where u and v are the components of the horizontal wind in the x and y direction. The
parameters x and y are the upwind and crosswind distances normalized by the tower
radius [38].
For the tower wake model of downwind turbines, which are of no interest within this
report, reference is made to [38] as well.
4 Structural dynamics
Especially the exible blades and the tower make a wind turbine is a highly dynamic
system. To model the behavior of the turbine, a dynamic structural model is necessary
for several reasons: Firstly to determine extreme loads for the certication process,
secondly for the time dependent load variations on the components for fatigue load
calculation, in the third place to calculate deections which inuence the aerodynamic
model and nally to analyze the stability of the design. Dierent kind of loads are
acting on the structure [42]:
Aerodynamic loads The aerodynamic loads are listed in Figure 3.1.
Gravitational loads The weight of the blades and the nacelle resulting in a force
pointing downwards. Blade mass imbalances cause additional periodic forces.
Inertial loads They include centrifugal and gyroscopic forces as well as acceleration
forces.
Operational loads Transient turbine operation loads, initialized by the control sys-
tem, such as starting up, pitching, breaking or yawing.
Currently lots of wind turbine analysis codes available. They all include an aerody-
namic model which is either a BEM method or other engineering models (3.1) and
a structural model, which describes the motions and deformations of the wind turbine.
Furthermore the simulation models include a representation of the exible drive train,
an electric model for the generator and interfaces for implementing control strategies
for the wind turbine operations. The latter are necessary for modeling the high dy-
namic stresses a wind turbine is suering from, during maneuvers like emergency stops
or power regulation. The main components of a state-of-the-art aeroelastic simulation
code for wind turbines is shown in Figure 4.1.
4 Structural dynamics 29
Blades
Tower
Foundation
Generator
Break
Drive Train
Control
Nacelle
Figure 4.1: Components of a HWAT structural model
4 Structural dynamics 30
The structural model of the turbine itself, can be generally described by Newton's
second law
Mx+ Cx+Kx = F (4.1)
with M being the mass matrix, C the damping matrix and K the stiness matrix. To
solve this set of equations there are dierent approaches used in wind turbine simulation
codes. There are mainly three types of models, which vary in their level of complexity:
Assumed mode shapes An (assumed) modal representation is used for the dynamic
modeling. The modal properties of the rotating blades and the non-rotating
tower are computed independently by using information like mass and stiness
distribution.
Multi Body System A multi body system describes the dynamic structure with only
a few rigid and eventually exible elements, which are coupled with joints. The
advantage of a MBS system is, that large displacements can be modeled.
Finite Element System The nite element method is used to approximately nd the
solutions to the partial dierential equations of the mechanical system. A large
but nite number of elements are used to mesh the structure, resulting in high
computational cost. The code is usually only used for layout design and stress
calculations but not for dynamic wind turbine models.
As mentioned above, there are several design codes available. Almost every major
research center has developed their own code. There are commercial products, like
GH Bladed with license costs of ¿30.000 [15] and free open-source codes like NREL
provides them. An overview of existing codes can be found in [36] and [42]. Some of
the most common design codes are:
GH BLADED from GL Garrad Hassan,
FLEX5 from the Technical University of Denmark,
HAWC2 from the Riso National Laboratory for Sustainable Energy DTU,
DUWECS from the Delft University of Technology,
ADAMS/WT from MSC Software in collaboration with the National Renewable Re-
search Laboratories
4 Structural dynamics 31
FAST-AD from the National Renewable Energy Laboratories
YawDyn from the National Renewable Energy Laboratories
Some of these codes were compared to each other in [39] and validated against mea-
surements in [47] and [48], showing sometimes big discrepancies between each other
and between simulation and wind tunnel experimental data.
As with the aerodynamic model, it is not intended to develop a new design code
within this project. Regarding the time constraint, a comparable level of complexity
could not have been reached. As mentioned above, the AeroDyn simulation routines
from NREL are used for the aerodynamic wind turbine calculation and for simplicity
the choice for the structural model to work with AeroDyn came down to YawDyn1.
Although it is a very simple model, it provides rst insight in the aeroelastic behavior of
wind turbines. As the focus of the project lies mainly in the preliminary comparison of
dierent AFC solutions rather than on the detailed investigation of a single approach
the provided model of YawDyn seems to be sucient, although it was mainly developed
for the investigation of yaw motions. In the following, the model used in YawDyn is
described.
4.1 YawDyn
YawDyn was developed in 1992 by the National Wind Technology Center (NWTC) at
NREL. It was preliminarily used to investigate the yaw dynamics of HAWTs. Next to
the aerodynamic models used in AeroDyn, it provides the structural response of the
wind turbine at xed rotational speed in a fully turbulent wind eld.
The following assumptions are made in the structural model of YawDyn:
• Only the yaw motion γ and the blade apping motion β are used in the devel-
opment of the equations of motions.
• The rotor can be either modeled as apping rotor with two or three blades, a
teetering rotor for two blades or completely rigid.
1Recently NREL stopped the support for YawDyn and does not recommend it any more. The useof FAST, which is certied from Germanischer Lloyd WindEnergie [7], is proposed instead.
4 Structural dynamics 32
Figure 4.2: The equivalent hinge-spring model for the blade ap degree of freedom[28].
• During simulation the system can operate at either a xed yaw angle, at a xed
yaw rate or with free yaw motion using parameters for yaw spring stiness, yaw
damper coecient and constant yaw friction moment.
• Upwind/Downwind rotor simulation with tilted rotor (τ) and precone blade angle.
• Blade pitch and lag motions are not considered, as they are not important to the
yaw response.
• The tower, the rotor shaft, the nacelle and the blades themselves are treated as
rigid bodies.
• The turbine can only be modeled at xed blade pitch and at a xed rotation rate
Ω. No controller input is implemented.
• The blades are described by a uniform mass distribution, their distance from the
hinge to the blade center of gravity, their mass moment of inertia about the hinge
axis and their torsional stiness of the blade root spring.
The hinge-spring model for the blade ap degree of freedom is shown in 4.2. The model
of the wind turbine is shown in Figure 4.3. For more information it is refered to [28]
and [17].
4 Structural dynamics 33
Figure 4.3: View of the HAWT dening selected terms and coordinate systems. Allangles are shown in their positive sense. The bold X,Y,Z axes are xedin space and are the coordinates in which the wind components aredened (VX, VY, VZ). Note that blade azimuth is zero when the bladeis at the 6 o'clock position [28].
5 Turbulent wind simulation
The wind inow seen by the wind turbine is everything else but uniform. The tur-
bulent atmospheric boundary layer leads to a vertical wind prole and high turbulent
structures in the ow eld. These unsteady turbulences seen by the wind turbine are
stochastic, but not completely random (e.g. white noise). There is a spatial and fre-
quency dependent correlation of the turbulence. Veers developed a model for turbulent
wind eld computation, which can be used as input for the aforementioned simulation
codes [54]. A disadvantage of this model is, that the time histories of the wind eld for
the three velocity components u, v and w are computed independently. In other words,
there is no correlation between them. The Mann model, which is based on lineralized
Navier-Stokes equations, takes this additional correlation into account [33].
QBladeAE provides two methods to generate turbulent wind eld input les for
the simulation. An internal windeld generator and a GUI for generating input les
for NREL's TurbSim [6], which works seamlessly together with the two other NREL
modules 1. Both modules are based on Veers' model. Figure 5.1 shows the hub height
wind speed and the inow velocity seen from the rotating blade at the blade tip position,
computed with TurbSim. The mean wind speed is 13m/s. Figure 5.2 shows the x-
component of a wind eld generated in QBladeAE. The eld has 20x20 points.
1Note, that the additional coherent structures in TurbSim are only compatible to AeroDyn version13, which is incompatible with YawDyn.
5 Turbulent wind simulation 35
9
10
11
12
13
14
15
16
17
0 10 20 30 40 50 60
v x [m
/s]
Time [s]
vhubvtip
Figure 5.1: Wind speed at hub height and inow velocity at blade tip includingrotation, wind shear and tower eect.
Figure 5.2: 3D wind eld from QBladeAE with 20x20 points.
6 Active Flow Control
The preliminary task of QBladeAE is to provide a possibility to investigate elements on
a wind turbine blade, which can actively inuence and control the ow eld: so called
active elements. Flow control provides a possibility to meet the growing problems on
large scale wind turbines due to their high blade mass. The blade mass increases with
the length of the blade by the power of three: mblade ∝ R3blade, whereat the power out-
put only increases with the power of two: P ∝ R2blade. To reduce the arising uctuating
loads on the blade and to develop mass optimized blades active (and passive) load
control becomes more and more interesting.
An obvious solution to meet the challenges of load control on wind turbines is to use
already existing systems. The pitch system on modern wind turbines was introduced
for power regulation, but it can used as well to alleviate the load uctuations either
in a cyclic or individual pitch motion. As can be seen in [26] or [24] these systems
have high potential, especially for the cyclic load compensation of the 1p frequency
and multiples of it. However, the pitch system has some inherent problems. Firstly
it is too slow to react on higher turbulent load uctuations. In addition to that, the
pitch acts always on the whole blade and can not cope for local disturbances on the
blade, like local wind gusts. Thirdly the actuator is located at the blade root, but the
highest potential for load reduction lies in the outer regions of the blade.
To meet these problems, other ways of active ow control (AFC) for load reduction
can be introduced. There is a variety of dierent solutions available. These include
• Trailing edge devices (Rigid Flaps, Split Flaps, Flexible Flaps)
• Leading edge devices (Slats, Flexible Leading Edge)
• Multi-Element devices
6 Active Flow Control 37
• Gurney Flaps / Micro Tabs
• Spoiler
• Boundary layer suction/blowing devices
A more detailed is found in [40] and [19]. AFC solutions inuence the ow eld around
the airfoil section of the blade and intent to either delay transition, decrease turbulence
or avoid ow separation. This usually entails drag reduction, lift enhancement, mixing
augmentation, heat transfer enhancement, and ow-induced noise reduction [19]. Un-
fortunately the benets of one eect usually include adverse eect on others. To nd
an optimized system, which might consist of several AFC elements, is the nal goal for
active ow control on wind turbines.
Active
Flow
Control
Triad
Devices & Actuators Controls & Sensors
Flow Phenomenon
LE / TE Flaps
MicroTaps
Vortex Generators
Synthetic Jets
Active Flexible Wall
Motor
Piezoelectric
MEMs
Fluidic
Conventional
Optical
MEMS
Neural Networks
Asaptive
Physical Model-Based
Dynamical Systems Based
Optimal Control Theory
Seperation Control
Adjust Sectional Lift
Drag Reduction
Noise Suppression
Figure 6.1: Feedback ow control triad (after [25]).
The main advantage over an intelligent blade pitch control is, that several AFC ele-
ments can be located on a blade independent from each other. This means, the ability
of individual AFC elements to mitigate fatigue loads or to reduce extreme loads is more
dierentiated. On the other hand, the control strategies get more complex. The use
of simple heuristic PID controlers might not be appropriate and more sophisticated
methods like neuro-fuzzy control approaches are necessary, in order to deal with the
non-linear aeroelastic wind turbine system. In addition to that, new sensors like strain
6 Active Flow Control 38
gauges or angle of attack senors to measure the control variables are necessary.
The use of AFC solutions on wind turbines is subject to current research. The focus
lies especially on trailing edge devices, either in form of active Flaps [4] [1] or in form
of Micro Tabs and active Gurney Flaps [57] [13].
In order to get further insight in the benets of active ow control concepts, it is
the overall goal of QBladeAE, to provide a simple method and a rst approximation
to investigate the inuences of dierent AFC solutions on wind turbines. The software
itself is presented in the following Part II.
7 QBladeAE
QBladeAE is used to investigate the behavior of wind turbine blades, which are
equipped with active ow control elements using an aeroelastic simulation. It is em-
bedded in the open-source software QBlade [34], which again is an extension of XFLR5
[11] (Figure 7.1). It is developed with the cross-platform C++ framework Qt, which
allows easy programming of applications with a graphical user interface. The features
of the program suite are:
XFLR5 Foil Design and XFoil analysis
• Direct geometric foil design
• XFoil direct analysis
• XFoil full and mixed inverse foil design
• (not included: wing and plane design)
QBlade Rotor and Turbine design:
• Blade design and optimization
• 360 polar extrapolation
• Turbine denition and simulation
• Rotor simulation
QBladeAE Active Flow Control Simulation
• Blade design with active elements
• Aerodynamic description of active elements
• Wind eld generator (beta)
• Aeroelastic simulation
7 QBladeAE 41
XFLR5Direct Foil Desing
XFoil Direct/Inverse Design
Figure 7.1: QBladeAE embedded in QBlade and XFLR5.
As mentioned above, QBladeAE works together with the aerodynamic routines Aero-
Dyn and the structural routines YawDyn. The original FORTRAN source code of
YawDyn was extended by an input/output handling and a control structure, in order
to simulate active elements. Therefore the name YawDynAE is chosen for the modied
version. QBladeAE allows the user to dene a blade structure, to handle aerodynamic
properties (with inherent XFoil calculations and 360-extrapolation) and to create all
necessary inputs for the NREL codes via a graphical user interface. It then automati-
cally calls the aeroelastic code externally and reads in the results when the simulation
is nished. All the results are visualized within QBladeAE and a binding to gnuplot
[56] allows the export of all graphs (currently beta status). Furthermore, QBladeAE
provides as well a GUI for creating TurbSim wind eld les, which works after the
same principle as described above. The workow of QBladeAE can be seen in Figure
7.2. After dening all necessary information, the output les for YawDynAE (yaw-
dyn.ipt),for AeroDyn (aerodyn.ipt), for all the used airfoils (airfoils.dat) and for the
active elements (active.ipt) are automatically generated. TurbSim les can be gener-
ated independently. Exemplary input les can be found in A.1.
7.1 Active Flow Control simulation
Two steps are necessary to simulate a blade which contains active ow control de-
vices. At rst, a standard blade has to be designed using the blade design module
7 QBladeAE 42
AeroDynYawDynAE
QBladeAEyawdyn.ipt
aerodyn.ipt
active.ipt
airfoils.dat
wind.wnd/.hh
yawdyn.plt
yawdyn.opt
element.plt
active.plt
TurbSim
wind.iptwind.wnd
wind.sum
Figure 7.2: Working principle of QBladeAE with input and output control to themodied NREL codes.
of QBlade. In a second step, user specied active elements can be added to the
blade. Dedicated blade sections can be declared active, which means these sections
have a variable aerodynamic representation during the runtime of the simulation. The
changing aerodynamic representation is realized by multiple airfoil polars for the active
elements. The multiple polars are used to describe the dierent operation points of
an active element. When using a trailing edge ap for example, the dierent polars
represent dierent ap angles or when using a boundary layer suction device, the polars
represent dierent suction rates. During the aeroelastic simulation, a controller is used
to determine the optimal actuator operation point and therewith the optimal polar.
Figure 7.3 shows a blade with several sections, where two of them are active.
By using the AeroDyn subroutines, this approach can be realized easily. AeroDyn
can already handle multiple airfoil polars. Originally this functionality is used for
Reynolds number dependent simulations, where the dierent polars represent the air-
foil characteristic under dierent ow conditions, but can as well be used for other
runtime variable airfoil characteristics like the eect of ailerons for example. Each
airfoil polar table has an ID, which is referred to as MulTabLocvariable (Multiple
7 QBladeAE 43
cl
Active Elements
α
cl
α
Figure 7.3: Blade with two active elements, which are represented by using severalairfoil polars.
Table Location) within AeroDyn. The multiple airfoil tables are stored in the airfoil
input les for AeroDyn. For each angle of attack there is more than one lift, drag and
eventually moment coecient. Each of the dierent polar sets represent a dierent
airfoil characteristic, dened by the MulTabLoc-value. This table ID is a numerical
value and by changing it, dierent airfoil tables can dynamically be selected within the
simulation. The table ID can be used to represent any variable airfoil characteristic,
from Reynolds numbers to aileron ap angles. If the desired MulTabLoc-variable is
not directly represented in the airfoil table, linear interpolation between the adjacent
tables is performed. If a desired MulTabLoc-value exceeds the range specied in the
airfoil table, there is no extrapolation but the most outer table is taken.
In order to determine the optimal actuator position during a simulation two control
approaches are implemented: a simple optimization loop and a PID controller. Figure
7.4 illustrates the used control circuit terminology. The aforementioned MulTabLoc is
synonymous with the actuator variable y(t) which is determined by the controller.
7.1.1 Optimization loop
In order to understand the working principle of the optimization loop, the communi-
cation between YawDyn and AeroDyn needs to be described rst. Figure 7.5 shows
the three main program loops of YawDyn: the time loop, the blade loop and the blade
element loop. The structural program YawDyn calls the aerodynamic program Aero-
7 QBladeAE 44
w(t) e(t) y(t)
z(t)
x(t)
Controler System
set point error
variableactuator
variable
disturbance
variable
control
variable
Figure 7.4: Schematic control circuit and control terminology.
Dyn once for each time step, blade and blade element. The AeroDyn routines return
the aerodynamic forces on a single element, which are the normal force on the element
(DFN), the tangential force (DFT) and the pitching moment (PMA). In order to deter-
mine the incident element velocity and the element forces, AeroDyn needs information
about the current state of the element. Within each call of AeroDyn, it calls four
routines in YawDyn again, which are GetVNVT (wind and blade element velocities),
GetRotorParams (rotor parameters like rotor speed, yaw angle, etc.), GetBladeParams
(blade parameters like azimuth angle) and GetElemParams (element parameters like
element pitch angle, radius, location and MulTabLoc). If the new parameters for all
the elements of all the blade are determined, the blade and rotor related parameters
are updated again.
Figure 7.5 shows as well the implementation of the optimization control loop. It is
positioned between the blade and the blade element loop. The loop's task is to nd the
optimal active element operation point, which is expressed by the MulTabLoc-variable.
The optimization loop runs the element loop several times, but each time with a dier-
ent MulTabLoc-value. If using a trailing edge ap for example, this is synonymous with
dierent ap angles. After the loop is nished, the best MulTabLoc-value (e.g. ap
angle) is determined, according to the control variable and the desired set point. These
parameters are given by the user, who has the possibility to choose the local blade ele-
ment forces DFN, DFT and PMA as control variables and to specify a desired set point
for this parameters (the set point needs to be determined in a previous simulation). It
is important to note, that one active ow control device (= active element) can spread
over several blade elements. Logically, there is only one MulTabLoc-variable for all the
7 QBladeAE 45
BLADE Loop
ELEMT Loop
TIME Loop
YawDynAE
MulTabLocCTRL Loop
Element dependant
control variables:
- DFN
- DFT
Figure 7.5: Implementation of the optimization loop for nding the optimal polarfor each active section (element dependent).
blade elements of one active element, like the trailing edge ap can only have one ap
angle for all the covered blade elements. It is important which blade element shall be
used to compare the control variable with the set point and only this blade element
will perform optimally. In reality this represents the sensor position.
The use of the control loop has disadvantages. First of all the computational eort
can be very high, especially when the step size of the loop is very small. Secondly only
element related parameters can be used as control variables (DFN and DFT). This
means blade or rotor dependent parameters, like blade ap deections or root bending
moments can not be controlled directly. Nevertheless one can argue, that the uctu-
ations of the element forces are sources of the load uctuations of the whole blade.
In addition to that, the use of an angle of attack sensor for example would only give
local element information as well. As the local element incident inow angle is directly
coupled to the aerodynamic performance of the element, this approach still has it's
7 QBladeAE 46
right to exist.
A major advantage of the control loop is that it does not require any controller
dependent information, as it always nds the optimal actuator operation point. This
is not realistic and the results of the optimization loop have to be seen as the optimal
potential of the active ow control device.
7.1.2 PID controller
Next to the optimization control loop, a simple PID controller is implemented as well.
Each active element on the blade has a PID controller with individual characteristic
parameters. The implementation is shown in Figure 7.6. Other than the control
loop, the PID controller is positioned in the time loop of YawDyn because the control
variables for the PID controller are no longer element dependent parameters, but blade
dependent. The user can choose the control variable to be either the blade ap rate, the
blade ap angle and the out-of-plane root bending moment. Unlike the optimization
loop, the computational eort using the PID controller is much smaller, as the desired
actuator variable (MulTabLoc-value) is directly determined by the PID controller only
once per time step.
The control circuit is shown in Figure 7.7. Next to the controller itself, there is
a rate limiter and a range limiter so that dierent actuator congurations can be
investigated (small actuator with big range or large actuator with small range). The
range limiter denes the minimal and maximal actuator operation point. Note that
this parameter depends as well on the multiple airfoil tables, specied in the airfoil
input le for AeroDyn. The actuator range should not exceed the MulTabLoc-values
in the airfoil les. The rate limiter controls the maximum speed of the actuator. In
order to simulate the time lag between sensor data acquisition and actuator control,
a time delay is included. For an easy implementation in the FORTRAN routines of
YawDyn, a simple array stores and keeps the desired actuator variables for a certain
time, before using them. Note that the controller delay can not be smaller than the
simulation time step and can only be expressed by integer multiples of the time step.
7 QBladeAE 47
BLADE Loop
ELEMT Loop
TIME Loop
YawDynAE
MulTabLoc
PID
Blade dependant
control variables:
- FlapRate
- FlapAngle
- AFMB
Figure 7.6: Implementation of the PID controller: one for each active element (bladedependent).
2
0
Blade Flap Rate β
1PID
Saturation
Rate limiter Delay
Actuator Angle
Set point
Figure 7.7: Exemplary control circuit for the PID controller using a trailing edgeap as actuator and the blade ap rate as control variable.
7 QBladeAE 48
Controller tuning
The PID uses three gain variables to determine the actuator variable, according to an
error between control variable and set point. The proportional gain Kp (more error →more reaction), the integral gain Ki (more error → faster reaction) and the dierential
gain Kd (faster error → more reaction):
yp(t) = Kp · e(t) (7.1)
yi(t) = Ki ·∫e(t)dt (7.2)
yd(t) = Kd ·de(t)
dt(7.3)
y = yp + yi + yd (7.4)
In case of only one active element per blade, the optimal gains can be found by applying
a tuning method like the Ziegler-Nichols tuning rule. The system reaction on a step
function is monitored while increasing the proportional gain Kp from zero (Ki = Kd =
0) until system gets instable. After the ultimate gain Ku and the oscillation period Tuare found, the gains of the PID controller gains are determined as follows:
Kp = 0.6Ku (7.5)
Ki = 2Ku
Tu(7.6)
Kd =KpTu8
(7.7)
A step system reaction can be simulated by using a inow wind speed drop dened in
a hub-height wind le for AeroDyn (see [27]).
A problem arises when using more than one active element per blade. As the control
variables are blade dependent, a change from a rst actuator, has an inuence on the
performance of a second. The heuristic tuning rule can not be applied anymore, as a
new set of gains for one active element would change the optimal gains for another.
To nd a global optimum for all the gains of all active elements, a simple sweep loop
could be applied, similar to the one described above. All the possible gain variations
of the controllers could be tested and compared with each other. This would yield in
an enormous computational eort. To reduce computational time, a more intelligent
7 QBladeAE 49
way, like an optimization strategy as described in [45]. In QBladeAE, there is no tuning
method built in and the user has to specify all controller gains manually.
7.2 Blade related simulation parameters
The geometrical blade representation in AeroDyn and the structural blade related
parameters for YawDyn need to be specied for a simulation. As described in 4.1 a
blade has one ap degree of freedom and is represented by a sti beam and a hinge-
spring model. QBladeAE automatically determines default values for the blade related
parameters, which are described in the following.
7.2.1 QBlade and NREL blade format
The blade representation in QBlade and in AeroDyn diers, as can be seen in Figure
7.8. In QBlade a blade is dened by sections, whose radius is measured from the
center of rotation. Additionally each section has a value for the chord length and
the twist angle. In AeroDyn the blade is represented by elements. The radius of the
element (RELM) is the distance between the blade-hub connection and the center of
the element. Each element then has a length (DR), a chord length and a twist angle.
DR5
RELM5
RH
ELM#5
x
y
x
yRleft,5
Rright,5
QBlade
AeroDyn
Figure 7.8: Dierent blade denition in QBlade and NREL format.
The single elements are two-dimensional extrusions and have constant properties
over their length. The values for the chord and the twist are interpolated between two
sections. The blade geometric conversion from QBlade in AeroDyn implicates that the
last section dened in QBlade can not be represented in AeroDyn.
7 QBladeAE 50
7.2.2 Blade mass
As mentioned above, the blade is assumed to have a uniform mass distribution. In order
to approximate the blade mass mblade as a function of the blade radius, the following
relation is used:
mblade(R) = 1.7 ∗R2.3 (7.8)
This empirical approximation is derived from real blade parameters shown in Figure
7.9.
0
5000
10000
15000
20000
25000
30000
20 30 40 50 60
blad
e m
ass
[kg]
blade length [m]
f(x) = 1.7*x2.3
real blade masses
Figure 7.9: Dierent blade masses over blade length and the used exponential ap-proximation function [51].
In YawDyn the blade mass is referred to as BM .
7.2.3 Blade center of gravity
The blade center of gravity is calculated by approximating the blade elements (QBlade
format) to be truncated cones, as illustrated in Figure 7.10. The two rectangles of the
cone have the same surface area An and An+1 as the two airfoils shapes and the same
chord lengths cn and cn+1. Hence the heights of the rectangles hn and hn+1 can be
7 QBladeAE 51
cn
cn+1
An
An+1
An
An+1
hn
cn+1
ELM#i COGi,real COGxi,approx
cn
hn+1
DRRnRnx
yz
Figure 7.10: Simplied geometric representation (rectangular cone) of a homoge-neous blade section for the calculation of the blade center of gravity.
determined. The element center of gravity cogx,i is referred to as RB within YawDyn
and can be calculated with
cogx,i = Rn +DR
2· cnhn + cnhn+1 + cn+1hn + 3cn+1hn+1
2cnhn + cnhn+1 + cn+1hn + 2cn+1hn+1
(7.9)
and the element volume is given as
Vi = DR · An + An+1
2(7.10)
By weighting the element center of gravity cogx,i with the element volume Vi, the total
blade center of gravity cogx,blade can be calculated using
cogx,blade =
∑NELMi=1 cogx,i · Vi∑NELM
i=1 Vi(7.11)
The uniform density of the blade can determined with
ρblade =mblade
Vblade(7.12)
Only the center of gravity in x-direction (spanwise direction) is computed. The dis-
tance between the center of gravity and the z-axis is neglected. Furthermore the COG
distance in y-direction is not needed, as YawDyn only reads the mass moment of inertia
for the rotation about the ap axis (y-axis).
7 QBladeAE 52
7.2.4 Blade mass moment of inertia
Similar to the assumptions made above, the mass moment of inertia Jflap about the
ap axis can be calculated using
Jflap =NELM∑i=1
mi · (cog2x,i + cog2z,i) =NELM∑i=1
Vi · ρblade · cog2x,i (7.13)
neglecting the distance of cogz,i. In YawDyn the mass moment of inertia is named
BLINER.
7.2.5 Torsional root spring constant
Another value needed by the structural model, is the torsional spring constant of the
equivalent apping hinge spring at the blade root. This value is named kβ in Figure
4.2 [28]. The parameter denes how sti or soft a blade is modeled. It is not possible
to give a general rule of thumb to determine this value and QBladeAE does not auto-
matically give a default value. In order to get a reasonable stiness value, the blade
deection under steady wind conditions can be taken as an indicator. According to
[29] a number of modern wind turbines show tip deections of about 8% of their blade
radius at wind speeds of 15m/s [2].
As YawDyn does not automatically generate an output for the tip deection, an
additional function is implemented in QBladeAE. The tip deection TDblade is derived
from the ap angle β of the blade and given as
TDblade = Rblade · sin β (7.14)
7.2.6 Dynamic stall parameters
Next to the structural blade related parameters, AeroDyn needs additional aerody-
namic information for the semi-empirical dynamic stall model (if enabled). QBladeAE
automatically calculates the necessary inputs, which are derived from the airfoil polars
and stored in the airfoil les for AeroDyn. These are
• Zero lift angle of attack: α0cl
• cn slope for zero lift: cn(α0cl)
7 QBladeAE 53
• cn at stall value for positive angle of attack: cn(α+stall)
• cn at stall value for negative angle of attack: cn(α−stall)
• Angle of attack for minimum cd: αcd,min
• Minimum cd value: cd,min
Most of the parameters can be directly derived from the polar tables. Only the
detection of the stall angles α+stall and α−stall needs further investigation. In order to
have a reliable stall indicator, the cl/cd over α is used. As can be seen in Figure 7.11
the rst stall eects occurring on the airfoil are visible as two peaks. This is when the
drag rapidly increases due to ow separation. This clear characteristic can be used
in an automated numerical algorithm to detect the positive and negative stall angles
reliably.
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
-180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180-100
-50
0
50
100
150
200
c n c l/c
d
α
cncl/cd
Figure 7.11: Dynamic stall related parameter using clcd-curve for automatically de-
tecting the critical static stall angle αstall.
To approximate the cn slope for zero lift cn(α0cl) a least square method is used, as
proposed in [27].
7 QBladeAE 54
Figure 7.12: QBladeAE active blade design module.
It shall be noted, that the active elements have no eect on the dynamic stall model.
As there is no proper aerodynamic formulation, it might be necessary to switch of the
dynamic stall model, when simulating a blade with active elements. If the actuators
are only positioned in the outer region of the blade, a negative blade pitch angle might
ensure, that there are no stall eects in the blade tip region.
7.3 Program modules
The following section gives a short overview of the four program modules used in
QBladeAE.
7.3.1 Blade design with active elements
In the active blade design menu (Figure 7.12), a new blade equipped with active el-
ements can be created. Therefore, a blade design from QBlade must be available.
Next to the 3D-visualization of the blade the geometric parameters are displayed in a
spreadsheet (QBlade format).
The Active Element group allows the denition of new active ow control devices.
The actuator parameters are specied in a separate dialog. These are
7 QBladeAE 55
Active Element Type The type of the active element species it's name and the type
of unit used for the representation of the MulTabLoc. For example a ap device
has the table ID unit [deg] for the ap angle, a boundary layer suction device has
the unit [m3/s] for the volumetric ow rate. Despite that, the correct setting of
the type and unit is not obligatory and does not change the calculation results.
Dimension The starting and end position of the active element expressed in the
QBlade blade format. Only dedicated sections can be used for the start and
end point of an active element.
Control Type The setting denes the control type of the active element. It can either
be LOOP for the optimization loop or PID for the PID controller. If there are
more than one active element, this setting is the same for all. The control types
can not be mixed.
Control Variable The control variable which can be selected depend on the controller
type. For the optimization loop the blade element dependent parameters DFN
and DFT can be selected. This is the element normal force and the element tan-
gential force. For the PID controller, the blade dependent parameters FlapRate,
FlapAngle and AFMB can be selected. This is the blade ap rate, the blade
deection angle β and the out-of-plane root bending moment (axial ap bending
moment). As can be seen in Equation 7.14 the blade deection angle and the
blade tip deection are synonymous.
Set point Set point for the controller. Depending on the selected control variable, the
units are either [kN ], [kNm], [deg] or [−].
Actuator speed This is the ± maximum actuator speed. The unit depends on the
used active element type. For a ap device the unit would be deg/s.
Min Max Range The operational range of the actuator can be limited. The range
specied here can not be bigger than the limits dened in the specic multiple
airfoil tables.
Kp, Ki, Kd The proportional, integral and derivative gain for the PID controller (only
PID).
7 QBladeAE 56
Delay The time delay in [ms] between sensor data acquisition and actuator response.
The values must be bigger than the simulation time step and is rounded to integer
multiples of the simulation time step (only PID).
Step size This is the step size for the optimization control loop expressed in active
element type dependent units. The smaller the step size the higher the controller
accuracy but the longer the simulation time (only LOOP).
Sensor position The sensor position denes the blade element from which the opti-
mization loop takes the control variable. If for example an active element covers
the blade element 10, 11 and 12 and the sensor position is set to element num-
ber 11, only this element will perform optimally. The control loop determines
the optimal actuator variable by comparing the control variable of this element
with the set point. The adjacent elements only follow the sensor element (only
LOOP).
7.3.2 Aerodynamic representation of active elements
In this module the variable aerodynamic characteristics of an active element are dened.
The inputs specied here are the basis for the AeroDyn airfoil les with multiple polar
tables. For each airfoil which is covered by an active element, a set of polars can be
added. Next to the multiple polars, the element needs a so-called mother-polar. This
is a 360-polar which was generated within QBlade. In addition to the mother-polar,
several child-polars can be added. These child polars can be either generated within
the XFLR5 module, or be imported. For the angle of attack range which is not specied
in the child-polar the coecients from it's 360-mother-polar are automatically taken
over. This assumption can be be made, as the airfoil characteristics are similar for very
high and very low angles of attack but the user can change all the polars manually as
well. Note that the polars can have dierent step sizes in the angle of attack range. By
adding child-polars to a mother-polar, the step size is automatically adapted, using a
liner interpolation. Once the child-polars are generated, the single values of the polar
can be edited manually, but single points shall not be deleted.
For each polar, which is added to the airfoil of the active element, a MulTabLoc-
values (table ID for the multiple airfoil tables) has to be specied.
7 QBladeAE 57
Figure 7.13: QBladeAE multiple aerodynamic polar module.
7.3.3 Wind eld simulation
The wind eld model of QBladeAE (Figure 7.14) is based on the Veers model [54] and
calculates correlated wind eld les for the x-direction. It includes an atmospheric
shear layer representation via a surface roughness length z0. For 3D-visualization the
the freely available Qt/OpenGl based C/C++ programming library QwtPlot3D is used.
As a turbulent inow eld can easily be generated with the more sophisticated Turb-
Sim, the wind eld module of QBladeAE is only implemented in a beta version and
will not be discussed any further.
7.3.4 Aeroelastic simulation
In the aeroelastic module of QBladeAE (Figure 7.15) all the necessary input for the
NREL routines AeroDyn ans YawDyn are provided. The user can select the available
(active) wings and dene the parameters for the aerodynamic and the structural model
in separate dialogs. The blade dependent parameters described above are automatically
updated once a new blade is selected.
After the inputs are dened, QBladeAE automatically starts YawDynAE as an ex-
ternal process and reads in the generated outputs, once the simulation is nished
7 QBladeAE 58
Figure 7.14: QBladeAE wind eld generator module (beta).
Figure 7.15: QBladeAE aeroelastic simulation module.
7 QBladeAE 59
successfully. QBladeAE automatically creates and organizes a folder structure on the
local hard drive.
For the visualization of the results, four dynamic 2D-graphs are used. They can be
used to plot the rotor parameters generated from YawDyn (yawdyn.plt) and the results
are stored in the output le from AeroDyn (element.plt). The blade element dependent
results can be either plotted over time (for a xed element) or over the radial blade
position (for a xed time step). In addition to that, the graphs can show the active
element dependent parameters as well. This is the actuator operation point and the
actuator speed over time.
A TurbSim dialog helps to dene TurbSim input les for the generation of turbulent
wind eld les.
8 Standard simulation
This part gives the results of an exemplary simulation of a wind turbine blade, which is
equipped with active elements. It is more meant to be a description of how to approach
the set up and use of a simulation with QBladeAE, rather than a complete scientic
investigation.
Before simulating a blade with active elements, a standard wind turbine simulation
is performed. The behavior of the aeroelastic model and the sensibility to changes of
specic key input parameters is shown. A nal test case is derived which will then be
used for a simulation using one and more than one active ow control elements.
8.1 Turbine and blade model
A ctive turbine of the 2.5MW class will be used for the simulation. For simplica-
tion the turbine operates at a xed rotational speed of 15rpm, as YawDyn can not be
used for a variable speed turbine simulation. The turbine is a conventional 3-bladed
upwind turbine and has a rotor diameter of 89m. The other parameters used for the
simulation are shown in Table 8.1.
8 Standard simulation 62
Table 8.1: Turbine parameters used for the simulation.
Parameter Symbol Value
Nominal power PN 2500kWRotor diameter D 89mBlade length lblade 43.3mHub radius Rhub 1.2mRotor tilt angle τ 4
Blade precone angle PC 2
Power regulation - pitchRotational speed (xed) n 15rpmTip speed ΩR 70m/sNumber of blades B 3Nominal wind speed vN 13m/sCut-in wind speed vin 3.5m/sCut-out wind speed vout 25m/sHub height HH 89m
The blades have a length lblade of 43.3m. The geometric shape of the blade created
with QBlade is shown in Figure 8.1. It has a limited twist of 12 in the root region
and is linearly tapered. The blade design is shown in Table A.2.
Figure 8.1: 3D view of blade
According to the assumptions described in 7.2, the structural blade parameters are
computed automatically by QBladeAE. The necessary blade parameters are the blade
mass, the distance of the center of gravity to the blade hinge axis and the mass moment
of inertia for the rotation about the ap axis. The values for the simulation are shown
in Table 8.2.
8 Standard simulation 63
Table 8.2: Blade structural parameters used for the simulation.
Parameter Symbol Value
Blade mass mblade 9840kgDistance to center of gravity cogx,blade 10.6mFlap mass moment of inertia Jflap 1.89 · 106 kgm2
Root spring stiness (no default value) FS 2.2 · 107 Nm/rad
As mentioned in 7.2.5, another blade structural input parameter is the torsional
spring constant of the equivalent apping hinge spring at the blade root. To determine
this value, a rule of thumb is used: as stated in [2], a usual blade deection at 15m/s
steady inow is about 8% of the blade length. As the turbine already pitches at this
wind speed, a parameter variation of the spring constant is performed at 13m/s. Figure
8.2 shows the blade tip deection for three dierent stinesses. The resulting tip
deection of about 3m (6.5% of blade length) seems reasonable and therefore a spring
constant of 2.2 · 107 Nm/rad is used.
0
1
2
3
4
5
6
7
8
9
0 10 20 30 40 50 60 70 80 90
Bla
de ti
p de
flect
ion
[m]
Time [s]
FS = 1.0 106 Nm/radFS = 2.2 106 Nm/radFS = 4.0 106 Nm/rad
Figure 8.2: Blade tip deection with three dierent root spring stinesses. Thewind inow is steady and constant over the whole rotor disk.
8 Standard simulation 64
8.2 Blade validation
The geometric blade denition in QBlade and QBladeAE diers, as described in 7.2.1.
To estimate the dierence between the two blade descriptions, a steady state rotor
performance simulation using QBlade is compared to the power output obtained with
several YawDynAE simulations. The comparison shows the inuence of the blade
translation from QBlade- to NREL-format. For the comparison, the turbine model in
YawDynAE has to be simplied: the tower and the blades are considered as rigid, so
no structural deections occur. The rotor tilt and pre-cone angle is set to zero. The
inow model in AeroDyn is set to the EQUIL-option (BEM), with a Prantl tip- and
hub-loss model. For the dierent wind speeds a totally uniform and constant inow
(no shear layer) is assumed. As can be seen in Figure 8.3 the results match up very
well.
0
500
1000
1500
2000
2500
3000
0 5 10 15 20 25
P [k
W]
v [m/s]
QBladeYawDynAE
Figure 8.3: Rotor power over wind speed calculated with QBlade and QBladeAE.
Note that modeled turbine is pitch regulated. QBlade automatically calculates the
necessary blade pitch angles θb to keep the power output at the rated nominal power
level of 2500kW if the wind speed exceeds the nominal wind speed of 13m/s. This pitch
angle is then used as well in YawDynAE. Figure 8.4 shows the used pitch angle over
the wind speed.
8 Standard simulation 65
0
5
10
15
20
25
30
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
θ b [d
eg]
v [m/s]
Blade pitch
Figure 8.4: Blade pitch over wind speed for power regulation.
8.3 Dynamic stall eects
As mentioned in 3.2.3 AeroDyn features a dynamic stall model. The use of active
elements change the aerodynamic behavior of the airfoil and thus the semi-empirical
derived dynamic stall model might not be valid any more. To avoid any stall related
phenomena for the active simulation, the blade is pitched towards lower angles of
attack. This ensures that the outer region of the blade (where the active elements will
be positioned) operate in the attached-ow region during the whole simulation. Figure
8.5 shows the inuence of the dynamic stall model with a pitch angle of θb = 0. The
bottom graph shows the angle of attack distribution over the blade at a specic instant
of time. It indicates that also the blade outer regions operate under stalled conditions,
as the angle of attack already exceeds the value for the static stall angle of attack.
Logically, the eect of the dynamic stall phenomenon can be seen in the top graph,
where the blade deection over the simulation time is shown. The eect of dynamic
stall leads to higher blade deections, as the maximum aerodynamic forces which occur
are higher. The slopes of the blade deection curve are also higher due to the rapid
lift break-in after maximum lift.
8 Standard simulation 66
0.0
5.0
10.0
15.0
20.0
25.0
30.0
35.0
40.0
45.0
0 5 10 15 20 25 30 35 40 45
/S
ymbo
l a
[deg
]
radial position [m]
STEADYBEDDOES
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
0 10 20 30 40 50 60 70 80 90
blad
e tip
def
lect
ion
[m]
Time [s]
STEADYBEDDOES
Figure 8.5: Inuence of the dynamic stall model on the blade tip deection overtime and the angle of attack over radial position for the blade with apitch angle of θb = 0.
On the contrary, a simulation with a blade pitch angle of θb = 5 almost eliminates
the eects of dynamic stall. As can be seen in Figure 8.6, the angles of attack along
the blade are reduced and the outer blade region now operates in the attached ow
region. The power output for the pitched blades is logically reduced.
The eect of a the apping exible ap is currently only modeled by jumping from
one static polar table to the other. As can not be foreseen which inuence the apping
has on the dynamic stall model, it is switched o for all further simulations. Although
the pitched blade ensures a exible ap operation in the attached ow region for most
of the time, it has to be noted, that the apping itself generates unsteady eects, just
like blade pitching in the attached ow region (3.2.3). These dynamic eects are not
modeled by the software.
8 Standard simulation 67
0.0
5.0
10.0
15.0
20.0
25.0
30.0
35.0
40.0
45.0
0 5 10 15 20 25 30 35 40 45
/S
ymbo
l a
[deg
]
radial position [m]
STEADYBEDDOES
-0.5
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
0 10 20 30 40 50 60 70 80 90
blad
e tip
def
lect
ion
[m]
time [s]
STEADYBEDDOES
Figure 8.6: Inuence of the dynamic stall model on the blade tip deection overtime and the angle of attack over radial position for the blade with apitch angle of θb = 5.
8.4 Yawed turbine
One of the key features of YawDyn is to model a yawed turbine and the opportunity
shall not be missed to show the inuence of dierent yaw conditions. YawDyn provides
dierent models, which include free yaw, xed yaw or yawing with a constant rate.
As the necessary parameters like yaw-stiness, -damping and -friction are dicult to
evaluate, the turbine is only modeled with a xed yaw model. Only dierent yaw angles
are regarded. Figure 8.7 shows the eect of three dierent yaw angles on the blade tip
deection. However, for the further simulations a yaw angle of 0 is used.
8 Standard simulation 68
1
1.5
2
2.5
3
3.5
4
0 10 20 30 40 50 60 70 80 90
blad
e tip
def
lect
ion
[m]
time [s]
0° yaw angle2° yaw angle
10° yaw angle
Figure 8.7: Blade tip deection for dierent yaw angles.
8.5 Wind eld
The GUI interface for TurbSim is used to create a turbulent inow eld for the further
simulations. An IEC Kaimal -spectrum is used, with a turbulence intensity of 15%. The
mean wind speed is set to 13m/s at a hub height of 89m. The shear layer is simulated
with a power law exponent PLExp of 0.2 and a surface roughness length of z0 = 0.03.
This wind le will be used for all the further simulations. The wind speed time series
at the hub height is shown in Figure 8.8.
8 Standard simulation 69
7
8
9
10
11
12
13
14
15
0 10 20 30 40 50 60 70 80 90
win
d sp
eed
[m/s
]
time [s]
Figure 8.8: Turbulent wind speed time series in x-direction at a hub height of 89mand a mean wind speed of 13ms.
8.6 Baseline simulation
After combing all the aforementioned model inputs and assumptions, a baseline test
case is derived. All the simulations will be performed using these default values, unless
stated otherwise. The basic results of the simulation can be seen in Figure 8.9. The
simulation at a constant rotational speed includes atmospheric turbulence and shear
layer, blade-tower interaction, skewed inow with an unsteady inow model (GDW),
tilted hub with pre-coned and pitched blades. Note that the out-of plane bending
moment has the same curve progression as the blade tip deection. This is due to the
simple structural blade representation in YawDyn (4.1).
8 Standard simulation 70
1
2
3
4
0 10 20 30 40 50 60 70 80 90
blad
e de
flect
ion
[m]
time [s]
1000
2000
3000
o-p
bend
ing
mom
ent [
kNm
]
1000
2000
3000
pow
er [k
W]
200
250
300
350
thru
st [k
N]
6
8
10
12
14
16
win
d sp
eed
[m/s
]
Figure 8.9: Results of baseline simulation.
9 AFC simulation
This chapter describes a simulation using a blade, which is equipped with active ele-
ments. As an example, a form-exible trailing edge device (exible ap) is used. As
mentioned above, this work does not claim to be a complete scientic investigation.
The following investigation shall rather demonstrate the performance of QBladeAE.
The exible airfoil structure was intensively investigated by Smart Blade GmbH in
the past and experimental wind tunnel data are available [41]. The aps are actuated
with pneumatic "muscles", which contract by applying air pressure. The aps can
be deected in both, upwards direction (negative deection) and downwards direction
(positive deection). The baseline airfoil is a DU-96-W-180 and the shapes for for
dierent ap deections are shown in Figure 9.1 and 9.2. The test wing in the wind
tunnel experiment had a chord of 0.6m with both the rigid and exible trailing edge
segments thus achieving a similar Re of 1, 300, 000.
Figure 9.1: Overlapping airfoil contours for positive deection. Red: The originalDU-96-W-180 airfoil; Green: slightly deected exible ap; Red: fullydeected exible ap [41].
The apping part has a chord length of 25% and the maximum apping range is
±20. The maximum actuator speed is 20/s.
9 AFC simulation 72
Figure 9.2: Overlapping airfoil contours for negative deection. Red: The originalDU-96-W-180 airfoil; Green: slightly deected exible ap; Red: fullydeected exible ap [41].
A ap deection induces changes in the circulation, and thus the aerodynamic charac-
teristics of the airfoil areb changed. In the wind tunnel measurements at the GroWiKa
of TU Berlin the performance of the exible-form airfoil was determined. The eect
on the lift and drag coecients can be seen in Figure 9.3. The overall change in the
lift coecient is about ∆cl = 2. On the other side, any ap deection increases the
drag value.
-2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-15 -10 -5 0 5 10 15 20
c l
α
0.0
0.1
0.1
0.2
0.2
0.2
-15 -10 -5 0 5 10 15 20
c d
α
neutral positionslight neg. deflection
full neg. deflectionslight pos. deflection
full pos. deflection
Figure 9.3: Lift and drag coecient over angle of attack cl(α) for exible ap atfour ap angles.
Using the airfoil extrapolation method of QBlade (3.2.2), 360-polars were derived,
which can be seen in Figure 9.4. The lift and drag values for a dierent ap deections
9 AFC simulation 73
are assumed to be the same for angles of attack higher than 28 and lower than −11.
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
-150 -100 -50 0 50 100 150
c l
α
neutral positionslight neg. deflection
full neg. deflectionslight pos. deflection
full pos. deflection
Figure 9.4: Extrapolated 360 cl-polar.
9.1 Flap parameter study
In this section a blade with dierent ap congurations is investigated. If not explicitly
mentioned, the PID controller is used for all simulations, with the ap rate as control
variable and a set point of zero. The default performance parameters of exible ap
are a range of ±20 and a maximum actuator speed of 20/s. A small parameter inves-
tigation shows the inuence of dierent ap positions, ap sizes and number of aps
as well as an investigation of the operational range, the maximum speed and the delay
of the actuator. All the controller gains are found manually for each conguration and
therefore an optimal controller tuning can not be ensured. However, the actuator rate
and range limit provides a robust controller tuning in most cases.
For the comparison of the dierent potentials for load reduction, the standard devi-
9 AFC simulation 74
ation σ is taken as indicator. It is dened as:
σ =
√√√√ 1
n− 1
n∑i=1
(xi − x)2 (9.1)
The load reduction is expressed as percentage of the standard deviation σ0 of the base-
line simulation:
χ = 1− σiσ0
(9.2)
For the parametric investigation, the outer part of the blade is divided in nine equidis-
tant parts of 1.5m, as can be seen in Figure 9.5. These sections, and combinations of
each other, will be modeled as active in the following simulations.
123456789
DU-96-W-180
43.25m1.5m
Figure 9.5: Blade with 9 equidistant outer sections.
According to the NREL blade format, the radial position of the sections is shown in
Table 9.1. According to the convention, r is the distance between blade-hub connection
and the center of the element and DR is the element length.
Some of the following simulations are as well performed to compare them with the
results obtained in [2], in which the performance of trailing edge aps was investigated.
The comparisson is used as well as a possibility to validate the obtained results.
9.1.1 Flap positions
Several simulations are performed using a single 1.5m ap at dierent spanwise posi-
tions. All the nine positions investigated independently. This investigation shows the
optimal spanwise position for the ap. The load reduction χ on the out-of plane bend-
ing moment can be seen in Table 9.2. The highest reduction with about 14% have the
9 AFC simulation 75
Table 9.1: Radial position of the 9 possible active elements.
AE# r [m] r/R [%] DR [m] DR/R [%]
1 41.75 96.5 1.5 3.52 40.25 93.1 1.5 3.53 38.75 89.6 1.5 3.54 37.25 86.1 1.5 3.55 35.75 82.7 1.5 3.56 34.25 79.2 1.5 3.57 32.75 75.8 1.5 3.58 31.25 72.3 1.5 3.59 29.75 68.8 1.5 3.5
active elements at position AE# 2 and AE# 3. This is around r/R = 92% of the blade
length. The more outer sections have a smaller chord and the more inner elements a
smaller radial velocity, which both lower their inuence on the total blade forces.
Table 9.2: Load reduction for 1.5m ap at dierent radial positions.
AE# σAFMB [kNm] χ [%]
Default 274.6 0.01 239.6 12.72 237.2 13.73 237.4 13.54 238.3 13.25 239.7 12.76 241.3 12.17 242.9 11.68 244.8 10.99 246.8 10.1
Figure 9.6 points out the relation between ap position and load reduction. Com-
pared to the investigation in [2], the results correspond well.
9 AFC simulation 76
9
10
11
12
13
14
65 70 75 80 85 90 95 100
redu
ctio
n in
out
-of-
plan
e be
ndin
g m
omen
t χ [%
]
r/R [%]
Figure 9.6: Load reduction of a single ap with a length of 1.5m at dierent radialpositions.
9.1.2 Flap size
This section investigates the inuence of the ap size. Three sections or more at a time
are modeled as one active element. This yields in a ap length of 4.5m, 9m, 10.5m,
12m and 13.5m respectively. The positions of the aps are varied again, and nally
all sections form AE# 1 − 9 are simulated as one single ap. The reduction in χ on
the out-of plane bending moment is listed in Table 9.3. It can be seen, that the load
reduction gets bigger, with increasing ap length, although the reduction potential is
not linearly growing. For small aps, the reduction is around 28%, the medium sized
aps reach 38%− 44% and the big aps up to 50%. In accordance with the results of
the previous section, an outer position of the ap is preferable.
Finally Figure 9.7 summarizes the results and shows the load reduction of the ap
wise root bending moment for all previous simulations. A comparison with [2] again
shows the correlation between the results. Although the potentials for load reduction
are of the same magnitude, it has to be noted that the denition of the load reduction
given in this document, diers from the one in the cited report.
9 AFC simulation 77
Table 9.3: Load reduction for dierent ap lengths and dierent radial positions.
AE# DR [m] σAFMB [kNm] χ [%]
Default - 274.6 0.0123 4.5 193.7 29.5456 4.5 196.9 28.8789 4.5 205.6 25.1
123456 9.0 152.8 44.4456789 9.0 169.1 38.41234567 10.5 135.2 50.82345678 10.5 136.9 50.13456489 10.5 145.8 46.912345678 12.0 138.4 49.623456789 12.0 137.6 49.9123456789 13.5 131.4 52.1
To have an understanding of how the out-of-plane bending moment is inuenced by
the active elements, Figure 9.8 shows the out-of-plane root bending moments for the
baseline simulations and the one for the single 13.5m ap, with the maximum load
reduction of about 52%.
9 AFC simulation 78
0
10
20
30
40
50
0 2 4 6 8 10 12 14
redu
ctio
n in
out
-of-
plan
e ro
ot b
endi
ng m
omen
t χ [%
]
flap length [m]
Figure 9.7: Load reduction of single aps with dierent lengths and dierent radialpositions.
600
800
1000
1200
1400
1600
1800
2000
2200
2400
0 10 20 30 40 50 60 70 80 90
out-
of-p
lane
ben
ding
mom
ent b
lade
1 [k
Nm
]
time [s]
baseline13.5m flap
Figure 9.8: Out-of-plane bending moment for baseline simulation and the single13.5m aps with maximum load reduction.
9 AFC simulation 79
9.1.3 Actuator speed and range
The software makes it easy as well, to change the operational range and the maximum
speed of the actuator. The standard exible ap described above has a range of ±20
and a speed of 20/s. Generally it can be said, that if a ap has a high range it needs
as well a certain speed, in order to make use of this range. If the ap is too slow, it
is not able to reach it's extreme positions within one rotor revolution. In this case the
actuator speed is the limiting parameter. On the other hand, if the ap is too fast, the
additional gain is low, as the ap range becomes the limiting factor. Figure 9.9 shows
the inuence of dierent actuator speeds for the the 4.5m ap (123) and the 9m ap
(123456). It shows, that for the specic simulation conditions (15rpm), the ap speed
of 20/s is sucient.
0
5
10
15
20
25
30
35
40
45
50
0 10 20 30 40 50 60
redu
ctio
n in
out
-of-
plan
e be
ndin
g m
omen
t /S
ymbo
l c
[%]
actuator speed [deg/s]
4.5m flap (123)9m flap (123456)
Figure 9.9: Inuence of dierent actuator speeds on the load reduction for two apcongurations.
The overall ap length has a big inuence on the parameters as well. The bigger the
ap, the the smaller the range of the ap has to be. In order to investigate the inuence
of the two actuator parameters in combination, the 9m ap (123456) was exemplary
modied in the range and speed. Figure 9.10 shows the dependencies of the absolute
ap range (it is assumed that the ap can be equally deected in positive as well as
negative direction) and the ap speed on the load reduction. It can be seen, that a
9 AFC simulation 80
ap with a limited ap range of ±15 can still achieve a reduction of 45% compared
to the maximal possible 49%, but provides a simpler structural design.
0 5 10 15 20 25 30 0 5
10 15
20 25
30 35
40 0 5
10 15 20 25 30 35 40 45 50
χ
efficient flap configuration
actuator speed [deg/s]
abs. actuator range [deg]
χ
Figure 9.10: Inuence of dierent actuator ranges and speeds on the load reductionfor a 9m ap conguration.
The simulation above makes the assumption that the range of deection in positive
and negative direction is the same. In other words, a ap with the proposed absolute
deection of 30, has a deection of ±15. For further structural simplications of
the ap, a non-equal deection range can be taken into account as well. Figure 9.11
indicates, that the ap deects slightly more in positive direction, than in negative
direction.
These results are only valid for the examined case and all the other relevant cases
have to be examined as well, but it shows the potential of the software. The short
computation time of about 10s per simulation makes it possible to investigate sev-
eral congurations and parameters and the software helps to nd an optimal actuator
conguration, where the actuator size is harmonized with the actuator rate and the
actuator speed.
9 AFC simulation 81
-15
-10
-5
0
5
10
15
0 10 20 30 40 50 60 70 80 90
flap
angl
e [d
eg]
time [s]
Figure 9.11: Flap angle for 9m ap over simulation time.
9.1.4 Sensor delay
In all the previous simulations it was assumed that the measured control variable is
processed by the controller without delay and the actuator can respond instantaneously.
The introduction of a delay between sensor data acquisition and the ap response has
a signicant inuence on the potential for load reduction. Dierent time delays were
simulated for the 9m ap with the sections AE# 1−9. Figure 9.12 shows the decrease
of load reduction potential χ with increasing time delay. A delay of only 5ms reduced
the potential by half.
9.1.5 Multiple aps
Until now, only adjacent sections have been apped. The question arises what happens,
if a ap is split into two or more parts. The simulation of multiple and independent
aps on one wing with the PID controller, is not scope of this work. The manual
controller tuning for two and more aps is demanding and using the gains from the
previous simulations did not yield in a higher load reduction than for a single ap.
This is not consistent with the results, which were found in other investigations, like
[2]. With the use of multiple aps, the total length of the apping sections could be
9 AFC simulation 82
5
10
15
20
25
30
35
40
45
0 5 10 15 20 25 30 35 40
redu
ctio
n in
out
-of-
plan
e be
ndin
g m
omen
t χ [%
]
controller delay [ms]
Figure 9.12: Load reduction over controller time delay.
reduced to obtain the same load reduction as for a big single ap.
9.2 Optimization loop
As described above, it is not trivial to nd the optimal controller gains for multiple
aps. Each actuator inuences the blade dependent control parameters individually
(ap rate, ap angle, out-of-plane bending moment) but they are not independent from
each other. A more sophisticated controller or a (intelligent) tuning sweep has to be
used.
Another approach to investigate the advantage of multiple active elements is the use
of the optimization loop. As described in 7.1.1, a local blade element control variable
is used for the controller. In this case it is the local blade element normal force DFN .
The loop does not represent a real controller, as it simply calculates all possible ap
angles for a time step and then chooses the optimal one. Another drawback is, that
keeping the local blade element force constant does not mean, that the uctuation of
the out-of plane bending moment is minimized. On the other hand, a concept of several
independent actuators, which only have local control variables has advantages as well.
9 AFC simulation 83
Multiple actuators with more simple integrated sensor-controller units can be used as
a modular concept for active ow control. The sensor could measure, for example, the
local angle of attack, which couples the local inow velocity directly to a trailing edge
ap angle. Then each ap works for it's own, yielding in a less powerful but also less
complex and modular system.
To show the potential of this concept, a simulation with all sections form AE# 1−9
being individually active is performed. This represents a conguration, in which each
actuator has it's own sensor and controller. The load reduction is 44%, as can be seen
in Table 9.4. In contrary, a single ap including the sections AE# 1− 9 with only one
sensor at section AE# 3, shows the same eect.
Table 9.4: Load reduction for individual sections using the optimization loop.
AE# Ctrl. Type DR [m] σAFMB [kNm] χ [%]
Default - - 274.6 0.0123456789 PID 13.5 131.4 52.1
1-2-3-4-5-6-7-8-9 LOOP (9 sensors) 13.5 155.1 43.6123456789 LOOP (1 sensor @ AE# 3) 13.5 155.5 43.7
In Figure 9.13 the element normal force for active element number AE# 3 is shown
for the three congurations listed above. It can be seen, that the optimization loop
keeps the element force at a constant value of about 2900kNm. Next to the normal
force of the baseline simulation, the element force from the PID controlled single ap
is shown as well. It has to be noted, that the optimization loop does not represent a
real controller. The results have to be seen as ideal. The reason for the uctuation of
the normal force in Figure 9.13 is the step size of the control loop. For the simulation
the step size is set to 1. If the step size of the loop is small enough, there would be
no uctuation at all.
The approach of keeping the local blade element force constant, shows less potential
in the reduction of the blade root bending moment than blade PID control approach.
Additionally, there is not a big gain in using multiple sensors for individually controlled
elements. This result has to be examined in more detailed investigations.
9 AFC simulation 84
1000
1500
2000
2500
3000
3500
4000
4500
5000
0 10 20 30 40 50 60 70 80 90
blad
e el
emen
t nor
mal
forc
e D
FN
of A
E#3
[kN
m]
time [s]
baselinePID: 13.5m flap
LOOP: 9x1.5m flap
Figure 9.13: Local blade element force DFN at AE# 3 for the baseline, the singlePID controlled 13.5m ap and the multiple individually optimizationloop controlled aps.
The parametric investigation presented here shall only demonstrate the potential
and the limitation of the software. It is the initial step for further, more detailed
investigations.
10 Suggestions for future research
The current software provides the possibility to estimate the performance of dierent
active ow control devices on wind turbine blades. Now dierent active ow concepts
have to be investigated in more detail and compared with each other.
In a second step, a better model representation has to be implemented. Especially the
oversimplied structural wind turbine and blade model of YawDyn has to be improved.
The following list gives an overview about possible further steps:
Structural model: To simulate the eects of AFC elements on wind turbines in more
detail, an advanced structural model has to be used. An enhanced structural
binding to AeroDyn is FAST [22]. The medium complexity structural model is as
well provided by NREL and features and the possibility to include wind turbine
operational control. Useful information on FAST can be found in [23], [20] and
[21]. An extension called CurveFAST adds a blade torsional degree of freedom
and is described in [30].
Aerodynamic model: As the implementation of actuators for active ow control
changes the physic eects on the ow eld, the aerodynamic models have to
be adapted accordingly. If trailing edge aps are used, for example, the imple-
mentation of an unsteady aerodynamic model for the attached ow region as
described in [1] or [4] might be necessary.
Control: A more advanced control strategy than the simple PID controller has to be
developed and/or the controller tuning has to be improved in order to simulate
multiple aps on a blade.
Noise: The additional noise, which created by the actuators has to be investigated.
Scope of investigations The whole range of the dierent wind turbine states of op-
eration and load cases has to be investigated.
10 Suggestions for future research 86
Validation: As only few other investigations have been made, the results have to be
validated against experimental data. Until now, there is no possibility to validate
the results properly.
11 Conclusion
A software for the preliminary investigation of dierent active ow concepts for wind
turbines has been developed. It is based on the open-source codes XFLR5 and QBlade.
The user friendly extension QBladeAE features a graphical user interface for the use of
the NREL aeroelastic design codes AeroDyn and YawDyn. The user can easily design
wind turbine blades on which dierent active elements can be positioned. The dier-
ent aerodynamic characteristics of the actuators for active ow control are represented
by their individual two-dimensional airfoil polars. The NREL FORTRAN routines
were modied and two approaches for the control of the dierent active elements were
implemented: a simple PID controller and an optimization loop.
The working principle and the limitation of the software was shown with the simula-
tion of a form exible trailing edge (exible ap). The results show once again the
potential of trailing edge devices for load reduction on wind turbine blades: an ideal
ap of 13.5m length on a 43m long blade reduced the standard deviation of the root
apwise bending moment by 52%. The inuence of dierent ap size, ap positions,
numbers, ap speed and ap range has been investigated as well.
Bibliography
[1] P. B. Anderson. Advanced Load Alleviation for Wind Turbines using Adap-
tive Trailing Edge Flaps. Sensoring and Control. PhD Report. Roskilde: Riso
National Laboratory for Sustainable Energy DTU, 2010.
[2] P. B. Anderson. Load Alleviation on Wind Turbine Blades using Variable Airfoil
Geometry (2D and 3D study). Sensoring and Control. M.Sc. Thesis. Lyngby:
Technical University of Denmark, 2005.
[3] C. Bak et al. Airfoil Characteristics for Wind Turbines. Roskilde: Riso National
Laboratory for Sustainable Energy DTU, 2010.
[4] T. Barlas and G. van Kuik. Aeroelastic Modelling and Comparison of Advanced
Active Flap Control Concepts for Load Reduction on the Upwind 5MW Wind
Turbine. Delft: Delft University of Technology, 2009.
[5] E. A. Bossanyi. GH Bladed Theory Manual. Hamburg: Garrad Hassan and
Partners Ltd, 2006.
[6] M. L. Buhl. WTPerf User's Guide. for Version 3.1. Golden, Colorado: National
Renewable Energy Laboratory, 2004.
[7] M. L. Buhl, Jr. and A. Manjock. A Comparison of Wind Turbine Aeroelastic
Codes Used for Certication. Golden, Colorado: National Renewable Energy Lab-
oratory, 2006.
[8] T. Burton et al. Wind Energy Handbook. Bans Lane, Chichester: John Wiley
& Sons, Ltd, 2001.
[9] L. W. Carr. Progress in the Analysis and Prediction of Dynamic Stall. In:
Journal of Aircraft Vol. 25 No. 1 (1988), pp. 617.
[10] P. E. H. Currin and J. Long. Horizontal Axis Wind Turbine Free Wake Model
for AeroDyn. Klamath Falls, Oregon: Oregon Institute of Technology, 2009.
Bibliography 89
[11] A. Deperrois. XFLR5 Gudieines. 2010. url:
http : / / sourceforge . net / projects / xflr5 / files / xflr5 % 20v6 . 03 % 20
Beta/Guidelines_Feb_2011.pdf/download (visited on 02/01/2011).
[12] M. Drela. XFOIL 6.94 User Guide. Cambridge, Massachusetts: Massachusetts
Institute of Technology, 2001.
[13] O. Eisele et al. Experimental Investigation of Dynamic Load Control Strategies
Using Active Microaps on Wind Turbine Blades. In: Conference Proceedings.
Ed. by European Wind Energy Association. Brussels 2011.
[14] European Wind Energy Association. Wind Energy The Facts. A guide to the
technology, economics and future of wind power. Sterling: Earthscan, 2009.
[15] Garrad Hassan and Partners Ltd. GH Bladed Price List. Silverthorne Lane
2009.
[16] R. Gasch and J. Twele.Wind Power Plants. Fundamentals, Design, Construction
and Operation. Berlin: Solarpraxis, 2002.
[17] A. C. Hansen. Yaw Dynamics of Horizontal Axis Wind Turbines. Final Report.
Golden, Colorado: National Renewable Energy Laboratory, 1992.
[18] M. O. Hansen. Aerodynamics of Wind Turbines. 2nd ed. London: Earthscan,
2008.
[19] S. J. Johnson, C. Dam, and D. E. Berg. Active Load Control Techniques for
Wind Turbines. Sandia Report SAND2008-4809. Albuquerque: Sandia National
Laboratories, 2008.
[20] J. M. Jonkman. Dynamics Modeling and Loads Analysis of an Oshore Float-
ing Wind Turbine. Technical Report NREL/EL-500-41958. Golden, Colorado:
National Renewable Energy Laboratory, 2007.
[21] J. M. Jonkman. Modeling of the UAE Wind Turbine for Renement of FAST.
Technical Report NREL/EL-500-34755. Golden, Colorado: National Renewable
Energy Laboratory, 2003.
[22] J. M. Jonkman and M. L. Buhl. FAST User's Guide. Technical Report NREL/EL-
500-38230. Golden, Colorado: National Renewable Energy Laboratory, 2005.
Bibliography 90
[23] J. M. Jonkman and M. Buhl. New Developments for the NWTCs FAST Aeroe-
lastic HAWT Simulator. Conference Report NREL/CP-500-35077. Golden, Col-
orado: National Renewable Energy Laboratory, 2004.
[24] S. Kanev and T. van Engelen. Exploring the Limits in Individual Pitch Control.
In: Conference Proceedings. Ed. by European Wind Energy Conference 2009.
Marseille 2009.
[25] L. D. Kral. Active Flow Control Technology. Technical Brief. St. Louis: ASME
Fluids Engineering Division, Washington University, 1998.
[26] M. A. Lackner and G. van Kuik. A comparison of smart rotor control approaches
using trailing edge aps and individual pitch control. In:Wind Energy 13 (2010),
pp. 117134.
[27] D. Laino and A. C. Hansen. User's Guide to the Wind Turbine
Aerodynamics Computer Software AeroDyn. NWTC Design Codes.
Salt Lake City, Utah: Windward Engineering, LC, 2002. url:
http :/ / wind . nrel . gov / designcodes / simulators / aerodyn/ (visited on
02/03/2011).
[28] D. Laino and A. C. Hansen. User's Guide to the Wind Turbine Dynamics Program
YawDyn. NWTC Design Codes. Salt Lake City, Utah: Windward Engineering,
LC, 2003. url: http://wind.nrel.gov/designcodes/simulators/yawdyn/
(visited on 05/26/2005).
[29] T. Larsen, A. Hansen, and T. Buhl. Aeroelastic eects of large blade deections
for wind turbines. In: Conference Proceedings. Ed. by Delft Conference The
Science of making Torque from Wind. 2004.
[30] S. M. Larwood. Dynamic Analysis Tool Development for Advanced Geometry
Wind Turbine Blades. Dissertation. Davis: University of California, 2009.
[31] J. G. Leishman and T. S. Beddoes. A Semi-Empirical Model for Dynmaic Stall.
In: Journal of the American Helicopter Society 34(3) (1989), pp. 317.
[32] J. G. Leishman. Challenges in Modeling the Unsteady Aerodynamics of Wind
Turbines. In: Conference Proceedings. Ed. by 21st ASME Wind Energy Sympo-
sium. Reno 2002.
[33] J. Mann. Wind eld simulation. In: Prob Engng Mech Vol. 13, No. 4 (1998),
pp. 269282.
Bibliography 91
[34] D. Marten et al. Integration of a wind turbine blade design tool in
XFOIL/XFLR5. In: Conference Proceedings. Ed. by DEWI 2011. 2010.
[35] D. Marten. Extension of an Aerodynamic Simulator for Wind Turbine Blade De-
sign and Performance Analysis. Diploma Thesis. Berlin: Technische Universitaet
Berlin, 2010.
[36] D.-P. Molenaar. Cost-eective design and operation of variable speed wind tur-
bines. Closing the gap between the control engineering and the wind engineering
community. Phd Thesis. Delft: Technische Universiteit Delft, 2003.
[37] B. Montgomerie.Methods for Root Eect, Tip Eects and Extending the Angle of
Attack Range to +-180deg with Application to Aerodynamics for Blades on Wind
Turbines and Propellors. Stockholm: Swedish Defence Research Agency, 2004.
[38] P. J. Moriarty and A. C. Hansen. AeroDyn Theory Manual. NWTC Design
Codes. Golden, Colorado: National Renewable Energy Laboratory, 2011. url:
http :/ / wind . nrel . gov / designcodes / simulators / aerodyn/ (visited on
02/03/2011).
[39] P. Passon et al. OC3Benchmark Exercise of Aero-elastic Oshore Wind Turbine
Codes. In: Journal of Physics (2007), Conference Series 75.
[40] G. Pechlivanoglou, C. Nayeri, and C. Paschereit. Performance optimization of
wind turbine rotors with active ow control. In: Conference Proceedings. Ed. by
Proceedings of ASME IGTI Turbo Expo 2011 ASME/IGTI June 6 - 10. Vancou-
ver 2011.
[41] G. Pechlivanoglou et al. Active aerodynamic control of wind turbine blades with
high deection exible aps. In: Conference Proceedings. Ed. by 48th AIAA
Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace
Exposition 4 - 7 January 2010. Orlando 2010.
[42] J. Peeters. Simulation of Dynamic Drice Train Loads in a Wind Turbine. Phd
Thesis. Arenbergkasteel: Katholieke Universiteit Leuven, 2006.
[43] K. G. Pierce. Wind Turbine Load Prediction Using the Beddoes-Leishman Model
for Unsteady Aerodynamics and Dynamic Stall. M.Sc. Thesis. Salt Lake City:
University of Utah, 1996.
Bibliography 92
[44] D. M. Pitt and D. A. Peters. Rotor Dynamic Inow Derivatives and Time Con-
stants from Various Inow Models. In: Conference Proceedings. Ed. by 9th Eu-
ropean Rotorcraft Forum. Stresa 1983.
[45] I. Rechenberg. Evolutionsstrategie 94. Werkstatt Bionik und Evolutionstechnik.
Stuttgart: frommann-holzboog, 1994.
[46] J. G. Schepers and H. Snel. Final Results of the EU JOULE Projects Dynamic
Inow. In: Conference Proceedings. Ed. by ASME Energy Week Wind Confer-
ence. Vennice 1995.
[47] J. Schepers and J. Heijdra. Verication of European Wind Turbine Design Codes.
VEWTDC Final report. Petten: Energy Research Center of the Netherlands,
2002.
[48] D. Simms et al. NREL Unsteady Aerodynamics Experiment in the NASA-Ames
Wind Tunnel: A Comparison of Predictions to Measurements. Technical Report.
Golden, Colorado: National Renewable Energy Laboratory, 2001.
[49] J. L. Tangler. The Nebulous Art of UsingWind-Tunnel Airfoil Data for Predicting
Rotor Performance. Golden, Colorado: National Renewable Energy Laboratory,
2002.
[50] J. L. Tangler. Wind Turbine Post-Stall Airfoil Performance Characteristics
Guidelines for Blade-Element Momentum Methods. Golden, Colorado: National
Renewable Energy Laboratory, 2004.
[51] Tembra GmbH. Internal database of blade masses over radius. 2010.
[52] Tembra GmbH. Strukturelle und aerodynamische Auslegung eines Rotorblattes
mit Active Flow Control (AFC) Elementen in modularer Hybridbauweise. 2010.
[53] T. Theodorsen. General Theroy of aerodynamic instability and the mechanism of
utter. NACA Report No. 496. 1935.
[54] P. S. Veers. Three-Dimensional Wind Simulation. Sandia Report SAND88â0152.
Albuquerque: Sandia National Laboratories, 1988.
[55] L. A. Viterna and R. Corrigan. Fixed pitch rotor performance of large HAWTs.
Brook Park, Ohio: NASA Lewis Research Center, 1981.
[56] T. Williams and C. Kelley. gnuplot 4.4 An Interactive Plotting Program. 2010.
Bibliography 93
[57] D. G. Wilson et al. Optimized Active Aerodynamic Blade Control for Load
Alleviation on Large Wind Turbines. In: Conference Proceedings. Ed. by AWEA
WINDPOWER 2008 Conference & Exhibition. Houston 2008.
A Appendix
A.1 QBladeAE input les for YawDynAE
1 AeroDyn . i p t f i l e from QbladeAE . Sim : FINAL . Wing : AD1 . 5 .
2 SI Units f o r input and output [ SI ]
3 STEADY Dynamic s t a l l model [BEDDOES or STEADY]
4 NO_CM Aerodynamic p i t ch ing moment model [USE_CM ro NO_CM]
5 DYNIN Inf low model [DYNIN or EQUIL ]
6 SWIRL Induct ion f a c t o r model [NONE or WAKE or SWIRL]
7 0 .005 Convergence t o l e r an c e f o r induct ion f a c t o r
8 PRAND Tip−l o s s model (EQUIL only ) [ PRANdtl , GTECH, or NONE]
9 PRAND Hub−l o s s model (EQUIL only ) [ PRANdtl or NONE]
10 wind\smaple_w ind f i l e
11 89 .00 Wind r e f e r e n c e (hub) he ight .
12 0 .10 Tower shadow c e n t e r l i n e v e l o c i t y d e f i c i t .
13 3 .00 Tower shadow ha l f width .
14 4 .00 Tower shadow r e f e r e n c e po int .
15 1 .2250 Air dens i ty .
16 1 .53 e−05 KinVisc − Kinematic a i r v i s c o s i t y
17 0.0010 Time i n t e r v a l f o r aerodynamic c a l c u l a t i o n s .
18 5 Number o f a i r f o i l f i l e s used . F i l e s l i s t e d below :
19 " a i r f o i l s \ c i r c l e _360_Polar . txt "
21 [ l e f t out f o r b r ev i ty ]
23 " a i r f o i l s \DU−96−W−180. txt "
24 25 Number o f blade elements per blade
25 RELM Twist DR CHOD F i l e ID Elem Data
27 0 .75 0 .00 1 .50 2 .19 1 PRINT
28 2 .50 0 .00 2 .00 2 .34 1 PRINT
30 [ l e f t out f o r b r ev i ty ]
32 42 .40 0 .00 0 .80 0 .68 9 PRINT
33 42.925 0 .00 0 .25 0 .41 9 PRINT
34 43.15 0 .00 0 .20 0 .17 9 PRINT
35 SINGLE
Listing A.1: Example aerodyn.ipt input le for QBladeAE simulation
A Appendix 95
1 NACA 63(3)−218, 360 Polar
2 QBladeAE , Reyn#: 3 .0 e6 Mach#: 0 .1
3 3 Number o f a i r f o i l t ab l e s in t h i s f i l e
4 −10 0 10 Table ID parameter
5 0 .00 No longe r used
6 0 .00 No longe r used
7 0 .00 No longe r used
8 0 .00 No longe r used
9 3 .59 −1.50 −6.46 Zero l i f t ang le o f attack ( deg )
10 6 .35 7 .03 6 .57 Cn s l ope f o r zero l i f t ( d imens i on l e s s )
11 0 .70 −0.55 1 .09 Cn at s t a l l va lue f o r p o s i t i v e angle o f attack
12 −1.02 −0.89 −0.60 Cn at s t a l l va lue f o r negat ive angle o f attack
13 0 .00 0 .00 0 .00 Angle o f attack f o r minimum CD ( deg )
14 0 .008 0 .005 0 .008 Minimum CD value
15 −180.0 −0.04 0 .006 −0.04 0 .001 −0.04 0 .006
16 −177.0 0 .25 0 .010 0 .25 0 .012 0 .25 0 .010
18 [ l e f t out f o r b r ev i ty ]
20 6 .5 0 .32 0 .009 0 .93 0 .001 1 .29 0 .016
21 7 .0 0 .37 0 .010 0 .98 0 .012 1 .33 0 .017
22 7 .5 0 .43 0 .010 1 .03 0 .012 1 .36 0 .018
24 [ l e f t out f o r b r ev i ty ]
26 179 .0 −0.24 0 .008 −0.24 0 .001 −0.24 0 .008
27 180 .0 −0.04 0 .006 −0.04 0 .001 −0.04 0 .006
Listing A.2: Example airfoil.ipt input le for QBladeAE simulation
1 YawDyn . i p t f i l e from QbladeAE .
2 90 .00 Time durat ion o f the s imula t i on ( sec )
3 1500 Number o f azimuth s e c t o r s used f o r i n t e g r a t i o n
4 5 .00 Decimation f a c t o r f o r output p r i n t i ng
5 0.0100 TOLER, Trim so l u t i on t o l e r an c e ( deg )
6 3 Number o f b lades
7 0 .00 0 .00 0 .00 I n i t i a l p i t ch ang l e s ( deg )
8 −4.00 Rotor hub s l i n g ( d i s t ance from yaw ax i s to hub ; p o s i t i v e downwind ) (m)
9 −4.00 Shaft t i l t ang le ( deg )
10 −2.00 Rotor precone angle ( deg )
11 15 .00 RPM, ro to r speed in r e vo l u t i on s per minute
12 0 .00 Ps i In i t , I n i t i a l r o to r po s i t i o n ( zero f o r Blade 1 down) ( deg )
13 FIXED Yaw Model : FREE or FIXED yaw system
14 2 .00 I n i t i a l yaw angle ( deg )
15 0 .00 I n i t i a l yaw rate ( deg/ sec )
16 0 .00 Mass moment o f i n e r t i a about yaw ax i s ( kg m^2)
17 0 .00 YawSti fstrong += f , s t i f f n e s s o f yaw spr ing (Nm/rad )
18 0 .00 YawDamp, yaw damping c o e f f i c i e n t (Nm sec )
19 0 .00 YawFriction , constant f r i c t i o n moment at yaw ax i s (Nm)
20 HINGE Hub model : HINGE, TEETER or RIGID
21 0 .00 0 .00 0 .00 I n i t i a l f l a p ang l e s ( deg )
22 0 .00 0 .00 0 .00 I n i t i a l f l a p r a t e s ( deg/ sec )
23 1 .2 RHinge , rad ius o f r o to r hub (m)
24 10 .57 RBar , d i s t ance from hinge to blade c . g . (m)
25 9845.00 Mass o f one blade ( kg )
26 1 .89 e6 Mass moment o f i n e r t i a o f blade about hinge ax i s ( kg m^2)
27 2 .20 e7 Tors iona l s t i f f n e s s o f blade root spr ing (Nm/rad )
28 0 .00 Teeter s l i n g d i s t ance o f t e e t e r ax i s upwind o f r o to r apex (m)
29 0 .00 Free t e e t e r ang le ( deg )
30 0 .00 Teeter s t i f f n e s s , f i r s t or l i n e a r c o e f f . (Nm/rad )
31 0 .00 Teeter s t i f f n e s s , c o e f f . o f d e f l e c t i o n (Nm/rad^2)
32 0 .00 Teeter damping c o e f f i c i e n t ( lb f−f t−sec )
33 1 ,2 ,3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 ,11 ,12 ,13 ,14 ,15 ,16 ,17 [ l e f t out f o r b r ev i ty ]
Listing A.3: Example yawdyn.ipt input le for QBladeAE simulation
A Appendix 96
1 TurbSim v1.50 Input F i l e . Generated in QBladeAE on Do Dez 16 2010 1 2 : 1 2 : 2 1 . . Wing : S88V3 .
3 −−−−−−−−−Runtime Options−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−4 123456 RandSeed1 − F i r s t random seed (−2147483648 to 2147483647)
5 RANLUX RandSeed2 − Second random seed f o r i n t r i n s i c pRNG, or an a l t e r n a t i v e
6 True WrBHHTP − Output hub−he ight turbu lence parameters in binary form?
7 False WrFHHTP − Output hub−he ight turbu lence parameters in formatted form?
8 False WrADHH − Output hub−he ight time−s e r i e s data in AeroDyn form?
9 False WrADFF − Output f u l l−f i e l d time−s e r i e s data in TurbSim/AeroDyn form?
10 True WrBLFF − Output f u l l−f i e l d time−s e r i e s data in BLADED/AeroDyn form?
11 False WrADTWR − Output tower time−s e r i e s data ?
12 False WrFMTFF − Output f u l l−f i e l d time−s e r i e s data in formatted form?
13 False WrACT − Output coherent turbu lence time s t ep s in AeroDyn form?
14 True Clockwise − Clockwise r o t a t i on look ing downwind?
15 0 ScaleIEC − Sca l e IEC turbu lence models to exact t a r g e t standard dev i a t i on ?
17 −−−−−−−−Turbine/Model Sp e c i f i c a t i o n s−−−−−−−−−−−−−−−−−−−−−−−18 25 NumGrid_Z − Ver t i c a l gr id−point matrix dimension
19 25 NumGrid_Y − Hor i zonta l gr id−point matrix dimension
20 0 .05 TimeStep − Time step [ seconds ]
21 600 .0 AnalysisTime − Length o f ana l y s i s time s e r i e s [ seconds ]
22 90 .0 UsableTime − Usable l ength o f output time s e r i e s [ seconds ]
23 89 .0 HubHt − Hub he ight [m] ( should be > 0.5∗GridHeight )
24 107 .0 GridHeight − Grid he ight [m]
25 107 .0 GridWidth − Grid width [m]
26 0 .00 VFlowAng − Ver t i c a l mean f low ( u p t i l t ) ang le [ degree s ]
27 0 .00 HFlowAng − Hor i zonta l mean f low ( skew ) angle [ degree s ]
29 −−−−−−−−Meteo ro l og i ca l Boundary Condit ions−−−−−−−−−−−−−−−−−−−30 IECKAI TurbModel − Turbulence model
31 1 IECstandard − Number o f IEC 61400−x standard
32 A IECturbc − IEC turbu lence c h a r a c t e r i s t i c
33 NTM IEC_WindType − IEC turbu lence type
34 de f au l t ETMc − IEC Extreme Turbulence Model "c" parameter [m/ s ]
35 d e f au l t WindProf i le − Wind p r o f i l e type
36 89 .00 RefHt − Height o f the r e f e r e n c e wind speed [m]
37 13 .00 URef − Mean ( t o t a l ) wind speed at the r e f e r e n c e he ight [m/ s ]
38 d e f au l t ZJetMax − Jet he ight [m]
39 0 .20 PLExp − Power law exponent [− ]
40 0 .03 Z0 − Sur face roughness l ength [m]
42 −−−−−−−−Non−IEC Meteo ro l og i ca l Boundary Condit ions−−−−−−−−−−−−43 de f au l t Lat i tude − S i t e l a t i t u d e [ degree s ] ( or " d e f au l t ")
44 0.0000 RICH_NO − Gradient Richardson number
45 de f au l t UStar − Fr i c t i on or shear v e l o c i t y [m/ s ] ( or " d e f au l t ")
46 de f au l t ZI − Mixing l ay e r depth [m] ( or " d e f au l t ")
47 de f au l t PC_UW − Hub mean u 'w' Reynolds s t r e s s ( or " d e f au l t ")
48 de f au l t PC_UV − Hub mean u ' v ' Reynolds s t r e s s ( or " d e f au l t ")
49 de f au l t PC_VW − Hub mean v 'w' Reynolds s t r e s s ( or " d e f au l t ")
50 de f au l t IncDec1 − u−component coherence parameters
51 de f au l t IncDec2 − v−component coherence parameters
52 de f au l t IncDec3 − w−component coherence parameters
53 de f au l t CohExp − Coherence exponent ( or " d e f au l t ")
55 −−−−−−−−Coherent Turbulence Sca l i ng Parameters−−−−−−−−−−−−−−−−−−−56 . . . CTEventPath − Name o f the path where event data f i l e s are l o ca t ed
57 LES CTEventFile − Type o f event f i l e s ("LES" , "DNS" , or "RANDOM")
58 true Randomize − Randomize the d i s turbance s c a l e and l o c a t i o n s ? ( t rue / f a l s e )
59 1 .00 Di s tSc l − Disturbance s c a l e ( r a t i o o f wave he ight to ro to r d i sk ) .
60 0 .50 CTLy − Frac t i ona l l o c a t i o n o f tower c e n t e r l i n e from r i gh t Randomize = true . )
61 0 .50 CTLz − Frac t i ona l l o c a t i o n o f hub he ight from the bottom of the datase t .
62 30 .0 CTStartTime − Minimum s t a r t time f o r coherent s t r u c t u r e s in RootName . c t s
64 ==================================================
65 NOTE: Do not add or remove any l i n e s in t h i s f i l e !
66 ==================================================
Listing A.4: Example full eld input le for TurbSim
A Appendix 97
1 Active Element . i p t f i l e from QBladeAE .
2 ON Act ivate a c t i v e s imula t i on
3 PID Def ine c on t r o l e r [ PID or LOOP]
4 25 Number o f blade elements per blade
5 1 Number o f a c t i v e e lements per blade
6 RELM AEID
7 1 0
8 2 0
10 [ l e f t out f o r b r ev i ty ]
12 23 1
13 24 0
14 25 2
15 General AE Se t t i ng s
16 1 Step s i z e f o r c on t r o l loop
17 0 Sensor delay
18 AE#1 Set t ing
19 −10 10 Min/Max range o f actuator
20 30 Maximum actuator speed
21 FlapRate Control parameter [LOOP: DFN, DFT; PID : AFMB[kNm] , FlapAngle , FlapRate ]
22 0 Set po int
23 23 Sensor po s i t i o n
24 20 12 5 Kp, Ki , Kd
Listing A.5: Example active.ipt input le for QBladeAE simulation
A Appendix 98
A.2 Geometric blade design
Table A.1: Blade geometric parameters in NREL format used for the simulation.
Radius r[m] Twist [deg] DR [m] Chord c[m] Airfoil
0.75 12.00 1.50 2.19 Cylinder2.50 12.00 2.00 2.34 Cylinder5.00 12.00 3.00 2.79 Morph17.25 11.50 1.50 3.15 Morph29.50 9.97 3.00 3.17 DU-00-W-40112.50 8.06 3.00 3.05 DU-00-W-35015.50 6.59 3.00 2.87 DU-97-W-30018.50 5.32 3.00 2.65 DU-91-W-25021.50 4.18 3.00 2.40 DU-93-W-21024.50 3.17 3.00 2.18 DU-93-W-21027.50 2.34 3.00 1.97 DU-93-W-21029.75 1.70 1.50 1.81 DU-96-W-18031.25 1.41 1.50 1.70 DU-96-W-18032.75 1.16 1.50 1.58 DU-96-W-18034.25 0.91 1.50 1.47 DU-96-W-18035.75 0.69 1.50 1.36 DU-96-W-18037.25 0.47 1.50 1.27 DU-96-W-18038.75 0.25 1.50 1.16 DU-96-W-18040.25 0.04 1.50 1.05 DU-96-W-18041.75 0.02 1.50 0.86 DU-96-W-18042.65 0.00 0.30 0.62 DU-96-W-18042.92 0.00 0.25 0.41 DU-96-W-18043.15 0.00 0.20 0.18 DU-96-W-180