Fluid Mechanics of Vertical Axis Turbines

poisonmammeringΜηχανική

24 Οκτ 2013 (πριν από 3 χρόνια και 11 μήνες)

269 εμφανίσεις

ACTA
UNIVERSITATIS
UPSALIENSIS
UPPSALA
2012
Digital Comprehensive Summaries of Uppsala Dissertations
from the Faculty of Science and Technology 998
Fluid Mechanics of Vertical Axis
Turbines
Simulations and Model Development
ANDERS GOUDE
ISSN 1651-6214
ISBN 978-91-554-8539-9
urn:nbn:se:uu:diva-183794
Dissertation presented at Uppsala University to be publicly examined in Polhemssalen,
Ångströmslaboratoriet, Lägerhyddsvägen 1, Uppsala, Friday, December 14, 2012 at 13:15 for
the degree of Doctor of Philosophy. The examination will be conducted in English.
Abstract
Goude, A. 2012. Fluid Mechanics of Vertical Axis Turbines: Simulations and Model
Development. Acta Universitatis Upsaliensis. Digital Comprehensive Summaries of
Uppsala Dissertations from the Faculty of Science and Technology 998. 111 pp. Uppsala.
ISBN 978-91-554-8539-9.
Two computationally fast fluid mechanical models for vertical axis turbines are the streamtube
and the vortex model. The streamtube model is the fastest, allowing three-dimensional modeling
of the turbine, but lacks a proper time-dependent description of the flow through the turbine. The
vortex model used is two-dimensional, but gives a more complete time-dependent description of
the flow. Effects of a velocity profile and the inclusion of struts have been investigated with the
streamtube model. Simulations with an inhomogeneous velocity profile predict that the power
coefficient of a vertical axis turbine is relatively insensitive to the velocity profile. For the
struts, structural mechanic loads have been computed and the calculations show that if turbines
are designed for high flow velocities, additional struts are required, reducing the efficiency
for lower flow velocities.Turbines in channels and turbine arrays have been studied with the
vortex model. The channel study shows that smaller channels give higher power coefficients
and convergence is obtained in fewer time steps. Simulations on a turbine array were performed
on five turbines in a row and in a zigzag configuration, where better performance is predicted
for the row configuration. The row configuration was extended to ten turbines and it has been
shown that the turbine spacing needs to be increased if the misalignment in flow direction is
large.A control system for the turbine with only the rotational velocity as input has been studied
using the vortex model coupled with an electrical model. According to simulations, this system
can obtain power coefficients close to the theoretical peak values. This control system study
has been extended to a turbine farm. Individual control of each turbine has been compared to a
less costly control system where all turbines are connected to a mutual DC bus through passive
rectifiers. The individual control performs best for aerodynamically independent turbines, but
for aerodynamically coupled turbines, the results show that a mutual DC bus can be a viable
option.Finally, an implementation of the fast multipole method has been made on a graphics
processing unit (GPU) and the performance gain from this platform is demonstrated.
Keywords: Wind power, Marine current power, Vertical axis turbine, Wind farm, Channel
flow, Simulations, Vortex model, Streamtube model, Control system, Graphics processing
unit, CUDA, Fast multipole method
Anders Goude, Uppsala University, Department of Engineering Sciences, Electricity,
Box 534, SE-751 21 Uppsala, Sweden.
© Anders Goude 2012
ISSN 1651-6214
ISBN 978-91-554-8539-9
urn:nbn:se:uu:diva-183794 (http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-183794)
To my family
List of papers
This thesis is based on the following papers,which are referred to in the text
by their Roman numerals.
I Goude,A.,Lundin,S.,Leijon,M.,“A parameter study of the influence
of struts on the performance of a vertical-axis marine current turbine”,
In “Proceedings of the 8th European wave and tidal energy conference,
EWTEC2009”,Uppsala,Sweden,pp.477–483,September 2009.
II Goude,A.,Lalander,E.,Leijon,M.,“Influence of a varying vertical
velocity profile on turbine efficiency for a Vertical Axis Marine Current
Turbine”,In “Proceedings of the 28th International Conference on
Offshore Mechanics and Arctic Engineering,OMAE 2009”,Honolulu,
USA,May 2009.
III Grabbe,M.,Yuen K.,Goude,A.,Lalander,E.,Leijon,M.,“Design of
an experimental setup for hydro-kinetic energy conversion”,
International Journal on Hydropower &Dams,15(5),pp.112–116,
2009.
IV Goude,A.,Ågren,O.,“Simulations of a vertical axis turbine in a
channel”,Submitted to Renewable Energy,October 2012.
V Goude,A.,Ågren,O.,“Numerical simulation of a farmof vertical axis
marine current turbines”,In “Proceedings of the 29th International
Conference on Offshore Mechanics and Arctic Engineering,OMAE
2010”,Shanghai,China,June 2010.
VI Dyachuk,E.,Goude,A.,Lalander,E.,Bernhoff,H.,“Influence of
incoming flow direction on spacing between vertical axis marine
current turbines placed in a row”,In “Proceedings of the 31th
International Conference on Offshore Mechanics and Arctic
Engineering,OMAE 2012”,Rio de Janeiro,Brazil,July 2012.
VII Goude,A.,Bülow,F.,“Robust VAWT control systemevaluation by
coupled aerodynamic and electrical simulation”,Submitted to
Renewable Energy,September 2012.
VIII Goude,A.,Bülow,F.,“Aerodynamic and electric evaluation of a
VAWT farmcontrol systemwith passive rectifiers and mutual
DC-bus”,Submitted to Renewable Energy,November 2012.
IX Goude,A.,Engblom,S.,“Adaptive fast multipole methods on the
GPU”,Journal of Supercomputing,DOI 10.1007/s11227-012-0836-0,
In Press,October 2012.
Reprints were made with permission fromthe publishers.
The author has also contributed to the following paper,not included in the
thesis:
A Yuen,K.,Lundin,S.,Grabbe,M.,Lalander,E.,Goude,A.,Leijon,
M.,“The Söderfors Project:Construction of an Experimental Hydroki-
netic Power Station”,In “Proceedings of the 9th European wave and
tidal energy conference,EWTEC2011”,Southampton,United Kingdom,
September 2011.
Contents
1 Introduction
................................................................................................
14
1.1 Different turbine types
...................................................................
14
1.2 Comparison between wind and marine current turbines
.............
16
1.3 Vertical axis turbine research at Uppsala University
...................
17
1.4 Extended studies within this thesis
...............................................
18
1.5 Outline of the thesis
.......................................................................
18
2 Theory for vertical axis turbines
...............................................................
20
2.1 Basic theory and the Betz limit
.....................................................
20
2.2 Extension to include channels
.......................................................
22
2.3 Theory of lift-based vertical axis turbines
....................................
24
2.4 Angle of attack including flow curvature
.....................................
26
3 Control strategy for vertical axis turbines
................................................
30
3.1 Control of a single turbine
.............................................................
30
3.2 Extension to multiple turbines
.......................................................
36
4 Simulation models
.....................................................................................
37
4.1 Streamtube models
.........................................................................
37
4.1.1 Description of model
.......................................................
39
4.1.2 Including struts
................................................................
41
4.1.3 Obtaining lift and drag coefficients
................................
43
4.1.4 Corrections due to flow curvature
..................................
45
4.1.5 Including flow expansion
................................................
46
4.2 Vortex models
.................................................................................
48
4.2.1 Implementing the turbine
................................................
49
4.2.2 Merging vortices
..............................................................
51
4.2.3 Calculation of velocity
....................................................
51
4.2.4 Numerical evaluation of the velocity field
.....................
53
5 Simulation results
......................................................................................
66
5.1 Evaluation of simulation tools
.......................................................
66
5.1.1 Strut modeling
.................................................................
68
5.1.2 Expansion model
.............................................................
69
5.1.3 Tip correction model
.......................................................
70
5.1.4 Curvature modeling
.........................................................
70
5.1.5 Vortex model
....................................................................
73
5.1.6 Concluding remarks about the simulation tools
............
74
5.2 Results from papers
........................................................................
74
5.2.1 The effects of struts
.........................................................
75
5.2.2 The effects of a velocity profile
......................................
77
5.2.3 Design of a turbine for use in a river
..............................
79
5.2.4 Turbines in channels
........................................................
80
5.2.5 Turbines in an array
.........................................................
83
5.2.6 Simulations of control systems
.......................................
87
5.2.7 Control of multiple turbines
............................................
90
6 Conclusions
................................................................................................
95
7 Suggestions for future work
......................................................................
97
8 Summary of papers
....................................................................................
98
9 Errata for papers
......................................................................................
102
10 Acknowledgments
...................................................................................
103
11 Summary in Swedish
...............................................................................
104
References
......................................................................................................
107
Nomenclature
A
m
2
Turbine cross-sectional area
A

m
2
Asymptotic area of streamtube enclosing turbine
A
c
m
2
Cross-sectional area of a channel
A
d
m
2
Area of turbine disc
A
e
m
2
Streamtube area far downstream/center of turbine
AR

Aspect ratio of a blade
C
D

Drag coefficient
C
Ds

Drag coefficient for strut
C
D∞

Drag coefficient for infinitely long blade
C
L

Lift coefficient
C
Ls

Lift coefficient for strut
C
L∞

Lift coefficient for infinitely long blade
C
N

Normal force coefficient
C
N0

Normal force coefficient without curvature corrections
C
Ns

Normal force coefficient for strut
C
P

Power coefficient
C
P
e

Power coefficient equivalent for extracted power
C
Pmax

Maximumpower coefficient for a given flowvelocity
C
T

Tangential force coefficient
C
Ts

Tangential force coefficient for strut
D
m
Turbine diameter
F

Velocity correction factor
F
D
N
Drag force
F
D0
N
Drag force at zero angle of attack
F
L
N
Lift force
F
N
N
Normal force
F
Nl
N/m
Normal force per meter
F
R
N
Force in radial direction
F
T
N
Tangential force
F
Ts
N
Tangential force on strut
F
x
N
Aerodynamic force from a blade on flow in a streamtube
F
xs
N
Aerodynamic force from a strut on flowin a streamtube
H
m
Channel height
J
kgm
2
Moment of intertia
L
m
Distance between struts
K
N
Constant to determine lift force
9
Ma

Mach number
N

Number of particles (FMM)
N
b

Number of blades
N
box

Number of boxes (FMM)
N
d

Number particles per box (FMM)
N
p

Number of panels
N
s

Number of struts
N
t

Number of turbines
N
v

Number of vortices
NF
i

Set of all boxes in near field of box i (FMM)
P
W
Power absorbed by the turbine
P
e
W
Power extracted fromthe turbine
P
e,tot
W
Total power extracted fromthe turbines of a farm
P
tot
W
Power available in flow
R
m
Turbine radius
R
inner
m
Strut inner attachment point
R
outer
m
Strut outer attachment point
Re

Reynolds number
T
s
Nm
Torque fromstrut
V
m/s
Flow velocity
V
0
m/s
Flow velocity far upstreamof turbine
V

m/s
Asymptotic flowvelocity
V
abs
m/s
Magnitude of incoming flow velocity at blade position
V
b
m/s
Blade velocity
V
d
m/s
Flow velocity at turbine disc
V
e
m/s
Flow velocity far downstream/center of turbine

V
i
m/s
Vortex velocity
V
r
m/s
Relative flow velocity for a blade (absolute value)
V
rs
m/s
Relative flow velocity for a strut (absolute value)
V
s
m/s
Flow velocity at strut position
V
re f
m/s
Reference flow velocity for estimating angle of attack
V
rel
m/s
Relative flow velocity for blade
V
relz
m/s
Relative flow velocity for blade in its own reference frame
V
s
m/s
Far downstreamvelocity of flowpassing outside turbine
V
s j
m/s
Flow velocity at strut segment
V
ω
m/s
Velocity due to vortices
W
m
2
/s
Complex velocity potential
W
b
m
2
/s
Complex velocity potential for blade velocity
a

Axial induction factor
a
i

Multipole coefficient i (FMM)
a
s

Slope of lift coefficient curve in C
L
/α plot
b
m
Circle radius used for conformal mapping
b
i

Local coefficient i (FMM)
10
b
y∞

Normalized asymptotic streamtube width
b
ye

Normalized streamtube width at the turbine center
b
z∞

Normalized asymptotic streamtube height
b
ze

Normalized streamtube height at the turbine center
c
m
Blade chord
c
s
m
Strut chord
c
sound
m/s
Speed of sound
c
0
m
Reference chord for struts
g
m/s
2
Gravitational acceleration
h
m
Turbine height
k
1
kgm
2
First control systemconstant
k
2
kgm
2
Second control systemconstant
k
3
kgm
2
Third control systemconstant
k
d1

Constant for time estimate of direct evaluation (FMM)
k
d2

Constant for time estimate of direct evaluation (FMM)
l
m
Blade length in streamtube
p

Number of multipole coefficients (FMM)
p
0
N/m
2
Pressure far upstreamof turbine
p
atm
N/m
2
Atmospheric pressure
p
d1
N/m
2
Pressure directly in front of turbine disc
p
d2
N/m
2
Pressure directly after turbine disc
p
e
N/m
2
Pressure far downstreamof turbine
r
0
m
Box center (FMM)
r
s
m
Radial position on a strut
r
m
Arbitrary position
r
i
m
Vortex position
s
m
Position on blade surface in transformed plane
t
s
Time
t
b
m
Blade thickness
t
d
s
Time estimate for direct evaluation (FMM)
u

Interference factor
x
m
Position in the x-direction
x
0
m
Blade attachment point
x
0r

Normalized blade attachment point
y
m
Position in the y-direction
y

m
Asymptotic streamtube position in y-direction
y
d
m
Streamtube position at turbine disc in y-direction
y
e
m
Streamtube position far downstreamin y-direction
Δy
m
Streamtube width
z
m
Position on blade surface in the blades reference frame
z
0
m
Position on blade surface in the turbines reference frame
z
b
m
Blade position
Δz
m
Streamtube height
11

Γ
p
m
3
/s
Circulation of three dimensional point vortex
Γ
m
2
/s
Circulation of two dimensional vortex
Δ
m
Cutoff radius used for vortex merging
Ω
rad/s
Turbine rotational velocity
Ω
i
rad/s
Rotational velocity of turbine i
Ω
1
rad/s
First control systemrotational velocity constant
Ω
2
rad/s
Second control systemrotational velocity constant
Ω
3
rad/s
Third control systemrotational velocity constant
α

Angle of attack
α
b

Corrected angle of attack
α
s

Angle of attack for strut
β

Direction of incoming wind
δ

Blade pitch angle
ε
m
Cutoff radius of Gaussian vortex kernel
η

Angle of blade relative to the vertical axis
η
s

Angle of strut relative to the horizontal plane
θ

Blade azimuthal position shifted 90 degrees
θ
b

Blade azimuthal position
λ

Tip speed ratio
λ
e

Equilibriumtip speed ratio
λ
max

The tip speed ratio that gives highest power coefficient
ν
m
2
/s
Kinematic viscosity
ρ
kg/m
3
Density of fluid
σ
N/m
2
Stress
ϕ

Angle of relative wind
φ
m
2
/s
Velocity potential
ω
1/s
Vorticity
12
Abbreviations
CPU
Central processing unit
CUDA
Compute unified device architecture
FEM
Finite element method
FMM
Fast multipole method
FVM
Finite volume method
GPU
Graphics processing unit
L2L
Local to local translation
L2P
Local to particle evaluation
M2L
Multipole to local translation
M2M
Multipole to multipole translation
NACA
National Advisory Committee for Aeronautics
P2M
Particle to multipole initialization
P2P
Particle to particle interaction (direct evaluation)
RANS
Reynolds-averaged Navier-Stokes equations
SIMD
Single instruction multiple data
SSE
Streaming SIMD Extensions
13
1.Introduction
Aturbine is used to convert the energy froma moving fluid into rotational mo-
tion,which in turn can drive an electric generator.The best suited turbine for
this energy conversion depends on the characteristics of the flow.One case,
which is present in gas turbines and traditional hydro power plants,is flowcon-
strained by walls,such as flowwithin a pipe.Here,the turbine cross-sectional
area usually covers the entire flowand the flowhas a large pressure difference
that drives the turbine.Asecond case is the free flow,where no confining walls
are present.This is the case for wind power and often a reasonable approxi-
mation for tidal power in the ocean.In free flow,the fluid can pass around the
turbine and the available energy is the kinetic energy in the flow.This thesis
mainly treat the free flow case,but also a hybrid case where there are confin-
ing walls,but the turbine does not cover the whole cross-sectional area of the
flow,is studied in this thesis.This situation occurs in a river,where the river
cross-sectional shape area usually prevents a turbine from covering the entire
cross-section.
1.1 Different turbine types
The typical turbine design for wind power is a horizontal axis turbine,where
the rotational axis of the turbine is parallel to the flowdirection [1,chapter 1].
In this thesis,however,the vertical axis turbine will be investigated.Here,the
rotational axis is perpendicular to the flow direction.This kind of turbine is
sometimes called “cross-flow turbine”,as the turbine in principle also can be
tilted 90 degrees to have a horizontal axis while still having its rotational axis
perpendicular to the flow.The traditional name “vertical axis turbine” will be
used here,even for situations where the rotational axis is tilted,since this is
the most commonly used name.
There are two different types of vertical axis turbines.The first type is
based on the drag force and is often called the Savonius rotor after the Finnish
inventor Sigurd Johannes Savonius,despite that Savonius only patented an im-
provement of older designs [2].This improvement is neither implemented on
all present drag-based turbines.Drag-based devices rely on variation of the
drag coefficient with respect to the orientation of the object.To create a rea-
sonably efficient drag-based turbine,the drag coefficient should be high in one
direction and low in the opposite direction,which gives a torque on the tur-
bine.Drag-based devices achieve lower power coefficients than the lift-based
14
Figure 1.1.Different types of vertical axis turbines.
devices described in the section below[3,chapters 2,7].Another drawback is
that the amount of construction material in drag devices is quite high (as can
be seen in figure 1.1).This cost is inhibiting the construction of large turbines,
as material usage is proportional to the volume,i.e.the cube of the character-
istic length of the turbine,while the power absorption is proportional to the
cross-sectional area,i.e.the square of the characteristic length.
The second type is the lift-based turbine,which was originally invented by
the French engineer George Jean Marie Darrieus [4] in the 1920’s (approxi-
mately one year after Savonius patented his design).The patent application of
Darrieus covers both the curved blade turbine and the H-rotor (see figure 1.1),
as well as turbines with varying pitch angle and ducted turbines.It is suggested
in the patent that the designs work both for wind and tidal energy.The aim
of the curved blade design is to reduce the bending stresses in the blades due
to centrifugal forces.The North American company Flowind commercialized
in the 1980s the Darrieus turbine with the curved blade design [5,chapter 1].
During that time,the curved blade turbine was also studied by Sandia National
Laboratories,which is the main reason why much of the published work on
Darrieus turbines is on the curved blade design.This thesis will instead focus
on the straight blade H-rotor design,which currently is in development at Up-
psala University.With recent progress for light materials,composites can be
used in the turbine construction,which reduces centrifugal forces due to the
15
lighter structure.This makes the H-rotor design more feasible.The straight
blade design has the advantages that straight blades are easier to manufacture
and by attaching the blades with struts,it is possible to place the upper bear-
ing much closer to the turbine center,reducing the bending moment on the
axis.In addition,the constant radius of the straight blade design gives a larger
cross-sectional area.Disadvantages compared to the curved blade design are
the addition of extra struts and the higher bending moments due to centrifugal
forces.
The main aerodynamic advantage of vertical axis turbines,compared to
standard horizontal axis turbines,is the independence of flow direction,re-
moving the need for a yaw mechanism.For water flow,an additional advan-
tage is that the cross-sectional area can be more flexibly chosen as both height
and diameter can be varied (and the diameter can vary with height).This can
be useful in shallowwater where a turbine with a large width and small height
can cover a larger area than a horizontal axis turbine,as the cross-sectional
area of a horizontal axis turbine is circular.Disadvantages of the vertical axis
turbines are the lower power coefficients and that the turbines are typically not
self-starting.
The vertical axis turbine can have its generator on the ground,which in
the wind power case simplifies maintenance,tower construction and makes
the weight of the generator less important.This is beneficial for direct driven
generators,which typically have large diameters.The use of direct driven gen-
erators further reduces the number of moving parts in the system.One major
concern for vertical axis turbines is the cyclic blade forces in each revolution,
which leads to torque oscillations and material fatigue.For further compar-
isons between horizontal and vertical axis turbines (and also between curved
and straight blade turbines) see e.g.[6].
1.2 Comparison between wind and marine current
turbines
Even though wind turbines operate in air (gas) while marine current turbines
operate in water (liquid),there are many similarities between the two.Tra-
ditionally,water is considered an incompressible fluid and can therefore be
modeled with the incompressible Navier-Stokes equations.For air,it is typi-
cally expected that compressibility effects can be neglected for Mach numbers
Ma within the range
Ma =
V
rel
c
sound
<0.3.
Here,V
rel
is the relative flow velocity (measured in the blades rest frame) and
c
sound
is the speed of sound in the fluid [7,chapter 9].Note that the major con-
tribution to V
rel
originates fromthe blades’ own motion for lift-based turbines.
16
The speed of the blades in a wind turbine is typically too low for the Mach
number to be above 0.3 and wind turbines can therefore also be modeled with
incompressible aerodynamics.Although both wind and marine turbines can be
studied with the incompressible Navier-Stokes equations,there are still some
characteristic differences.One difference is that for marine current turbines,
there is both a sea bed and a free surface that bounds the flow.Another dif-
ference is the risk for cavitation at too high flow velocities.Cavitation would
modify the flow characteristics and can cause damage to the turbine [8].
The energy absorbed by a turbine is proportional to the fluid density and
the cube of the flow velocity.As the density of water is 800 times higher than
the density of air,comparatively low fluid velocities are adequate for marine
current power generation.For equal cross-sectional area,a wind speed of
10 m/s has the same incoming kinetic power as a water flowspeed of 1.1 m/s.
However,as the forces only are proportional to the square of the flowvelocity,
the marine current turbine experiences approximately 9.3 times higher fluid
mechanical forces than a wind turbine at the same conditions (assuming that
the turbines are identical) and rotates 9.3 times slower.The increased forces
for marine current turbines require both stronger blades and support structure.
Many of the experimental vertical axis turbines for water have used relatively
large blades and thereby low optimal tip speed ratios [9–12].
One important parameter for the effectiveness of the turbine is the Reynolds
number,which for a blade is defined as
Re =
cV
ν
,
where V is the flowvelocity,c is the blade chord and ν is the kinematic viscos-
ity.Ahigher Reynolds number usually decreases the drag losses and increases
the stall angle,which is beneficial for vertical axis turbines.For 20

C,the
kinematic viscosities are 15.1 μm
2
/s for air and 1.00 μm
2
/s for water [13,ap-
pendix A].Under the conditions of equal power extraction mentioned above
(i.e.9.3 times higher flow velocity for the wind turbine),this would give a
63 %higher Reynolds number for the marine current turbine,which is within
the same order of magnitude as the wind turbine.
1.3 Vertical axis turbine research at Uppsala University
At the Division of Electricity at Uppsala University,three vertical axis wind
turbines have been built.The first turbine had a cross-sectional area of 6 m
2
and was later followed by a turbine with the cross-sectional area 30 m
2
and the
rated power 12 kW[14–16].This larger turbine is used for most of the experi-
ments.A10 kWturbine for telecomapplications has also been built [17].Fur-
ther,a 200 kWturbine has been constructed by the spin-off company Vertical
Wind AB [18].Additionally,a marine current turbine (described in paper III)
is scheduled to be deployed by the end of 2012.
17
Several simulation tools for turbine simulations have previously been de-
veloped at the division.A two-dimensional inviscid vortex model based on
conformal mappings for the blades has been created by Deglaire et al.[19].In
the turbine implementation,each blade is solved independently [20],allowing
for coupling to an elastic method developed by Bouquerel et al.,see paper IV
in [21].A multibody version for simulating turbines has been developed by
Österberg et al.,see [22] and paper III in [21].Two streamtube models have
also been implemented by Deglaire and Bouquerel.
1.4 Extended studies within this thesis
This thesis focuses on the fluid mechanical modeling of the vertical axis tur-
bine,and two different simulation tools have been developed.The first simu-
lation tool uses the streamtube model and the development of this tool started
from the basic streamtube model implemented by Bouquerel,which is based
on the model of Paraschivoiu [3].All additional modeling and code develop-
ment have been developed within this thesis.
The second simulation tool uses a vortex model,and this tool has been
developed from scratch within this work.This model is based on empirical
data for lift and drag coefficients instead of the conformal mapping method
by Deglaire,which is based on inviscid theory.The computational speed is
crucial for the developed vortex model and large efforts have been put into
this.The existing implementation of the fast multipole method by Stefan Eng-
blom[23] has been significantly improved and ported to a GPU (paper IX)
Several studies have been carried out with the two simulation models.The
streamtube model has been used to study losses due to struts (paper I),the
effects of a velocity profile (Paper II) and to design a turbine for deployment
in a river (paper III).The more computationally demanding vortex model has
been used to study turbines in channels (paper IV) and turbine arrays (paper V
and VI).The vortex model is also coupled to an electrical model to study con-
trol systems for a single turbine (paper VII) and extended simulations analyze
control systems for a turbine farm (paper VIII).
1.5 Outline of the thesis
After the introduction,theory for vertical axis turbines is presented in chap-
ter 2.This is followed by an introduction to control systems in chapter 3.
The theory and implementation for the simulation models are then presented
in chapter 4,which also includes the GPU implementation of the fast multi-
pole method.The results from the simulations are given in chapter 5,where
the first part evaluates the accuracy of the simulation models and the second
part summarizes the results from the articles.The thesis ends by conclusions,
18
suggestions for future work,summary of papers,errata for papers and ac-
knowledgments.
19
2.Theory for vertical axis turbines
Given a cross-sectional area A perpendicular to a homogenous flow of a fluid,
the kinetic power that passes through this area is given by
P
tot
=
1
2
ρAV
3
,(2.1)
where ρ is the density and V is the flow velocity.If the flow is not confined
by any surrounding boundaries,the kinetic power is the available power for a
wind/current turbine.The efficiency (i.e.outgoing power divided by incoming
power) would be one possible measure of how good the energy conversion
is.However,adding a turbine will change the velocity and force parts of the
flow to pass outside the turbine area and thereby change the kinetic energy
that passes through this area A.Moreover,some kinetic energy is left in the
flow and can possibly be used later.Therefore,turbine performance is usually
measured with the power coefficient instead,which is defined as
C
P
=
P
1
2
ρAV
3

,(2.2)
where P is the power absorbed by the turbine and V

is the asymptotic up-
stream flow velocity.With this expression,the absorbed power is compared
to the power that would have passed through the cross-sectional area,if the
turbine would be absent,instead of compared to the power that actually passes
through the area.Since this expression is normalized against an expression
that does not change with the turbine characteristics,it is a better measure than
efficiency.Improving the power coefficient will give higher power absorption,
which is not always the case with efficiency.
2.1 Basic theory and the Betz limit
One of the most basic approximations of a turbine is the one used in the tradi-
tional Betz theory [24],where the turbine is approximated as a single flat disc
with a constant pressure drop over the whole turbine surface.All flowpassing
through the disc is encapsulated in a streamtube that starts far ahead of the
turbine and ends far behind.By making the assumption that the pressure at
both ends of the streamtube is the atmospheric pressure p
atm
,and by using the
20
p
atm
p
atm
V

A

V
d
p
d1
p
d2
V
e
A
e
A
d
Figure 2.1.The streamtube used in the Betz limit derivation.The dashed line shows
the control volume used for momentumconservation.
Bernoulli equation before and after the turbine
p
atm
+
1
2
ρV
2

= p
d1
+
1
2
ρV
2
d
,(2.3)
p
d2
+
1
2
ρV
2
d
= p
atm
+
1
2
ρV
2
e
,(2.4)
combined with continuity
A

V

=A
d
V
d
=A
e
V
e
(2.5)
and momentumconservation for the control volume in figure 2.1 (marked with
the dashed line)
ρA
e
V
2
e
−ρA

V
2

=A
d
(p
d2
−p
d1
),(2.6)
it is possible to showthat the velocity at the turbine disc V
d
is equal to
V
d
=
V

+V
e
2
.(2.7)
Considering that the Betz theory assumes no losses,the power absorbed by
the turbine is given as the difference between incoming and outgoing power in
the fluid
P =
1
2
ρA

V
3


1
2
ρA
e
V
3
e
.(2.8)
If the axial induction factor a,defined as
V
d
=(1−a)V

,(2.9)
is combined with expression with equations (2.5) and (2.8),it can be shown
that the power will reach its maximum value for a = 1/3,and the optimal
power is given by
P =
16
27
·
1
2
ρA
d
V
3

(2.10)
where 16/27 is the traditional Betz limit,limiting the power coefficient to
approximately 59.3 %.
21
p
0
p
0
V
0
V
0
V
0
A
0
A
c
V
d
p
d1
p
d2
V
e
A
e
A
d
p
e
V
s
Figure 2.2.Illustration of a streamtube confining the flow that passes through the
turbine disc for a channel of cross-sectional area A
c
.
2.2 Extension to include channels
The Betz theory assumes no outer boundaries in the system.For a turbine op-
erating in water,it is more common that boundaries are present.One example
is a river,where the flow is limited by the width and depth.The outer walls
prevent the flowfromexpanding,pushing more flowthrough the turbine.This
occurs in a traditional hydro power plant,where the entire flow is forced to
pass through the turbine,which results in much higher power absorption than
the Betz limit [25].
To analyze this case analytically,assume that the flow upstream of the tur-
bine has constant velocity V
0
and pressure p
0
(see figure 2.2).Note that in
this case,the pressure upstream and the pressure downstream are not equal.
Instead,there will be a drop in pressure,which,for open channel flow,would
correspond to a drop in the surface level.Due to the continuity of the pressure,
the pressure inside the streamtube and outside has to be the same downstream
(p
e
).The cross-sectional area of the channel is A
c
and the cross-sectional area
of the turbine is A
d
.In this case,the Bernoulli equation gives
p
0
+
1
2
ρV
2
0
= p
d1
+
1
2
ρV
2
d
,(2.11)
p
d2
+
1
2
ρV
2
d
= p
e
+
1
2
ρV
2
e
,(2.12)
p
0
+
1
2
ρV
2
0
= p
e
+
1
2
ρV
2
s
,(2.13)
the continuity equation gives
A
0
V
0
=A
d
V
d
=A
e
V
e
,(2.14)
(A
c
−A
0
)V
0
=(A
c
−A
e
)V
s
(2.15)
22
and momentum conservation for a control volume that encloses the entire
channel gives
ρA
e
V
2
e
+ρA
s
V
2
s
−ρAV
2
0
=A
d
(p
d2
−p
d1
) +A
c
(p
0
−p
e
).(2.16)
The velocity at the turbine can be derived from equations (2.11) – (2.16) as
V
d
=
V
e
(V
s
+V
e
)
V
s
+2V
e
−V
0
,(2.17)
which in the free flowlimit (V
s
→V
0
) reduces to equation (2.7).The force on
the turbine becomes
F
x
=A
d
(p
d1
−p
d2
) =
1
2
A
d

V
2
s
−V
2
e

,(2.18)
giving the power as
P =F
x
V
d
=
1
2
A
d
V
e
(V
s
+V
e
)

V
2
s
−V
2
e

V
s
+2V
e
−V
0
.(2.19)
From equation (2.19),it can be found that the highest power absorption is
obtained when V
e
=V
0
/3,which actually is the same as for the free flow.The
highest power is thus given by
P =
16
27
·
1

1−
A
d
A
c

2
·
1
2
ρA
d
V
3
0
(2.20)
and the pressure drop is
p
e
−p
0
=
4
A
d
A
c

3−
A
d
A
c

9

1−
A
d
A
c

2
ρV
2
0
.(2.21)
From these results,it can be seen that the maximum theoretical power co-
efficient increases with the factor (1−A
d
/A
c
)
−2
for a channel.In the limit
A
d
→A
c
,the power coefficient diverges,along with the pressure drop in equa-
tion (2.21).Considering that an infinite drop in pressure is unfeasible,when
the turbine area is almost as large as the channel,the available pressure dif-
ference will start limiting the maximumpower coefficient,which will prevent
infinite energy extraction.This is the case for hydro power turbines,where the
power is limited by the difference in water elevation.
The model above assumes that the cross-section of the channel is constant.
An open channel will have a drop in surface level over the turbine and an
extension of the model to include this drop is given by Whelan et al.[26].
This correction has not been included in the present work,as the model is
used for comparisons with the two-dimensional vortex simulations where no
free surface is modeled.
23
R
θ
b
Ω
x
y
Figure 2.3.Illustration of a vertical axis turbine.
2.3 Theory of lift-based vertical axis turbines
The H-rotor is a lift-based design,which means that the aerodynamic torque is
generated by the lift force of the blades.Therefore,airfoil profiles are typically
used for the blades.Airfoils are generally designed to operate at relatively low
angles of attack,where the lift force increases approximately linearly with
the angle of attack and the drag force remains low.When the angle of attack
increases above the stall angle,which for a typical airfoil occurs for angles
around 10 – 15 degrees,the lift force is reduced,and the drag starts to increase
substantially.This work is focused towards blades with a fixed pitch angle.To
keep the angle of attack lowwithout pitching the blades,the blades must move
with a high velocity if the wind comes fromthe side.To illustrate this,assume
that a blade is located at angle θ
b
(see figure 2.3).Using complex notation,
this gives the blade position
z
b
=Re

b
(2.22)
and the velocity of the blade is therefore
V
b
=i
˙
θ
b
Re

b
.(2.23)
The rotational velocity
˙
θ
b
is commonly denoted Ω.The incoming wind V,if
complex,represents wind fromany direction.The blade will nowsee a relative
wind of
V
rel
=V −i ΩRe

b
.(2.24)
24
V
rel
α
δ
ϕ
V
b
Figure 2.4.Definitions of angles and velocities.The positive direction for angles is
counter-clockwise,hence α and ϕ are negative for the directions of V
b
and V
rel
in the
figure.
To obtain the angle of relative wind,rotate V
rel
with the angle ie
−iθ
b
,aligning
the blade motion with the negative real axis,
V
relz
=Vie
−iθ
b
+ΩR.(2.25)
The angle of relative wind will be the argument of this complex number.In
the special case of V being real,the angle of relative wind is
ϕ=arctan
cosθ
b
ΩR
V
+sinθ
b
(2.26)
and the absolute value of the relative velocity is
|V
rel
| =V


ΩR
V
+sinθ
b

2
+(cosθ
b
)
2
.(2.27)
The angle of attack is given by
α=ϕ+δ,(2.28)
where δ is the blade pitch angle (see figure 2.4).Equation (2.26) shows how
the angle of attack varies during a turbine revolution and how it decreases as
the rotational velocity increases.As an example,ΩR/V =4 gives a maximum
angle of attack of around 14 degrees,approximately where stall begins to oc-
cur.It should be noted that in equations (2.24) – (2.27),V is the flowvelocity
at the blade position.This can be compared with the tip speed ratio,which is
defined as
λ =
ΩR
V

(2.29)
where the asymptotic velocity is used.Due to the energy extracted,the veloc-
ity at the blade will generally be lower than the asymptotic velocity.
25
The turbine torque is calculated from the tangential force,which is given
by
F
T
=F
L
sinϕ−F
D
cosϕ.(2.30)
When the angle of relative wind is low,the approximations sinϕ ≈ ϕ and
cosϕ ≈ 1 can be applied.Assume that the pitch angle is zero,hence α =
ϕ.For symmetric blades,the blade forces can be approximated as F
L
≈Kα,
where K is a constant,and F
D
≈F
D0
,where F
D0
is a constant.The tangential
force can therefore be estimated as
F
T
=Kϕ
2
−F
D0
,(2.31)
showing that when the angle of relative wind decreases,the drag force be-
comes dominating.The conclusion is that for high tip speed ratios,drag will
give a more significant contribution,reducing the power coefficient.At too
low tip speed ratios,the turbine will enter stall,where lift decreases and drag
increases,which also should be avoided.For these reasons,the turbine should
be designed to operate with a tip speed ratio close to the stall limit in order to
obtain the highest possible power coefficient.
2.4 Angle of attack including flow curvature
The expression (2.26) is only valid for infinitely small symmetric blades.The
blade performs a rotational motion,which leads to additional curvature effects,
changing the effective angle of attack.To conclude the theory section,a more
proper derivation will be performed using a rotating flat plate instead.
By the use of conformal mappings,a circle can be transformed into a flat
plate with the Joukowski transformation.The s-plane represents the circle
with radius b and the z-plane a flat plate extending between −2b to 2b giving
the blade chord c as c = 4b.The blade coordinates z in its own frame of
reference is given by
z =s +
b
2
s
.(2.32)
Using the same transformation as Deglaire [20],the z
0
plane can be defined as
z
0
=

(z +x
0
)e
−iδ
+iR


e

=(z +D)e
i(θ−δ)
,(2.33)
where D = x
0
+iRe

.By assuming that the blade only rotates around the
center,the blade velocity is given by
V
b
= i
˙
θ

(z +x
0
)e
−iδ
+iR


e

=i Ωz
0
(2.34)
with Ω=
˙
θ.Introduce a complex velocity potential W,with complex conju-
gate
W,such that
d
W
d
z
0
=V,(2.35)
26
hence the potential can be used to calculate the flow velocity V.Assume that
the potential is given on the form
W(s) =V
abs
e
i(−β+θ−δ)
s +V
abs
e
−i(−β+θ−δ)
b
2
s



log(s) +W
1
(s),(2.36)
where V
abs
e

is the flow velocity at the blade position.Now,construct a
potential such that
d
W
b
d
z
0
= V
b
,(2.37)
which gives
dW
b
ds
=
V
b
dz
0
ds
=−i Ω
z
0
dz
0
ds
=−i Ω

z +
D

dz
ds
.(2.38)
Note that on the boundary,
z =z.Integrated,this is
W
b
= −i Ω

1
2
z
2
+
Dz

= −i Ω

1
2

s
2
+2b
2
+
b
4
s
2

+
D

s +
b
2
s

.(2.39)
The no-penetration boundary condition states that the stream function should
be constant (possibly time-dependent) on the boundary.Given that the bound-
ary is moving,the condition becomes
Im[W(s)−W
b
(s)] =C.(2.40)
The first part of W(s) already fulfills the condition,while W
1
(s) remains to be
determined,hence
Im[W
1
(s) −W
b
(s)] =C.(2.41)
The boundary condition at infinity states
dW
1
ds




s→∞
=0,(2.42)
and on the boundary,
s =
b
2
s
(2.43)
applies.Write equation (2.41) in terms of complex conjugates
iC = W
1
(s)−
W
1
(s)+i Ω

1
2

s
2
+2b
2
+
b
4
s
2

+
1
2

s
2
+2b
2
+
b
4
s
2

+
D

s +
b
2
s

+D

s +
b
2
s

= W
1
(s)−
W
1
(s)+
i Ω

b
4
s
2
+
b
4
s
2
+
b
2
s
b+
b
2
s
D+
b
2
s
D+
D
b
2
s
+2b
2

.(2.44)
27
Note that constants can be excluded from the potential.Therefore,W
1
(s) can
be identified as
W
1
(s) = −i Ω

b
2
s

D+
D

+
b
4
s
2

= −iΩ

b
2
s

iR

e

−e
−iδ

+2x
0


+
b
4
s
2

.(2.45)
The Kutta condition [27] states that the velocity has to be finite at s =b where
ds/dz diverges,giving
dW
ds




s=b
= 0 ⇒



1
b
= −V
abs
e
i(−β+θ−δ)
+V
abs
e
−i(−β+θ−δ)

i Ω

iR

e

−e
−iδ

+2x
0
+2b


= −V
abs
2i sin(−β+θ−δ) −2i Ω(−Rsinδ+x
0
+b).(2.46)
Use the reference case of a static wing



1
b
= V
ref
2i sinα⇒
sinα =
−V
abs
sin(−β+θ−δ) −Ω(−Rsinδ+x
0
+b)
V
ref
,(2.47)
and that for small values of x
0
V
re f


(V
abs
cos(θ−β) +ΩR)
2
+V
2
abs
sin
2
(θ−β).(2.48)
Now,redefine x
0
in terms of the chord as x
0
=x
0r
c,which means that x
0r
=
0.25 is the quarter chord position and use that 4b =c
sinα =
−V
abs
sin(θ−β−δ) −Ω

−Rsinδ+x
0r
c+
c
4

V
ref
=
(V
abs
cos(θ−β) +ΩR)sinδ−V
abs
cosδsin(θ−β)
V
ref

Ω

x
0r
c+
c
4

V
ref
= sin

δ+arctan
−V
abs
sin(θ−β)
V
0
cos(θ−β) +ΩR


Ω

x
0r
c+
c
4

V
ref
.(2.49)
Assuming small angles of attack,one can approximate
arcsin(α+β) ≈ α+β (2.50)
28
which gives the simplifications
α = δ+arctan
−V
abs
sin(θ−β)
V
abs
cos(θ−β) +ΩR

Ωx
0r
c
V
re f

Ωc
4V
ref
.(2.51)
Note that at the position θ = 0,the blade is at the position θ
b
= π/2.The
substitution θ =θ
b
−π/2 gives
α=δ+arctan
V
abs
cos(θ
b
−β)
V
abs
sin(θ
b
−β) +ΩR

Ωx
0r
c
V
re f

Ωc
4V
ref
,(2.52)
which with β = 0 would correspond to equations (2.26) and (2.28),but in-
cludes mounting position x
0r
and flowcurvature.As an example,for a turbine
with chord 0.25 mand radius 3 m(i.e.the experimental turbine in Marsta [16]),
at very high rotational velocities

V
ref
≈ΩR

,the change in angle of attack
due to flow curvature is approximately -1.2 degrees.This gives higher angles
of attack upstreamand lower downstream.
29
3.Control strategy for vertical axis turbines
Optimizing the power from a wind/marine current turbine does not only re-
quire that the turbine is designed with the highest possible power coefficient.
Another important factor is to make sure that the turbine actually runs at the
tip speed ratio associated with the peak power coefficient.Therefore,a con-
trol strategy which keeps the turbine tip speed ratio near this optimal value is
preferable.
Wind turbines in general are either controlled by pitch or stall regulation,
where the most common design today is a horizontal axis turbine with pitch
regulation [28].The advantage with pitch control is that it introduces an ad-
ditional parameter that can be controlled,allowing for a more flexible control
system.Pitch control is mainly used in the region above rated wind speed (see
figure 3.1) to keep a smoother power and reduce mechanical loads,as the pitch
angle can be changed to reduce the blade forces [1,Chapter 8].Stall regula-
tion,instead,reduces the tip speed ratio,which increases the angle of attack.
This will eventually cause stall,which reduces the lift force and increases the
drag force.This will be most prominent for the tangential force,and thereby
the turbine torque,due to the significant increase in the drag force.Pitch con-
trol has been used for vertical axis turbines,mainly to improve performance at
lowtip speed ratios [29,30] where stall is avoided by actively altering the pitch
angle to reduce the angle of attack.The angle of attack oscillates between
positive and negative values as the blade moves between the upstream and
downstream section of the turbine.Hence,reducing the angle of attack with
active pitch requires a change in pitch angle during each revolution.An active
pitch mechanismwould complicate the turbine further.No pitch mechanisms
are included in the turbines studied here to reduce the sources of mechanical
failure.
Without a pitch mechanism,the remaining parameter to control is the ro-
tational velocity,where the turbine power is controlled by regulating the tip
speed ratio and thereby the power coefficient.Even with a pitch mechanism
installed,it is common for horizontal axis wind turbines to use a fixed pitch
angle in the variable rotational speed region illustrated in figure 3.1 [1,Chap-
ter 8].The following sections will focus on the variable rotational speed re-
gion,where the aimis to maximize the extracted power.
3.1 Control of a single turbine
One way to control a turbine is to perform real time flow velocity measure-
ments and adjust the tip speed ratio to optimal values.However,this would,
30
0
2
4
6
8
10
12
14
16
0
0.5
1
1.5
Fractionofratedpower
Wind speed (m/s)
Cut in
Variable rotational speed
region,constant C
P
Constant
power
Constant
rotational
speed
Rated wind
speed
Figure 3.1.Example of the different control strategies for different flow velocities.
The turbine starts operating at wind speed 3 m/s and operates at optimal power coeffi-
cient (approximately constant tip speed ratio) up to rated rotational velocity (wind
speed 10 m/s),where the rotational velocity is kept constant until rated power is
achieved (wind speed 12 m/s).At higher wind speeds,power is kept constant.There
is also a cut-out wind speed where the turbine is stopped (not shown in the figure).
Rotational velocity
Power
Speedincreases
Turbine stops Stable region
Equilibrium
Pe
P
P
e
Speeddecreases
Speeddecreases
Figure 3.2.Illustration of a control system,where the extracted power only depends
on the rotational velocity.
31
rely on the accuracy of the flow measurements.An alternative approach is
to let the extracted power P
e
be a function only of the rotational velocity of
the turbine.This type of control has been used for horizontal axis wind tur-
bines [31],but due to the low power coefficient at low tip speed ratios for
vertical axis turbines,some care has to be taken when transferring this control
system to a vertical axis turbine to avoid that the turbine ceases to rotate for
low tip speed ratios.
An example of a control system only using the rotational speed as input
parameter is illustrated in figure 3.2.The angular acceleration
˙
Ωof the turbine
is
˙
Ω=
P−P
e
J Ω
,(3.1)
where J is the moment of inertia of the system.Here,extracted power P
e
in-
cludes both power from the generator and electrical and mechanical losses in
the system.If the extracted power P
e
is less than the turbine power P,the tur-
bine accelerates,while if the extracted power is larger,the rotational velocity
decreases.In Figure 3.2,the turbine would stop if the turbine has too low
rotational velocity as the turbine power at this low rotational velocity is very
limited (P <P
e
).This region where the turbine will stop is characterized by
a low power coefficient at low tip speed ratios,which may apply to a vertical
axis turbine.The existence of such a region depends on the power absorption
characteristics of the turbine and for high enough wind speeds,this region
typically becomes smaller.
Figure 3.2 shows that there will be an equilibrium where extracted power
equals turbine power (P =P
e
).If this equilibrium occurs at the peak of the
power curve in figure 3.2,maximumenergy is extracted.If the extracted power
is normalized the same way as turbine power,the extracted power coefficient
is
C
P
e
=
P
e
1
2
ρAV
3

.(3.2)
Combining equations (3.1) and (3.2) gives for the equilibrium

where
˙
Ω=0

C
P
e
=C
P
.(3.3)
Equation (3.2) can be written in terms of tip speed ratio and rotational velocity
as
C
P
e
=
P
e
λ
3
1
2
ρA(RΩ)
3
.(3.4)
Denote λ
max
the tip speed ratio with the peak power coefficient C
Pmax
.By
choosing λ
max
as the desired equilibrium,equations (3.3) and (3.4) give
C
Pmax
=
P
e
λ
3
max
1
2
ρA(RΩ)
3
⇒P
e
=
1
2
ρAC
Pmax


λ
max

3
=k
2
Ω
3
,(3.5)
32
0
1
2
3
4
5
6
7
8
-0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Tip speed ratio
Powercoefficient
C at 3 m/s
P
C at 6 m/s
P
C at 12 m/s
P
C
P
e
Figure 3.3.The control strategy A from equation ( 3.5) illustrated for three different
flow velocities.
where the constant k
2
is related to λ
max
and C
Pmax
according to equation (3.5).
Under the conditions that C
Pmax
and λ
max
are constants,maximumpower will
be extracted if the extracted power is chosen to vary with the cube of the rota-
tional velocity and the constant k
2
is chosen according to equation (3.5).This
control strategy will be denoted “strategy A”.With C
Pmax
and λ
max
constant,
equation (3.5) is independent of flow velocity.However,due to the increased
Reynolds number at higher flow velocities,the airfoil performance will in-
crease,as drag is reduced and stall angle is increased.Therefore,the values
for λ
max
and C
Pmax
have a small dependence on the Reynolds number,see
figure 3.3.The performance of the strategy in equation (3.5) is plotted in fig-
ure 3.3 for three flow velocities (for specifications on the simulated turbine,
see section 5.2.6).The change in power coefficient with respect to the flow
velocity causes the obtained equilibrium tip speed ratio λ
e
to be slightly dis-
tanced to λ
max
,although the obtained power coefficient is very similar toC
Pmax
(see figure 3.3).For a flowvelocity of 3 m/s,there is an unstable region for tip
speed ratios below2.2,while for 6 m/s and 12 m/s,the strategy is stable in the
entire interval as long as the turbine power is larger than the mechanical and
electrical losses.Even though the strategy is stable for 6 m/s and 12 m/s,the
difference between turbine power and extracted power is low at low tip speed
ratios,causing a slow acceleration,cf.equation (3.1).
Two modifications to the devised strategy in equation (3.5) have been made,
with the intention to increase the rotational velocity at low tip speed ratios in
33
0
1
2
3
4
5
6
7
8
-0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Tip speed ratio
Powercoefficient
C at 3 m/s
P
C at 6 m/s
P
C at 12 m/s
P
C at 3 m/s
P
e
C at 6 m/s
P
e
C at 12 m/s
P
e
Figure 3.4.Control strategy B from equation (3.6) illustrated for three different flow
velocities.
0
1
2
3
4
5
6
7
8
-0.05
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Tip speed ratio
Powercoefficient
C at 3 m/s
P
C at 6 m/s
P
C at 12 m/s
P
C at 3 m/s
P
e
C at 6 m/s
P
e
C at 12 m/s
P
e
Figure 3.5.Control strategy C from equation (3.7) illustrated for three different flow
velocities.
34
order to create a more stable and faster control strategy.Strategy B is given as

P
e
=0 Ω≤Ω
0
,
P
e
=k
1
Ω
2
(Ω−Ω
0
) Ω>Ω
0
,
(3.6)
where the rotational velocity Ω
0
is required before the strategy starts to extract
energy.This strategy obtains higher differences between extracted and turbine
power for low rotational velocities at the cost of having λ
e
further distanced
from λ
max
.
Strategy C combines a higher torque at low rotational velocities with the
good performance of strategy A by dividing the rotational velocities in three
regions according to









P
e
=0 Ω≤Ω
0
,
P
e
=k
1
Ω
2
(Ω−Ω
0
) Ω
0
<Ω≤Ω
1
,
P
e
=k
2
Ω
3
Ω
1
<Ω≤Ω
2
,
P
e
=k
3
Ω
2
(Ω−Ω
2
) +k
2
Ω
2
Ω
2
Ω
2
<Ω,
(3.7)
where
k
1
=
Ω
1
Ω
1
−Ω
0
k
2
due to continuity at Ω=Ω
1
.This control strategy extracts minimum energy
up to Ω
0
,where Ω
0
for both strategy B and C is chosen to correspond to the
rotational velocity at tip speed ratio 4 and flowvelocity 3 m/s.These particular
values are chosen because 3 m/s is considered a reasonable cut-in velocity and
tip speed ratio 4 is the optimal value for the chosen turbine.
The control strategy C is designed to operate according to equation (3.5)
between Ω
1
and Ω
2
.The region between Ω
0
and Ω
1
is the increase from 0 to
the curve in equation (3.5).Hence,with increasing rotational flowvelocity,the
tip speed ratio decreases in the region between Ω
0
and Ω
1
until it has reached
its desired value for high power capture.At rotational velocities above Ω
2
,
the tip speed ratio is intentionally decreased to reduce mechanical loads on
the turbine,which creates a softer transition to the constant speed region in
figure 3.1.This region is reasonable to implement on a real turbine,but it has
not been studied in detail in this work and for simplicity,k
1
=k
3
was chosen.
The control systemis designed to reach Ω
2
at 9 m/s.This strategy is illustrated
in figure 3.5,where it is seen that at 3 m/s,λ
e
is higher that the peak value,
causing a small drop in power coefficient,although the difference is smaller
than for strategy B.For 12 m/s,λ
e
is slightly lower than the peak value due
to the intentional decrease in rotational velocity above Ω
2
.Note that the value
for the constant k
1
differs between strategy B and C,while the value for the
constant k
2
is equal for strategy A and C.Simulations of the three control
strategies are found in section 5.2.6.
35
3.2 Extension to multiple turbines
It is common to locate several turbines in close proximity to each other in a
farmconfiguration,as this can give economical benefits due to synergy effects,
since some parts of the systemcan be commonly used.Therefore,two differ-
ent electrical topologies are suggested,where the second “mutual topology”
is designed to reduce the number of electric components [32].For additional
possible topologies,see e.g.[33–35].
One way to control multiple turbines is to apply the model described in sec-
tion 3.1 for each turbine in the farm.This can be accomplished by using an
individual electric system for each turbine,where each turbine has a passive
diode rectifier and an inverter (“separate topology”).An alternative approach
is to connect all turbines to the same inverter,but to obtain individual control.
This would for a permanent magnet synchronous generator require some addi-
tional electrical components such as an active rectifier or a DC-DC converter.
One simplification is to connect all turbines to the same DC-bus with pas-
sive rectifiers,which makes the inverter the only parameter to control (“mutual
topology”).This reduces the number of electric components and it is therefore
of interest to study how performance is affected by this simplification.
For the mutual topology,the total power extracted from the entire turbine
farm is chosen as
P
e,tot
=
N
t

i=1
P
e

i
),(3.8)
where P
e

i
) is calculated according to one of the strategies A – C in sec-
tion 3.1.With the total extracted power chosen according to equation (2.1),
the separate and mutual topologies extract the same power if all turbines ex-
perience identical flow velocities.
Other studies have previously been performed for horizontal axis turbines,
with a global control strategy for the entire farm [36,37].Farm effects have
been included in the control system in [38,39] and a method for estimating
the flow velocity within the farm with only rotational velocities and extracted
power as input parameters is presented in [40].
36
4.Simulation models
The flow through a turbine can be simulated by solving Navier-Stokes equa-
tions coupled with the continuity equation.Navier-Stokes equations are valid
for all Newtonian fluids,i.e.both air and water,although in water,cavitation
would require additional modeling.The validity of Navier-Stokes equations
is well established,and the most accurate simulation model would be a direct
solution of these equations.However,Navier-Stokes equations are non-linear
partial differential equations and due to the computational complexity,accu-
rate direct solutions of these equations are numerically very hard to obtain,
even for the two-dimensional case.Considering that the fluid-flow in a large
turbine generally is turbulent,one common simplification is to model the tur-
bulence by e.g.Reynolds-averaged Navier-Stokes equations (RANS),which
include additional turbulence models.This has been done by e.g.[41,42].
The most common methods to solve Navier-Stokes equations are the Finite
Element Method (FEM) and the Finite Volume Method (FVM),which are
both based on a mesh over the entire simulated volume.For this reason,
these methods are suited for confined regions.There are commercial FEM
and FVMsimulation software available,but due to the long simulation times,
other methods have been chosen in this work.
Another method for solving Navier-Stokes equations is to use a vortex
method (see section 4.2).This method can also directly solve Navier-Stokes
equations,and is designed for open flows with small boundaries,which is the
case for a vertical axis turbine.Directly solving Navier-Stokes equations with
this method requires much computational time,and has not been done in this
work.An advantage of the vortex method is that there are many simplifica-
tions available,and by including models of blade forces,the simulation times
can be reduced substantially.This simplification is used within this work.
One even more simplified and very fast method is the streamtube model,
which assumes static flow,and does not give a full description of the flow
through the turbine.Due to the high computational efficiency,this model can
be used to quickly performmany simulations,and it is one of the models used
in this work,even though the simplifications of the flowlimit the applicability
of the model.
4.1 Streamtube models
Streamtube models are among the fastest models used for simulating vertical
axis turbines.The first application of streamtube models,usually credited to
37
V

V
d
V
e
A
d
R
θ
b
Δθ
Ω
c
Figure 4.1.Illustration of the double multiple streamtube model of Paraschivoiu with-
out flow expansion,where the streamtubes are illustrated by the horizontal lines and
the vertical line in the center illustrates the transition fromthe upstreamto the down-
stream part.Flow velocities are given for the upstream part.It is assumed that the
height of the streamtube Δz is included in the streamtube area A
d
.
Templin [43],used one streamtube for the entire turbine.The disadvantage
of this method is that the flow velocities at the blade positions are constant
throughout the entire turbine.A vertical axis turbine will interact with the
flow two times,one upstreamand one downstream.To overcome the problem
with equal flowvelocity when the blade is up- and downstream,Lapin [44] in-
troduced the double streamtube model,consisting of one streamtube upstream
and another downstream,allowing for lower flow velocities downstream.For
a real turbine,the velocity does not only vary between the upstreamand down-
streampart of the turbine,but could also vary over the entire swept area of the
turbine.To account for this,Strickland introduced the idea of using multiple
streamtubes [45].In this case,the systemhas to be solved for each individual
streamtube.
In Strickland’s implementation,the double streamtube model of Lapin was
not included.Later,both these two ideas were combined into the double
multiple streamtube model.Two different versions of this model were de-
veloped around the same time,one by Paraschivoiu [46] and one by Read and
Sharpe [47].The difference between the models is that the implementation
by Paraschivoiu assumes straight streamlines including only a velocity in the
x-direction (i.e.the direction of the asymptotic flow),while the model by Read
and Sharpe assumes that the streamlines expand linearly through the turbine.
This gives the difference that in the model by Paraschivoiu,all streamtubes are
independent of each other,while in the model of Read and Sharpe,they are
coupled.Both models have shown quite good agreement with experiments for
the curved blade Darrieus turbine [3,48].
38
4.1.1 Description of model
The double multiple streamtube model implementation used here is based on
the implementation of Paraschivoiu [3] and an illustration of this model is
given in figure 4.1.In the double multiple streamtube model,the turbine is
separated into two discs,one upstreamand one downstream.In the middle of
the turbine,the pressure is assumed to be the same as the asymptotic pressure.
With this assumption,the velocity at the upstream disc will be the average of
the velocity at the middle and the asymptotic velocity,similar to equation (2.7)
in the Betz theory.In the same way,the velocity at the downstream disc be-
comes equal to the average of the velocity at the middle,and the velocity far
behind the turbine.In the implementation by Paraschivoiu,the upstream disc
is solved first,independently of the downstream disc.The downstream disc
uses the velocity in the middle of the turbine,calculated when solving the
upstream disc,as input.Except from this calculation of velocity,both discs
are solved in the same way,and most equations will only be presented for the
upstreampart of the system.
This description will followthe same notation as Paraschivoiu for the direc-
tion of the flow velocity.Hence,the flow direction illustrated in figure 4.1 is
considered positive direction.The definitions of the blade angle θ
b
and the an-
gle of attack α will still remain the same as in section 2.3.Therefore,negative
angles of attack will be obtained when the blades are upstream,while the de-
scription by Paraschivoiu [3,Chapter 6] has positive angle of attack upstream.
The difference in the definition of the flow direction is the reason why the
expressions for flow velocity and angle of attack are different in this section,
compared to section 2.3.
The main principle behind streamtube models is the use of momentumcon-
servation in each streamtube.Similar to the Betz derivation,momentum con-
servation gives
ρA
ei
V
2
ei
−ρA
∞i
V
2
∞i
=ρA
di
V
di
(V
ei
−V
∞i
) =F
xi
,(4.1)
where A
∞i
,A
di
and A
ei
are the cross-sectional areas of a streamtube and i
denotes the streamtube index (cf.figure 2.1).The difference,compared to
the Betz theory,is that in this model,the force F
xi
is calculated from blade
section data.If the streamtubes are discretized to give each streamtube the
same angular distance Δθ and height Δz,the size of each streamtube becomes
A
di
=R
i
ΔzΔθ|cosθ
bi
|.(4.2)
To calculate the force,the first step is to obtain the lift and drag coefficients,
which depend on the angle of attack and the Reynolds number.In the model
of Paraschivoiu,it is assumed that the direction of the flow does not change.
This gives the relative velocity as
V
ri
=V
di


R
i
Ω
V
di
−sinθ
bi

2
+cos
2
θ
bi
cos
2
η
i
(4.3)
39
and the angle of attack as
α
i

i

i
=arctan

cosθ
bi
cosη
i
sinθ
bi

R
i
Ω
V
di


i
,(4.4)
where η is the angle of the blade relative to the vertical axis.
When the coefficients have been determined from the empirical data,the
lift and drag forces are obtained from
F
Li
=
1
2
C
Li
ρc
i
l
i
V
2
ri
,(4.5)
F
Di
=
1
2
C
Di
ρc
i
l
i
V
2
ri
,(4.6)
where c
i
is the chord and l
i
is the blade length in the streamtube.Generally,
it is the normal and tangential forces that are of interest.These forces can be
determined from the corresponding normal and tangential force coefficients
C
Ni
= C
Li
cosϕ
i
+C
Di
sinϕ
i
,(4.7)
C
Ti
= C
Li
sinϕ
i
−C
Di
cosϕ
i
(4.8)
and therefore,the forces are
F
Ni
=
1
2
C
Ni
ρc
i
l
i
V
2
ri
,(4.9)
F
Ti
=
1
2
C
Ti
ρc
i
l
i
V
2
ri
.(4.10)
Note that by this definition,the normal force is perpendicular to the blade,and
the force in the radial direction is
F
Ri
=F
Ni
cosη
i
.(4.11)
To calculate how the velocity decreases in the streamtube,the forces in the
x-direction (parallel to the flow) are of interest.These are obtained as
F
xi
=F
Ni
cosθ
bi
cosη
i
−F
Ti
sinθ
bi
.(4.12)
Considering that there are N
b
blades on the turbine,and that a streamtube only
have a blade inside it a fraction of the time,the mean force is given by
F
xi
=
N
b
Δθ

(F
Ni
cosθ
bi
cosη
i
−F
Ti
sinθ
bi
).(4.13)
If the length of the blade in a streamtube l
i
=Δz/cosη
i
is inserted into equa-
tions (4.9) and (4.10),equation (4.13) becomes
F
xi
=
N
b
ρc
i
V
2
ri
ΔθΔz


C
Ni
cosθ
bi
−C
Ti
sinθ
bi
cosη
i

.(4.14)
40
V
s
V
di
V
ei
c
s
Δr
s
R
outer
R
inner
θ
b
θ
Streamtube i
Ω
η
s
r
s
Δz
Figure 4.2.Illustration of the parameters used for the strut model.
Combining this with equation (4.1) and (4.2) gives
ρR
i
ΔzΔθ|cosθ
bi
|V
di
(V
ei
−V
∞i
) =
N
b
ρc
i
V
2
ri
ΔθΔz


C
Ni
cosθ
bi
−C
Ti
sinθ
bi
cosη
i

(4.15)
which with the substitutions
V
di
=u
i
V
∞i
(4.16)
and
V
di
=
V
∞i
+V
ei
2
(4.17)
gives
1−
1
u
i
=
N
b
c
i
V
2
ri
8πR
i
|cosθ
bi
|V
2
di

C
Ni
cosθ
bi
−C
Ti
sinθ
bi
cosη
i

.(4.18)
This is the equation to be solved to obtain the interference factors u
i
(defined
in equation (4.16)).Since V
r
,C
N
and C
T
depend on u,this equation has to be
iterated.It should be noted that from equations (4.16) and (4.17),the velocity
in the center is
V
ei
=(2u
i
−1)V
∞i
.(4.19)
4.1.2 Including struts
The above theory only includes the blades,but for an H-rotor,it is necessary
to attach the blades to the main shaft with struts.These struts interact with
the flow,affecting turbine performance.One strut model was implemented by
Moran [49],who derived a correction coefficient for the drag.This is a fast
solution suitable for horizontal struts,but for struts with an angle,a lift force
41
will also be generated.To increase the flexibility of the code,it was chosen to
calculate the forces directly fromthe lift and drag coefficients instead.
To calculate the forces,the relative velocity has to be known.Let r
s
be a
position along the strut
R
inner
<r
s
<R
outer
,(4.20)
where the strut starts at R
inner
and ends at R
outer
,see figure 4.2.The struts
are discretized into Δr long segments in the radial direction.In the streamtube
model,the velocity is only known at the turbine disc and at the center of the
turbine.The simple approximation that the velocity varies linearly between
these positions gives the velocity V
s j
of strut segment j as
V
s j
=V
ei
+
V
di
−V
ei

R
2
i
−r
2
s j
sin
2
θ
b
r
s j
cosθ
b
.(4.21)
The position of the blade,to which the strut is attached,is given as θ
b
.This
expression involves the velocity in streamtube i,in which the segment is lo-
cated.To determine which streamtube this is,the corresponding angle at the
circumference of the turbine is given by
θ
i
=arcsin

r
s j
R
j
sinθ
b

.(4.22)
In the vertical direction,the position of segment j will be equal to the vertical
position of the corresponding streamtube i.The relative velocity can be calcu-
lated basically in the same way as in equation (4.3),but since the discretization
is performed in the radial direction for the struts,the angle is to be given rel-
ative this direction as well (the reason for not using vertical discretization for
struts is that
F
xi
in equation (4.14) would diverge for horizontal struts).With
the strut angle given as η
s
,the velocity is given by
V
rs j
=V
s j


r
s j
Ω
V
s j
−sinθ
b

2
+cos
2
θ
b
sin
2
η
s j
(4.23)
and the angle attack α
s
for the strut is
α
s j
=arctan
cosθ
b
sinη
s j
sinθ
b

r
s j
Ω
V
s j
,(4.24)
where zero pitch angle is assumed for the strut.By defining the chord of the
strut as c
s
and following the same derivation as with equations (4.7) – (4.14),
the expression for the mean force in the flow direction for a section (Δr
s
,Δθ)
of the strut is obtained as
F
xs j
=
N
b
c
s j
ρV
2
rs j

ΔθΔr
s

C
Ns j
cosθ
b
tanη
s j
−C
Ts j
sinθ
b
cosη
s j

,(4.25)
42
where C
Ns j
and C
Ts j
are the normal and tangential force coefficients of the
strut,which are defined in a similar way as for the blades
C
Ns j
= C
Ls j
cosα
s j
+C
Ds j
sinα
s j
,(4.26)
C
Ts j
= C
Ls j
sinα
s j
−C
Ds j
cosα
s j
.(4.27)
Here,C
Ls j
and C
Ds j
are the lift and drag coefficients for the struts.The differ-
ences in equations (4.24) and (4.25),compared to equations (4.4) and (4.14),
originate fromthe different definition of η
s j
,and that the discretization is done
in the radial direction.
The interference factor u
i
can be calculated according to
1−
1
u
i
=
1
2R
i
|cosθ
b
|V
2
di
ρΔθΔz

F
xi
+

j∈i
F
xs j

(4.28)
where ∑
j∈i
F
xs j
is the sum of the forces from all points (r
si

b
) that corre-
sponds to streamtube i.Similar expressions can be derived for the downwind
part of the turbine.The torque from the struts can be calculated in a similar
way as for the blades in a curved blade Darrieus turbine
T
s

b
) =

j
r
s j
F
Ts j

b
)
cosη
j
,(4.29)
where
F
Ts j

b
) =
1
2
C
Ts j

b
)ρc
s j
Δr
s
V
rs j

b
)
2
.
4.1.3 Obtaining lift and drag coefficients
In the streamtube model,it is necessary to have the lift and drag coefficients
for a given angle of attack.Experimental data can be used to obtain these if
the Reynolds number and angle of attack are known.For several symmetrical
NACA profiles,such data can be obtained from [50],which has been used in
this work.The problemwith this kind of data is that it usually is valid for very
long blades and a static angle of attack.For a vertical axis turbine,the blades
will change their angles of attack during the revolution.
When the angle of attack for a blade increases above the stall angle,the flow
will separate fromthe blade surface,which reduces the lift force and increases
the drag.The flow separation takes time to develop and if the blade rapidly
increases its angle of attack,the lift force can be higher than in the static case
for a short period of time.During this time,a vortex starts to form at the
surface of the blade and when this vortex detaches,the lift force is greatly
reduced.This phenomenon is called dynamic stall.When the turbine has a
low tip speed ratio,the angles of attack become high enough for the blades to
experience this in each revolution.This occurs for tip speed ratios around 4.
43
To model dynamic stall,the Gormont model has been used.This model was
originally developed by Gormont [51],later modified by Massé [52] and ad-
justed by Berg [53].The advantages of this model is that it only requires data
for lift and drag coefficients,flow velocity,blade thickness to chord ratio,an-
gle of attack and rate of change of the angle of attack,making it easy to apply
to any blade,which is the reason why this model was chosen.Other mod-
els have been shown to give better results [3].One example is the Beddoes-
Leishman model [54,55],but it has more parameters that have to be calibrated,
which makes it harder to apply the model to an arbitrary airfoil.
Tip effects
For a blade to obtain a lift force,there has to be circulation around the blade
(cf.Kutta Joukowski lift formula,equation (4.57)).Since the vorticity is di-
vergence free

ω=∇×

V ⇒∇·

∇×

V

=0

,there are vortices generated
from the blade tips with the same circulation as around the blade.These vor-
tices are the source of induced drag,and the corresponding losses.It should be
noted that the curved blade Darrieus turbine does not have any distinct blade
tips (the blades are attached to the main axis).Here,the tip vortices have
to leave the blades more evenly distributed over the length of the blade,and
in the model of Paraschivoiu,no tip losses are applied to curved blade Dar-
rieus turbines.For straight blade turbines,one correction model,derived from
Prandtl’s theory for screwpropellers,used by e.g.Sharpe [48],uses a velocity
correction factor F,which for the upstreamdisc is
F
i
=
arccos

e

πa
i
s
i

arccos

e

πh
2s
i
,(4.30)
where
s
i
=
πV
ei
N
b
Ω
,(4.31)
a
i
=
h
2
−|z
ai
|,(4.32)
z
ai
is the altitude (with zero defined in the center of the turbine) and h is the
turbine height.A corresponding expression exists for the downstream disc.
The new expressions for velocity and angle of attack are
V
ri
= V
di


R
i
Ω
V
di
−sinθ
bi

2
+F
2
i
cos
2
θ
bi
cos
2
η
i
,(4.33)
α
i
= arctan
F
i
cosθ
bi
cosη
i
sinθ
bi

R
i
Ω
V
di

i
.(4.34)
44
In addition,Paraschivoiu [3] suggests using finite wing theory as well to cal-
culate induced drag and reduced angle of attack.This is given by
C
Li
=
C
L∞i
1−
a
0i
πAR
i
,(4.35)
a
0i
= 1.8π

1+
0.8t
bi
c
i

,(4.36)
C
D
i
= C
D∞i
+
C
2
Li
πAR
i
,(4.37)
α
bi
= α
i

C
Li
πAR
i
,(4.38)
where AR is the aspect ratio,t
b
the blade thickness,C
L∞
and C
D∞
lift and
drag coefficients for an infinitely long blade and α
b
the corrected angle of
attack.Paraschivoiu suggests that equations (4.35) – (4.38) should not be ap-
plied to dynamic stall calculations,but since this would lead to a discontinuity
in the limit when dynamic stall starts to have an effect,the current implemen-
tation uses equations (4.35),(4.36) and (4.38),then calculates the dynamic
stall corrections to the lift coefficient,and then calculates the induced drag
using equation (4.37) with the new lift coefficient.It should be noted that
equations (4.35) – (4.38) are derived using the assumption that the angle of at-
tack is constant,the tip vortices extend along a straight line behind the blade,
and the wing is elliptical,which is not the case for a vertical axis turbine.This
has to be taken into account when evaluating the accuracy of the model for
straight blade turbines.
4.1.4 Corrections due to flow curvature
As was illustrated in section 2.4,there are effects due to flow curvature that
should be taken into account.With the definition of the velocity as in the
streamtube model (V →−V),the assumption of straight streamlines (β =0),
using that the reference velocity V
ref
is the relative velocity V
r
and using the
velocity at the turbine disc V
d
instead of V
abs
,equation (2.52) gives
α=δ+arctan
V
d
cosθ
b
ΩR−V
d
sinθ
b

Ωx
0r
c
V
r

Ωc
4V
r
.(4.39)
This is included in the adaptation of Sharp [48],but with the assumptions that
the curvature force only acts in the normal direction and that the mounting
point is in the center of the blade (x
0r
=0).With the assumption
C
N
=a
s
α,(4.40)
the correction to the normal coefficient becomes
C
N
=C
N0

a
s
Ωc
4V
r
,(4.41)
45
z
y
a
y∞
a
ye
aze
az∞
bze
1
bz∞
b
y∞
1
b
ye
Figure 4.3.Illustration of the cross-section for a streamtube at the three positions:Far
ahead of turbine (black),at the upstreamdisc (white) and in the center (hatched).
where C
N0
is the original normal force without corrections.
4.1.5 Including flow expansion
In contrast with the model of Read and Sharpe [47],the model of Paraschi-
voiu does not include streamtube expansion in the basic implementation.One
method for implementing the expansion was used by Paraschivoiu [56] where
the streamtubes were allowed to expand,but the direction of the velocity was
not changed.In this implementation,similar to the one by Read and Sharpe,
it was assumed that streamtube expansion only occurs in the horizontal plane,
and that the streamline that intersects the turbine center is unchanged.In this