# Gundtoft2009-WindTurbines.pdf

by user

on
Category: Documents
5

views

Report

#### Transcript

Gundtoft2009-WindTurbines.pdf
```Wind turbines
1
Foreword
Wind Turbines
- Design of an optimal rotor i.e. pitch angle and chord length of the
blades and how to calculate the power production
- Including a spread sheet model in Excel and an EES-model for the
calculations
This paper (now 2. edition) deals with the general procedures for the design of the rotor blades –
pitch angle and chord length – for horizontal axis wind turbines (HAWT) and for the calculation of
power production in the rotor.
The paper is used in the course Fluid Dynamics (code: MIFLD1) at the Department of Mechanical
Engineering at the University College of Århus.
If you find errors in this report, you should send them to me, mail to: [email protected]
In this paper a comma is used as decimal point, i.e. π = 3,1416.
Compared to the 1. Edition (2004) Chapter 3 and Chapter 6 have been totally rewritten.
Also an EES-model of the BEM method is included – see appendix C
June 2009 / Søren Gundtoft
w
r
u= r
v1
R
v = (1-a)v1
FL
c
Rotor plane
Rotor axis
FD
2. Edition
Søren Gundtoft
University College of Aarhus
June 2009
Content
1. Introduction................................................................................................................................... 2
2. Power in the wind ......................................................................................................................... 2
3. Rotor design .................................................................................................................................. 5
3.1. Air foil theory – an introduction............................................................................................. 5
3.2. Pitch angle, β, and chord length, c, after Betz ........................................................................ 7
3.3. Pitch angle, β, and chord length, c, after Schmitz................................................................. 11
4. Characteristics of rotor blades ..................................................................................................... 16
5. The blade element momentum (BEM) theory ............................................................................. 20
6. Efficiency of the wind turbine..................................................................................................... 23
6.1. Rotor .................................................................................................................................... 23
6.2. Gear box, generator and converter........................................................................................ 26
7. Example, BEM............................................................................................................................ 28
8. Distribution of wind and annual energy production .................................................................... 30
9. Symbols ...................................................................................................................................... 33
10. Literature................................................................................................................................... 34
App. A: Conservation of momentum and angular momentum ........................................................ 35
App. B: Formulas, spread sheet calculations................................................................................... 39
App. C: Formulas, EES-programme ............................................................................................... 41
App. D: Formulas, EES-programme, to solve the integral in (6.4).................................................. 43
Wind turbines
2
1. Introduction
It is assumed that the reader knows some basic fluid mechanics, for example the basic theories of
fluid properties, the ideal gas law, Bernoulli’s equation, turbulent and laminar flow, Reynolds’s
number, etc. Important for the understanding of the theory is mastering the conservation laws for
momentum and angular momentum, for which reason these theories are presented in Appendix A.
Wind turbines
1
ρ v 2 + p = p tot
2
[Pa]
A3
A
v1
v3
Speed
v1
Pressure
All important calculations are demonstrated by examples. The calculations after the BEM method
can be done by simple spread sheet calculations and in Appendix B the formulas are printed. Also
the code for an EES-model of the BEM method is presented, see Appendix C.
Notation: In the following we will use index 1 for states “far up stream” the rotor plane, index for
the states in the rotor plane and index 3 for states far downstream. For simplicity index 2 – states in
the rotor plane - will be omitted in most cases.
A1
v1
Most of the theory is taken from ref./1/ and /4/. Data for rotor blade sections are taken from ref./2/
and /3/.
The question is: How much energy can be taken from the wind? The wind turbine decelerates the
wind, thereby reducing the kinetic energy in the wind. But the wind speed cannot be reduced to zero
– as a consequence, where should the air be stored? As first time shown by Betz, there is an
optimum for the reduction of the wind speed, and this is what is to be outlined in this chapter.
Figure 2.1 shows the streamlines of air through a wind turbine
(2.1)
where ptot is the total pressure, which is constant. That means, if the speed of flow goes up, the
pressure goes down and vice versa.
In this paper you will find:
• Chapter 2: How much energy that can be taken from the wind in an idealized wind turbine
(proof of Betz’ law) and how to design the rotor
• Chapter 3: How to design a an optimal rotor – pitch angle and chord length
• Chapter 4: Characteristic of rotor blades (coefficients of lift and drag)
• Chapter 5: How to calculate the power of a given rotor (the BEM theory)
• Chapter 6: Efficiency of a wind turbine
• Chapter 7: An example – BEM method in Excel
• Chapter 8: Distribution of the natural wind and calculation of annual energy production
2. Power in the wind
3
v
v3
p
p1
p3=p1
p - p
Figure 2.1: Interaction between wind and wind turbine
Assumption: The pressure changes are relatively small compared to the pressure in the ambient
(about 1 atm = 101325 Pa) therefore we assume the density to be constant.
If we use (2.1) for the flow up-stream of the rotor, we get
p1 +
1
1
ρ v1 2 = p + + ρ v 2
2
2
[Pa]
(2.2)
If we use (2.1) down stream of the rotor plane, we get
p + − ∆p +
1
1
ρ v 2 = p1 + ρ v3 2
2
2
[Pa]
(2.3)
[Pa]
(2.4)
Long in front of the rotor, the wind speed is v1. After passing the wind turbine rotor (called the rotor
in the following), the wind speed would be reduced to v3. The pressure distribution is as follows.
The initial pressure is p1. As the air moves towards the rotor, the pressure rises to a pressure p+ and
by passing the rotor, the pressure suddenly falls by an amount of ∆p i.e. the pressure is here p- = p+ ∆p. After passing the rotor, and far down stream the pressure again rises to p3 = p1. Curves for wind
speed and pressure are shown in figure 2.1.
Subtracting (2.3) from (2.2) we get
Bernoulli’s equation: If we look at the air moving towards the rotor plane, we can use the
Bernoulli’s equation to find the relation between the pressure p and the speed v, while we can make
the assumption that the flow is frictionless:
Change of momentum: This differential pressure can also be calculated on the basis of”change in
momentum”. (For information see Appendix A). If we look at one square meter of the rotor plane,
∆p =
(
1
ρ v1 2 − v3 2
2
)
Wind turbines
the mass flow equals ρ v. Momentum equals mass times velocity, with the unit N. Pressure equals
force per surface, then the differential pressure can be calculated as
T = ½ ρ v1 A C T
[N]
∆p = ρ v (v1 − v3 )
In figure 2.2, curves for CP and CT are shown.
Wind turbines
[Pa]
4
(2.5)
2
Now (2.4) and (2.5) give
1
(v1 + v3 )
2
[m/s]
(2.6)
This indicates that the speed of air in the rotor plane equals the mean value of the speed upstream
and down stream of the rotor.
Power production: The power of the turbine equals the change in kinetic energy in the air
P=
(
)
1
ρ v v1 2 − v3 2 A
2
[W]
(2.7)
Here A is the surface area swept by the rotor.
[N]
(2.8)
We now define ”the axial interference factor” a such that
v = (1 − a ) v1
[m/s]
(2.9)
Using (2.6) and (2.9) we get v3 = (1 – 2a) v1 and (2.7) and (2.8) can be written as
P = 2 ρ a (1 − a ) v1 A
[W]
(2.10)
T = 2 ρ a (1 − a )v1 A
[N]
(2.11)
2
3
2
We now define two coefficients, one of the power production and one of the axial forces as
C P = 4a(1 − a )
[-]
(2.12)
C T = 4a (1 − a )
[-]
(2.13)
[W]
(2.14)
2
Then (2.10) and (2.11) can be written as
P=
1
ρ v13 A C P
2
0,8
0,6
0,4
C_T
C_P
0,2
=16/27
0,0
0
0,1
0,2
0,3
0,4
0,5
a [-]
Figure 2.2: Coefficient of power CP and coefficient of axial force CT for an idealized wind turbine.
As shown, CP has an optimum at about 0,593 (exactly 16/27) at an axial interference factor of 0,333
(exactly 1/3). According to Betz we have
The axial force (thrust) on the rotor can be calculated as
T = ∆p A
(2.15)
1,0
C_P & C_T [-]
v=
5
PBetz = C p,Betz ½ ρ v1 A
3
with C P,Betz =
16
27
[W]
(2.16)
Example 2.1
Let us compare the axial force on rotor to the drag force on a flat plate? If a = 1/3 the CT = 8/9 ≈
0,89. Wind passing a flat plate with the area A would give a drag on the plate of
1
2
[N]
(2.17)
FD = C D ρ v1 A
2
where CD ≈ 1,1 i.e. the axial force on at rotor – at maximal power – is about 0,89/1,1 = 0,80 = 80%
of the force on a flat plate of the same area as the rotor!
3. Rotor design
3.1. Air foil theory – an introduction
Figure 3.1 shows a typical wing section of the blade.
The air hits the blade in an angle αA which is called the “angle of attack”. The reference line” for
the angle on the blade is most often “the chord line” – see more in Chap. 4 for blade data. The force
on the blade F can be divided into two components – the lift force FL and the drag force FD and the
lift force is – per definition – perpendicular to the wind direction.
Wind turbines
6
Wind turbines
7
3.2. Pitch angle, β, and chord length, c, after Betz
F
FL
Figure 3.3 shows the velocities and the angles in a given distance, r, from the rotor axis. The rotor
shown on the figure is with two blades, i.e. B = 2. To design the rotor we have to define the pitch
angle β and the chord length c. Both of them depend on the given radius, that we are looking at
therefore we sometimes write β(r) and c(r).
FD
Chord line
w
Figure 3.1: Definition of angle of attack
v1
(3.1)
Similar for the drag force
1
(3.2)
FD = C D ρ w 2 (bc )
2
The coefficient of lift and drag both depend of the angle of attack, see figure 3.2.
c
where CL is the “coefficient of lift”, ρ is the density of air, w the relative wind speed, b the width of
the blade section and c the length of the chord line.
FL
Rotor axis
Rotor plane
The lift force can be calculated as
1
FL = C L ρ w 2 (bc )
2
r
u= r
w
R
v = (1-a)v1
FD
Figure 3.3. Velocities and angles
Angles, that all depends on the given radius
• γ(r) = angle of relative wind to rotor axis
• φ(r) = angle of relative wind to rotor plane
• β(r) = pitch angle of the blade
For angles of attack higher than typically 15-20° the air is no longer attached to the blade, a
phenomenon called “stall”.
The ratio CL/CD is called the “glide ratio”, i.e. GR = CL/CD. Normally we are interested in at high
glide ratio for wind turbines as well as for air planes. Values up to 100 or higher is not uncommon
and the angles of attack giving maximum are typical in the range 5 – 10°.
The blade, as shown on the figure is moving up wards, thus the wind speed, seen from the blade, is
moving down wards with a speed of u. We have
w2 = v 2 + u 2
GR
0
15
30
45
alpha [°]
60
75
90
1,8
Betz does not include rotation of the wind, i.e. a’ = 0 (see definition of a’ later – formula (3.23)).
Therefore
180
C_L
1,6
160
C_D
1,4
140
GR
1,2
120
1,0
100
0,8
80
0,6
60
0,4
40
0,2
20
0,0
GR [-]
C_D
C_L & C_D [-]
180
160
140
120
100
80
60
40
20
0
C_L
GR [-]
C_L & C_D [-]
1,8
1,6
1,4
1,2
1,0
0,8
0,6
0,4
0,2
0,0
(3.3)
NACA 23012
NACA 23012
0
0
2
4
6
8
10
12
14
16
18
20
u =ωr
[m/s]
Here ω is the angular speed of the rotor given by
ω = 2π n
where n is the rotational speed of the rotor in round per second.
alpha [°]
Figure 3.2: Coefficient of lift and drag as a function of the angle of attack (left: 0<α<90°; right:
0<α<20°)
(3.4)
Now we define the “tip speed ration” i.e.
(3.5)
X =
v tip
v1
=
ωR
[-]
v1
8
Wind turbines
(3.6)
dF L
dUL
Wind turbines
dF L
9
dF L
Combining these equations we get
Torque
(3.9)
Chord length, c(r):
For the rotor plane (torque) we have
1
dU = ρ w 2 c dr C x
2
with
[N]
(3.12)
C x = C L sin (ϕ ) − C D cos(ϕ )
[-]
(3.13)
1
ρ w 2 c dr C y
2
[N]
(3.14)
C y = C L cos(ϕ ) + C D sin (ϕ )
[-]
(3.15)
For the rotor axis (thrust) we have
If we look at one blade element in the distance r from the rotor axis with the thickness dr the lift
force is, see formula (3.1) and (3.2)
1
ρ w 2 c dr C L
2
and the drag force
dFL =
[N]
1
ρ w 2 c dr C D
2
[N]
(3.10)
(3.11)
dT =
with
Now, in the design situation, we have CL >> CD, then (3.12) and (3.13) becomes
1
ρ w2 c d r CL cos(γ )
2
and then the power produced
R
dU =
[N]
(3.16)
dP = dU r ω
[W]
(3.17)
If we have B blades, (3.16) including (3.17) gives
1
dP = B ρ w 2 c d r C L cos(γ ) r ω
2
[W]
According to Betz, the blade element would also give
Thrust
y
Figure 3.5. Forces on the blade element decomposed on the rotor plane, dU (torque), and in the
rotor axis, dT (thrust)
2R
−αD
3r X
r dr
x
12.07.2008/SGt
where αD is the angle of attack, used for the design of the blade. Most often the angle is chosen to
be close to the angle, that gives maximum glide ration, see figure 3.2 that means in the range from 5
to 10°, but near the tip of the blade the angle is sometimes reduced.
dFD =
Rotor plane
(3.8)
dTD
dF D
dF D
dUD
(3.7)
dF D
and then the pitch angle
β (r )Betz = arctan
dTL
Rotor axis
3r X
γ (r ) = arctan
2R
or
2R
ϕ (r ) = arctan
3r X
(3.18)
Wind turbines
dP =
16 1 3
ρ v1 (2π r d r )
27 2
[W]
10
16π R
9 B C L,D
1
2
4
⎛r⎞
X X ⎜ ⎟ +
9
⎝R⎠
[m]
11
(3.19)
Using v1 = 3/2 w cos(γ) and u = w sin(γ), then (3.18) and (3.19) gives
c(r )Betz =
Wind turbines
(3.20)
3.3. Pitch angle, β, and chord length, c, after Schmitz
Schmitz has developed a little more detailed and sophisticated model of the flow in the rotor plane.
The torque M in the rotor shaft can only be established because of the rotation of the wake, cf.
Appendix A which is a result of the conservation law for angular momentum
v
2
u
where CL,D is the coefficient of lift at the chosen design angle of attack, αA,D.
Example 3.1
What will be shown later is that a tip speed ration of about X = 7 is optimal (see fig. 6.2). Further
more 3 blades seem to be state of the art. Figure 3.6 and 3.7 shows the results of formula (3.20)
concerning the chord length i.e. according to Betz.
[m/s]
v
Figure 3.6. Chord length as function of radius for X = 7 and for different numbers of blades
u
Figure 3.8. Down stream rotation of the wake – The wake rotates in the opposite direction to the
rotor
The power can be calculated as
P=Mω
[W]
(3.21)
where M is the torque in the rotor shaft and ω is the angular speed. According to the conservation
rule of angular momentum, the torque in the rotor shaft can only be established because of a swirl
induced in the slipstream in the flow down stream of the rotor. As for the axial speed v it can be
shown theoretically that the change in the tangential speed in the rotor plane is half of the total
change, i.e. we have in the rotor plane
u = r ω + ½ ∆u
Figure 3.7. Chord length as function of radius for three blades B = 3 and for different tip speed
ratios
or
[m/s]
(3.22)
As mentioned previously index 1 is used for the upstream situation, index 2 and 3 for rotor plane
and downstream respectively. In the following index 2 is some times omitted – for simplicity.
w
1
which defines the “tangential interference factor a’ “
v1
Rotor plane
r =u 1
(3.24)
(3.26)
[m/s]
(3.27)
b3)
v
b4)
1
Rotor plane
From figure 3.9 we further have
r
½ u
(3.28)
1
w
w
w3
w3
c1)
v3
[m/s]
far down stream
∆w = 2 w1 sin (ϕ1 − ϕ )
v
v1
w
1
Combining (3.25) and (3.26) we get
v = w1 cos(ϕ1 − ϕ )sin (ϕ )
½ u
v3
v = w sin (ϕ )
½ u
½ w
½ v
[m/s]
and from figure 3.9. – b2)
r
w
(3.25)
b2)
Rotor plane
r
½
[m/s]
b1)
Rotor plane
w = w1 cos(ϕ1 − ϕ )
w
The change in w1 is because of the air foil effect. If we assume that the drag is very low (compared
to lift, i.e. CD << CL => CD ≈ 0) then the ∆w vector is parallel to the lift force vector dFL (because
of the conservation law of momentum) and we – per definition of the direction of lift force – also
have that the ∆w vector is perpendicular to w – see figure 3.9-b4). Based on these considerations we
have the following geometrical relations
w
[m/s]
a)
1
Now look at the flow in the rotor plane, see figure 3.9. What is important here is the relation
r r
r
w = w1 + ½ ∆w
13
(3.23)
w
[m/s]
Wind turbines
up stream
u = r ω (1 + a')
12
w
½
Wind turbines
c2)
Rotor plane
r
r
u
u
12.07.2008/SGt
Figure 3.9. Speed in the rotor plane a) far upstream; b) in the rotor plane and; c) far down stream
Now, let us look at the power! From the conservation of momentum we have
dFL = ∆w dq
[N]
(3.29)
where dq is the mass flow through the ring element in the radius r with the width dr, i.e.
Wind turbines
dq = 2 ρ π r dr v
[kg/s]
14
Wind turbines
15
Optimal pitch angle
(3.30)
X = 5; B = 3; alfa_D = 7,0°; C_L = 0,88
Power equals “torque multiplied by angular velocity” and (neglecting drag) then
80
70
dP = dM ω
= ∆w dq sin (ϕ ) r ω
[ 2 ρ π r dr ) w1 cos(ϕ1 − ϕ )sin (ϕ )]sin (ϕ )rω
= {2 w1 sin (ϕ1 − ϕ )}(
[kg/s]
(3.31)
= r ω ρ 2π dr w sin [2 (ϕ1 − ϕ )]sin (ϕ1 )
2
2
1
beta [°]
= dFL sin (ϕ ) r ω
60
beta(Betz)
50
beta(Schmitz)
40
30
20
2
10
0
In the bottom transaction above we have used the relation sin(x) cos(x) = sin(2x).
-10
We have now a relation for the power of the ring element as a function of the angle φ but we do not
know this angle? The trick is now to solve the equation d(dP)/dφ = 0 to find the angle that gives
maximum power. Doing this for (3.31) we get
d (dP )
= r 2ω ρ 2π dr w12 − 2 cos[2(ϕ1 − ϕ )]sin 2 ϕ + 2 sin[2(ϕ1 − ϕ )]sin ϕ cos ϕ
dϕ
(
)(
(
(
)
)
= r 2ω ρ 2π dr w12 2 sin{sin[2(ϕ1 − ϕ )]cos ϕ − cos[2 (ϕ1 − ϕ )]sin ϕ }
)
[W/°] (3.32)
v1
ωr
=
2
R
arctan
3
Xr
(3.33)
(3.34)
2
3
R
−αD
rX
1,0
Note, that only for small r/R the two theories differ. And here the power produced is small because
of the relatively small swept area. At the tip (r/R = 1) the optimal angle is approx. 0,5° for both.
⎛ ϕ ⎞ ⎛ ϕ ⎞ ⎛ 2ϕ ⎞
2
= 2 w1 2 ρπ r dr sin ⎜ 1 ⎟ cos⎜ 1 ⎟ sin ⎜ 1 ⎟
⎝ 3⎠ ⎝ 3⎠ ⎝ 3 ⎠
[N]
(3.36)
⎛ϕ ⎞
⎛ϕ ⎞
2
= 2 w1 2 ρπ r dr sin 2 ⎜ 1 ⎟ cos 2 ⎜ 1 ⎟
3
⎝ ⎠
⎝ 3⎠
where we again use sin(2x) = 2 sin(x)cos(x).
and the for pitch angle
β (r )Schmitz = arctan
0,8
= 2 w1 sin (ϕ1 − ϕ ) 2 ρπ r dr (w1 cos(ϕ1 − ϕ )sin (ϕ ))
or
2
3
0,6
r/R [-]
dFL = ∆w dq
ϕ max = arctan
0,4
Using the result of (3.27), (3.28) and (3.33) in (3.29) we get
From d(dP)/dφ = 0, it follows
ϕ max
0,2
Figure 3.10: Optimal pitch angel
= r 2ω ρ 2π dr w12 2 sin ϕ {sin (2ϕ1 − 3ϕ )}
2
= ϕ1
3
0,0
From the air foil theory we have
(3.35)
Example 3.2
Let’s compare Betz’ and Schmitz’ formulas for the design of the optimal pitch angle. Assuming X =
7; B = 3; αD =7,0°; CL = 0,88 one gets
dFL = ½ ρ w 2 B c dr C L
⎛ϕ ⎞
2
= ½ ρ w1 B c dr C L cos⎜ 1 ⎟
⎝ 3⎠
[7]
(3.37)
where we have used (3.25) and φ = 2/3φ1.
Combining (3.37) and (3.36) we get
1 16 π r
⎛ϕ ⎞
c(r )Schmitz =
sin 2 ⎜ 1 ⎟
B CL
⎝ 3⎠
[m]
(3.38)
Wind turbines
16
or
c(r )Schmitz
⎛1
⎛ R ⎞⎞
1 16 π r
⎟⎟ ⎟
=
sin 2 ⎜⎜ arctan⎜⎜
⎟
B CL
3
⎝ X r ⎠⎠
⎝
[m]
(3.39)
Wind turbines
QM* = ½ ρ w 2 c 2 C M
[Nm]
17
(4.3)
The density of air is at a nominal state, defined as 1 bar and 11°C, 1,225 kg/m3.
The curves in figure 4.1 are given at different Reynolds’s number, defined as
Example 3.3
Let’s again compare Betz’ and Schmitz’ formulas for the design of the optimal pitch angle.
Assuming X = 7; B = 3; αD =7,0°; CL = 0,88 one gets
Optimal chord ratio
X = 5; B = 3; alfa_D = 7,0°; C_L = 0,88
c/R [-]
0,7
0,6
c/R(Betz)
0,5
c/R(Schmitz)
0,3
0,2
0,1
0,0
0,2
0,4
0,6
0,8
[-]
µ/ρ
(4.4)
For PC-calculation it is convenient to have the curves as functions. For the NACA 23012 profile
one can use the following approximation: CD,L = k0 + k1α + k2α2 + k3α3 + k4α4, with the following
constants
NACA 23012
CL
CD
1,0318e-1
6,0387e-3
k0
1,0516e-1
–3,6282e-4
k1
k2
1,0483e-3
5,4269e-5
7,3487e-6
6,5341e-6
k3
–6,5827e-6
–2,8045e-7
k4
Table 4.1: Polynomial constants – for 0 < α < 16°
0,4
0,0
cw
Re =
As shown in figure 4.2, the data are given in the range of α < 20°. For wind turbines it is necessary
to know the data for the range up to 90°. In the range from αst <α < 90° we can use the following
assumptions, see ref./2/
1,0
r/R [-]
Figure 3.11: Optimal chord length
Lift:
Note, near the tip there are no difference between Betz’ and Schmitz’ theory.
C L = A1 sin (2 α ) + A2
cos 2 (α )
sin (α )
[-]
(4.5)
where
Wing profiles are often tested in wind tunnels. Results are curves for coefficient of lift and drag and
moment. Data for a lot of profiles can be found in ”Theory of Wing Sections, Ira H. Abbott and A.
E. Doenhoff, ref./3/.
A1 =
B1
2
( B1
see (4.8) below!)
A2 = (C Ls − C D ,max sin (α st )cos(α st ))
sin (α st )
cos 2 (α st )
[-]
(4.6)
Drag:
Figure 4.1 shows data for the profile NACA 23012.
C D = B1 sin 2 (α ) + B2 cos(α ) + C Ds
Lift, drag and torque (per meter blade width) are defined by the equations
[-]
FL* = ½ ρ w 2 c C L
[N]
(4.1)
where CDs is the coefficient of drag at the beginning of stall, αstall, and
B1 = C D,max
FD* = ½ ρ w 2 c C D
[N]
(4.2)
B2 =
(
)
1
C Ds − C D,max sin 2 (α st )
cos(α st )
[-]
(4.7)
(4.8)
Wind turbines
18
Wind turbines
CDmax can be set at 1. For the NACA 23012 profile, the angle of stall is a little uncertain, but could
in practice be set at 16°. Figure 3.2 show the result of the formulas above.
Figure 4.1 show some typical data for an air foil.
Location of
max. thickness
Max thickness
edge
edge
Uppe
r surf
ace
Mean camber line
Max camber
Location of
max. camber
rface
Lower su
Chord line
Trailing
edge
Chord
Figure 4.1: Definition of typical air foil data
•
•
•
•
The chord line is a straight line connecting the leading and the trailing edges of the air foil.
The mean camber line is a line drawn halfway between the upper and the lower surfaces. The
chord line connects the ends of the mean camber lines.
The frontal surface of the airfoil is defined by the shape of a circle with the leading edge radius
The center of the circle is defined by the leading edge radius and a line with a given slope of the
Data for the NACA 23012 profile is given by the table (upper left corner) on figure 4.2.
Figure 4.2: Data for NACA 23012 (Ref./3/)
19
Wind turbines
20
Wind turbines
21
Here we have used
5. The blade element momentum (BEM) theory
w=
In the blade element momentum (BEM) method the flow area swept by the rotor is divided into a
number of concentric ring elements. The rings are considered separately under the assumption that
there is no radial interference between the flows in one ring to the two neighbouring rings.
or
Figure 3.3 shows the profile and the wind speeds in one ring. The angle of attack α is given by
w=
α =ϕ − β
(5.1)
From figure 3.3 we get
tan (ϕ ) =
1 − a v1
1 + a' r ω
[-]
(5.2)
dT = ½ ρ w 2 c B C y dr
and
dU = ½ ρ w c B C x r dr
2
a' =
[Nm]
(5.4)
where Cy and Cx are given by (3.15) and (3.13)
dT = 2 π r ρ v2 (v1 − v3 )dr
[N]
(5.5)
dU = 2π r 2 ρ v2 u 3 dr
[Nm]
(5.6)
In (5.6) we are using u3 for the tangential speed far behind the rotor plane, even though there is
some tangential rotation of the wind. This can be shown to be an allowable approximation, because
the rotation of the wind normally is small.
a=
[m/s]
(5.10)
[-]
(5.11)
1
4 sin 2 (ϕ )
+1
σ Cy
1
4 sin (ϕ )cos(ϕ )
−1
σ Cx
[-]
(5.12)
[-]
(5.13)
1
4 F sin 2 (ϕ )
+1
σ Cy
[-]
(5.14)
[-]
(5.15)
[N]
(5.16)
and
a' =
Combining (5.3) and (5.5) - and - (5.4) and (5.6) we get
c B Cx
a'
=
a '+1 8π r sin (Φ ) cos(Φ )
ω r (1 + a')
cos(ϕ )
For rotors with few blades it can be shown that a better approximation of a and a’ is
If we now use the laws of momentum and angular momentum, we get
c B Cy
a
=
a − 1 8π r sin 2 (Φ )
(5.9)
and solve the equation (5.7) and (5.9) we get
a=
(5.3)
[m/s]
If we now define the solid ratio as
cB
σ =
2π r
If the number of blades is B, we can calculate the axial force dT and the torque dU on a ring element
with the radius r and the width dr and the torque as
[N]
v1 (1 − a )
sin (ϕ )
1
4 F sin (ϕ )cos(ϕ )
−1
σ Cx
where
[-]
(5.7)
F=
[-]
(5.8)
2
π
⎛
⎛ B R − r ⎞⎞
⎟⎟ ⎟⎟
arccos⎜⎜ exp⎜⎜ −
⎝ 2 r sin (ϕ ) ⎠ ⎠
⎝
Wind turbines
22
This simple momentum theory breaks down when a becomes greater than ac = 0,2. In that case we
replace (5.14) by
a = ½⎛⎜ 2 + K (1 − 2 a c ) −
⎝
(K (1 − 2ac ) + 2)2 + 4(K ac 2 − 1) ⎞⎟
⎠
Wind turbines
P = ω B ∫ r U * (r ) dr
[N]
R
0
[-]
23
(5.22)
(5.17)
6. Efficiency of the wind turbine
where
K=
4 F sin 2 (ϕ )
σ Cy
[-]
(5.18)
6.1. Rotor
Betz has shown that the maximum power available in the wind is given by (2.16). Let us define this
power as
Calculation procedure
16 1
ρ v13 A
27 2
We can now calculate the axial force and power of one ring element of the rotor by making the
following iteration:
Pmax =
For every radius r (4 to 8 elements are OK), go through step-1 to step-8
where we have used Cp = Cp,Betz =16/27.
Step-1: Start
In (6.1) A is the swept area of the rotor, and in the following we define this area as A = π/4 D2 i.e.
we do not take into account, that some part of the hub area is not producing any power!
Step-2: a and a’ are set at some guessed values. a = a’ = 0 is a good first time guess.
(6.1)
We can now define the rotor efficiency as
Step-3: φ is calculated from (5.2)
Step-4: From the blade profile data sheet (or the polynomial approximation) we find CL and CD
[W]
η rotor =
Protor
Pmax
[-]
(6.2)
Step-5: Cx and Cy are calculated by (3.13) and (3.15)
where Protor is the power in the rotor shaft.
Step-6: a and a’ are calculated by (5.14) and (5.15). Or if a > 0,2 then a is calculated from (5.17).
The rotor efficiency can be calculated on the basis of a BEM-calculation of the power production in
a real turbine – see the example in Chapter 7.
Step-7: If a and a’ as found under step-5 differ more than 1% from the last/initial guess, continue at
step-2, using the new a and a’.
Step-8: Stop
The rotor efficiency is divided into three parts
When the iterative process is ended for all blade elements, then the axial force and tangential force
U (r ) = ½ ρ w c C x
[N]
(5.19)
T * (r ) = ½ ρ w 2 c C y
[N]
(5.20)
*
2
and then the total axial force and power as
T = B ∫ T * (r ) dr
R
0
Another model will be presented here:
[N]
(5.21)
η rotor = η wakeη tipη profile
[-]
(6.3)
where “wake” indicates the loss because of rotation of the wake, “tip” the tip loss and “profile” the
profile losses.
Wake loss:
The wake loss can be calculated on the basis of Schmitz’ theory. Integrating (3.31) over the whole
blade area and using (3.8) and (3.33) gives.
Wind turbines
⎛2 ⎞
sin 3 ⎜ ϕ1 ⎟
π 2 3
⎛r⎞
⎝ 3 ⎠ d⎛ r ⎞
= ½ ρ D v1 ∫ 4 X ⎜ ⎟
⎜ ⎟
2
4
⎝ R ⎠ sin (ϕ1 ) ⎝ R ⎠
0
2
1
PSchmitz
24
[W]
(6.4)
This can be solved numerically, see an example in Appendix D. Based on this we can define
C p,Schmitz
P
= Schmitz
1
ρ v13 A
2
[-]
25
Profile loss:
From the power calculation after (3.12) and (3.13) we can see, that the power is proportional to Cx.
For an ideal profile, i.e. with no drag, the power would the be higher, from which we can define the
profile efficiency to
η profile (r ) =
(6.5)
C L cos(γ ) − C D sin (γ )
C
= 1 − D tan (γ )
C L cos(γ )
CL
[-]
(6.8)
Using (3.7) we get
3r X
[-]
(6.8)
2 R GR
Assuming the angle of attack to be the same over the entire blade length the glide ratio is constant
too and then (6.8) can be integrated over the blade length to give
η profile (r ) = 1 −
0,7
0,6
0,5
Cp [-]
Wind turbines
Cp(Betz)
0,4
0,3
η profile = 1 −
Cp(Schmitz)
0,2
0,1
0
0
2
4
6
8
X
GR
[-]
(6.9)
Example 6.1
Assuming the glide ration to be GR = 100 and the blade number to B = 3 then the rotor efficiency
can be calculated as function of the tip speed ratio, see figure 6.2.
10
X [-]
Figure 6.1. Coef. of power according to Betz and Schmitz
Rotor efficiency
The difference between Betz and Schmitz is, that Schmitz takes the swirl loss into account and
therefore we can define swirl loss or the wake loss as
1,0
0,9
C p,Schmitz
[-]
C p,Betz
0,8
(6.6)
Tip loss:
In operation there will be a high negative (compared to ambient) pressure above the blade and a
(little) positive pressure under the blade. Near the tip of the blade, this pressure difference will
induce a by pass flow from the high pressure side to low pressure side – over the tip end of the
blade – thus reducing the differential pressure and then the power production!
profile
0,7
eta_rotor [-]
η wake =
Based on: GR = 100; B = 3
w ake
0,6
tip
0,5
rotor
0,4
0,3
0,2
0,1
0,0
0
2
4
The model of Betz – see ref. /4/, page 153-155 – results in a tip efficiency of
⎛
η tip = ⎜⎜1 −
⎝
⎞
⎟
⎟
2
B X + 4/9 ⎠
0,92
2
[-]
6
8
10
X [-]
(6.7)
Figure 6.2. Rotor efficiency
Most modern wind turbines have tip speed ration at nominal wind speed and power around x = 7,
and from the curve it is obvious, that this is close to optimal!
Wind turbines
26
Wind turbines
Example 6.2
Most modern wind turbines have glide ratios around 100 and three blades. Figure 6.3 shows the
rotor efficiency for 2,3 and 4 blades and with the glide ratio as parameter.
27
Generator
Protor
Pmax
PLS
Gear
box
Pgen
Converter
Pgrid
Figure 6.4: Main components in a wind turbine
The total efficiency of such a turbine can the be defined as
Pgrid
η total =
= η rotorη gearboxη genη conv
[-]
Pmax
Figure 6.3. Rotor efficiency
For X = 7 and for a glide ratio GR = 100 it can bee seen, that the number of blades have the
following influence on the rotor efficiency
3 and 4 blades are more efficient than 2 blades, but also more expensive. When a 3 blade rotor in
spite of that has become a de facto standard it is due to a more dynamical stable rotor.
6.2. Gear box, generator and converter
Most wind turbines have the following main parts, a rotor, a gear box a generator and an electric
converter, see figure 6.4. Each of these components has losses.
(6.10)
where
=
Protor
Pmax
η gearbox =
PHS
Protor
η rotor
η gen
η conv
=
=
Pgen
[-]
(6.11)
PHS
Pgrid
Pgen
where the indices stand for “LS” = low speed (shaft); “gen” = generator; “conv” = frequency
converter and “grid” = grid net.
Typical values for the efficiencies are – at nominal power
Gearbox:
0,95-0,98
Generator: 0,95-0,97
Converter: 0,96-0,98
At part load, the lower values can be expected.
Cooling:
The cooling of the components can be calculated as “power input minus power output”. As an
example for the gear box: Φgearbox = Protor - PLS.
Wind turbines
28
7. Example, BEM
The programming can be done in a spread sheet (Excel). Let us look at an example.
We would like to design a wind turbine with 3 blades, a nominal tip speed ratio of X = 5 and to use
the NACA profile 23012. For this profile we have an optimum glide number (=CL/CD) at an angle
of attack on 7,0° and here the coefficient of lift is 0,88. The result is shown in figure 7.1.
Next we want to calculate the power and axial force on the turbine with a radius of 5 m at a wind
speed of 10 m/s and at a rotary speed of 88 rpm. The calculation is shown in figure 7.2.
Figure 7.3 shows the performance of the turbine at a fixed number of revolutions but varying wind
speed.
B
C
D
E
F
G
H
I
J
K
Design procedure for an optimal rotor
2008.07.05/SGt
Blue numbers are input to the program
Outputs are teta and c/R, see figures
Tip speed ratio
Angle of attack
Coeff. of lift
Table 1: Input
X
B
alpha
C_L
°
-
5
3
7
0,88
r/R
Speed ratio
x
Angle, optimal
phi
°
Pitch
beta
°
Rel. chord length
c/R
Table 2: Output (according to Schmitz)
8 element-model where there is no blade
at the inner element!
0,188
0,94
31,2
24,2
0,259
0,313
1,56
21,7
14,7
0,212
0,438
2,19
16,4
9,4
0,169
0,563
2,81
13,0
6,0
0,138
Optimal pitch angle and chord ratio
60
0,3
beta
50
c/R
40
0,2
c/R [-]
beta [°]
A
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
30
20
0,1
10
0
0,0
0,0
0,2
0,4
0,6
0,8
1,0
r/R [-]
Figure 7.1: Design of the rotor (Formulas, see App. B)
0,688
3,44
10,8
3,8
0,116
0,813
4,06
9,2
2,2
0,100
0,938
4,69
8,0
1,0
0,087
Wind turbines
1
A
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
B
B
C
C
D
D
E
E
F
F
G
G
29
H
H
I
I
J
J
Power calculation
2008.07.05/SGt
R_o
Wind speed
v_1
Rotational speed
n
Density of air
rho
B
Angular speed
omega
Thickness, 1 ring
dr
R_i
Swept surface
A_s
Max. power
P_max
Tip speed
V_TIP
Tip speed ratio
X_act
Table 1: Rotor size and speed
m
m/s
min^-1
kg/m3
s^-1
m
m
m2
W
m/s
-
5
10
88
1,225
3
9,2
0,625
0,625
78,54
28507
46,1
4,61
Ring no.
Pitch
Chord
Table 3: From "Design"
N_r
r/R
r
beta
c
m
°
m
1
0,19
0,94
24,232
1,293
Pitch, no pitch control
Chord
Pitch angle
Solid ratio
beta_0
c
beta
sigma
r*omega
°
m
°
m/s
Axial int.factor
Tang. int. factor
Angle of rel. wind
Angle of attack
Coeff. of lift
Coeff. of drag
y-component
x-component
Factor
Factor
Axial int.factor (1)
Axial int.factor (2)
Axial int.factor
Tang. int. factor
a and a' must be
guessed until errors
< 2%
Rel. speed
Tang. force
Axial force
Power
Swept area
Axial force
Table 4: Calculation
Wing performance
k0
0,10318
k1
0,10516
k2
0,001048
k3
7,35E-06
k4
-6,58E-06
alfa_s
16
C_ls
1,652801
C_ds
2,25E-02
B1
1,00
B2
-0,0556
A1
0,500
A2
0,414
Table 2: Wing data
2
0,31
1,56
14,746
1,059
3
0,44
2,19
9,378
0,845
4
0,56
2,81
6,049
0,691
5
0,69
3,44
3,813
0,581
6
0,81
4,06
2,219
0,500
If you want to use the pitch and cord values from the optimal rotor design y
copy the values in rows 24 and 25 to rows 29 and 30!
a
a'
PHI
alpha
C_L
C_d
C_y
C_x
F
K
a_1
a_2
a
a'
error_a
error_a'
°
°
-
%
%
24,232
1,293
24,2
0,66
8,6
14,746
1,059
14,7
0,32
14,4
9,378
0,845
9,4
0,18
20,2
6,049
0,691
6,0
0,12
25,9
3,813
0,581
3,8
0,08
31,7
2,219
0,500
2,2
0,06
37,4
Use the macro "Ctrl + p" to made the iteration ! You might have to use th
more times to reduce the errors in rows 51 and 52 to 0
0,316
0,243
32,5
8,2
1,0154
0,0091
0,861
0,538
1,000
2,032
0,330
0,316
0,316
0,243
0
0
0,319
0,099
23,3
8,5
1,0460
0,0095
0,965
0,405
1,000
2,001
0,333
0,319
0,319
0,099
0
0
0,315
0,052
17,9
8,5
1,0461
0,0095
0,998
0,313
0,999
2,051
0,328
0,315
0,315
0,052
0
0
0,310
0,031
14,5
8,4
1,0348
0,0093
1,004
0,250
0,994
2,107
0,322
0,310
0,310
0,031
0
0
0,309
0,021
12,1
8,3
1,0159
0,0091
0,995
0,203
0,976
2,122
0,320
0,309
0,309
0,021
0
0
0,320
0,015
10,1
7,9
0,9794
0,0087
0,966
0,164
0,911
1,990
0,334
0,320
0,320
0,015
0
0
w
F_x
F_y
m/s
N/m
N/m
12,7
69,0
110,6
17,2
77,9
185,6
22,3
80,4
256,6
27,6
80,6
324,2
33,1
79,2
387,5
38,6
74,8
440,7
P_r
A_i
A_i
F_yi
W
m2
%
N
1118
3,682
5
207,3
2102
6,136
8
348,1
3037
8,590
11
481,1
3915
11,045
14
607,8
4705
13,499
17
726,6
5248
15,953
20
826,3
°
m/s
min^-1
kW
%
Nm
N
°
0
10
88
24,85
87,17
2,697
4039
4,61
8,1
Pitch control
dbeta
Wind speed
v_1
Rotational speed
n
Power
P
Efficiency
eta_r
Torque
M
Axial force
T
Tip speed ratio
X_act
Mean angle of attack alpha_m
Table 5: Summary of results
Figure 7.2: Calculation of power and axial force (Formulas, see App. B)
5
10
15
20
25
30
0
Wind speed [m /s]
2
4
6
8
10
Tip speed ratio [-]
Figure 7.3: Power as function of wind speed (left) and efficiency as function of tip speed ratio
(right)
24,5
0
0
22,5
10
0
20,5
20
10
18,5
30
20
16,5
40
14,5
50
12,5
30
60
10,5
40
70
8,5
50
6,5
60
31
Distribution of wind A=8 & k = 2
0,14
0,12
0,10
0,08
0,06
0,04
0,02
0,00
4,5
70
80
2,5
90
Wind turbines
0,5
100
80
30
p(v(i)<v<v(i+1))
90
Efficiency [%]
Power [kW]
Wind turbines
v_m [m/s]
Figure 8.1: Distribution of the natural wind – for ∆v = 1 m/s
Production curve: Let us imagine a power curve for a given stall controlled wind turbine given by
P(v) = min ( kP v3, PN )
[W]
(8.3)
with kP = 0,1 kW/(m/s)3 and PN = 200 kW.
8. Distribution of wind and annual energy production
The nominal power 200 kW would be reached at a wind speed of 12,6 m/s.
Weibull distribution
This would result in a power curve as given by figure 8.2.
The wind is distributed close to the Weibull distribution curve. For practical purposes one can
calculate the probability for the wind being in the interval vi < v < vi+1
Power curve
250
(8.1)
200
P(v) [kW]
⎛ ⎡⎛ v ⎞ k ⎤ ⎞
⎛ ⎡⎛ v ⎞ k ⎤ ⎞
p (vi < v < vi +1 ) = exp⎜ − ⎢⎜ i ⎟ ⎥ ⎟ − exp⎜ − ⎢⎜ i +1 ⎟ ⎥ ⎟ [-]
⎜ ⎢⎝ A ⎠ ⎥ ⎟
⎜ ⎢⎝ A ⎠ ⎥ ⎟
⎦⎠
⎦⎠
⎝ ⎣
⎝ ⎣
where A and k are found for a given site on the basis of measurements.
150
100
50
0
Annual production
0
5
10
15
20
25
30
Wind speed, v [m/s]
If the power for the turbine at a given wind speed is P(u), the annual production can be calculated as
E ann = ∑ {8766h ⋅ p (vi < v < vi +1 ) P(v m )}
Figure 8.2. Power curve
[J]
(8.2)
where vm is the mean value of vi and vi+1 i.e. vm = ( vi + vi+1)/2.
Example
Typical values of A and k could be A = 8 m/s and k = 2. This will give the following distribution
If we combine the data from figure 8.1 and 8.2, we will get the energy production curve as shown in
figure 8.3.
Wind turbines
32
24,5
22,5
20,5
18,5
16,5
14,5
12,5
8,5
10,5
6,5
4,5
2,5
30.000
20.000
10.000
0
0,5
E [kWh]
70.000
60.000
50.000
40.000
v_m [m/s]
Figure 8.3: Annual distribution
What would be an optimal maximum power, Pmax ? From figure 8.1 we can see that wind speed
above 15 – 20 m/s are very rare. On the contrary, the power production of af wind turbine rises with
a power of 3.
Figure 8.4 shows a calculation of the annual production as a function of the maximum power, Pmax.
200
400
600
0,7
0,6
0,5
0,4
E_ann
0,3
CF
0,2
0,1
0,0
800
1000
CF [-]
E_ann [MWh]
Annual energy production
0
9. Symbols
Annual production
700
600
500
400
300
200
100
0
Wind turbines
P_N [kW]
Figure 8.4: Annual energy production and capacity factor as function of nominal power
The capacity factor is definded as CF = Eann / (PN 8766h).
Conclusion:
To answer the question we must know the price of the turbine, including tower and foundations, but
more than about 200-300 kW does not seems reasonable.
a
a’
A
A
A1
A2
B1
B2
B
CD
CD,max
CDst
CL
CLst
CP
CF
Cy
Cx
c
Eann
F
FL
FD
k
K
M
n
p
ptot
P
PN
Pmax
QM *
r
Re
T
T*
U
U*
u2 = u
v2 = v
v1
v3
vTIP
w
X
m/s
m2
m
Wh
N/m
N/m
Nm
1/s
Pa
Pa
W
W
W
N/m
m
N
N
N
N
m/s
m/s
m/s
m/s
m/s
m/s
-
Axial interference factor
Tangential interference factor
Wind speed, distribution curve
Area, swept area of the rotor
Constant
Constant
Constant
Constant
Coefficient of drag
Coefficient of drag, max value, at α = 90°
Coefficient of drag, where stall begins
Coefficient of lift
Coefficient of lift, where stall begins
Power production factor
Axial force factor
Coefficient of axial forces
Coefficient of tangential forces
Chord length
Annually produced energy
Calculation value
Lift force (per length of blade)
Drag force (per length of blade)
Constant
Factor
Torque
Rotational speed of rotor
Pressure
Total pressure (Bernouilli’s equation)
Power
Power, nominal wind
Max power of a given turbine
Reynold’s number
Axial force (thrust) on the rotor
Axial force per width of the blades
Tangential force on the rotor
Tangential force per width of the blades
Tangential wind speed in the rotor plane
Axial wind speed in the rotor plane
Wind speed, upstream the rotor
Wind speed, down-stream the rotor
Relative wind speed
Tip speed ratio
33
Wind turbines
x
-
Local speed ratio
αA
αst
°
°
°
°
°
°
s-1
kg/(m s)
Pa
m/s
m/s
m/s
kg/m3
Angle of attack, relative wind in relation to blade chord
Angle of attack, where stall begins
Pitch angle of the blade to rotor plane
Relative wind to rotor axis
Efficiency
Angle of relative wind to rotor plane
Angular velocity
Dynamic viscosity
Differential pressure, over the rotor
Change of relative wind speed
Change of tangential wind speed
Change of wind speed
Density of air (here 1,225 kg/m3)
β
γ
η
φ
ω
µ
∆p
∆w
∆u
∆v
ρ
34
Wind turbines
35
App. A: Conservation of momentum and angular momentum
Momentum
Momentum of a particle in a given direction is defined as
p = mu
(A1)
where m is mass and u is speed of the particle
According to the Newton’s 2nd law we have
F=
dp
dt
(A2)
where F is the force acting on the particle
10. Literature
/1/
/2/
/3/
/4/
If the mass is constant, we have (Newton’s 2nd law)
Andersen, P. S. et al
Basismateriale for beregning af propelvindmøller
Forsøgsanlæg Risø, Risø-M-2153, Februar 1979
Guidelines for design of wind turbines
Wind Energy Department, Risø, 2002, 2nd edition
ISBN 87-550-2870-5
F =m
(A3)
where a is the acceleration of the particle
If we have a flow of particles with the mass flow qm we can calculate the force to change the
velocity for u1 to u2 as
Abbott, I. H., Doenhoff, A. E.
Theory of wing sections
Dover Publications, Inc., New York, 1959
Gasch, R; Twele, J.
Wind power plants - Fundamentals, Design, Construction and Operation
James and James, October 2005
du
=ma
dt
F = q m (u 2 − u1 )
(A4)
Force equals differential pressure, ∆p, times area, A, i.e. (A4) can be written as
∆p =
q m (u 2 − u1 )
A
(A5)
Example
For a wind turbine we have a wind speed up-stream the turbine of u1 = 8 m/s and a wind speed
down stream of u2 = 2,28 m/s. In the rotor plane the wind speed is just the mean value of these to
values, i.e. u = 5,14 m/s. The blade length is R = 25 m. Find the axial force on the rotor and the
differential pressure over the rotor.
First we calculate the mass flow as
q m = qV ρ = π R 2 u ρ = π ⋅ 25 2 ⋅ 5,14 ⋅1,225 = 12369 kg/s
Wind turbines
36
Wind turbines
37
and
Using (A4) we get F = 12369 ( 2,28 – 8,0 ) = -70,7 kN. The negative sign tells us that the force is in
the opposite direction to the flow. The differential pressure is calculated by (5) giving ∆p = 36 Pa.
M 2 = m r22
ω2
(A9)
t
The torque to produce the change of angular momentum can then be calculated as
Angular momentum:
M = M 2 − M1 =
ut
(
m
ω 2 r22 − ω1 r12
t
)
(A10)
This formula applies equally to a stream of fluid moving in a curved path, since m/t is the mass
flowing per unit of time, qm. Thus the torque which must be acting on a fluid will be
m
r
M = q m (ω 2 r22 − ω1 r12 )
(A11)
M = q m (u 2 t r2 − u1 t r1 )
(A12)
or
Figure A1: Rotating mass
Example
Figure A1 shows a particle of mass m rotating a radius r with a tangential velocity of ut. The
angular momentum, L, is given by
Figure 2 shows a wind turbine with 2 blades. The blade length is R = 25 m and the rotational speed
is n = 25 rpm which gives an angular velocity of ω = 2,62 s-1.
L = m r ut = m r 2
ut
= mr2 ω = Iω
r
(A6)
R
where I is the moment of inertia and ω is the angular velocity.
The torque of the particle is given by
dL
dω
M =
=I
= Iα
dt
dt
u 2t
dr
r
(A7)
where α is the angular acceleration
Now consider a particle moving in a curved path, so that in time t it moves from a position at which
it has an angular velocity ω1 at radius r1 to a position in which the corresponding values are ω2 and
r2. To make this change we must first apply a torque, M1, to reduce the particle’s original angular
momentum to zero, and then apply a torque, M2, in the opposite direction to produce the angular
momentum required in the second position, i.e.
M 1 = m r12
ω1
t
(A8)
u1
u
u2
Figure A2: Wind turbine with 2 blades
Let us calculate the power for the annular element given by radius r = 17 m and with a thickness of
dr = 10 cm. In a calculation concerning the BEM theory, one can find the axial velocity in the rotor
plane at u = 5,14 m/s (a = 0,357) and at tangential velocity of the air after pasing the rotor plane at
u2t = 0,65 m/s (a’ = 0,0072)
Wind turbines
38
The mass flow through the annular element is
(
)
App. B: Formulas, spread sheet calculations
(
)
q m = qV ρ = π (r + dr ) − r 2 u ρ = π (17 + 0,1) − 17 2 ⋅ 5,14 ⋅1,225 = 68,0 kg/s
2
Wind turbines
2
In formula (A12) we have u1t = 0 because there is no rotation of the air before the rotor plane and
u2t = 0,65 m/s and r1 = r2 = r = 17 m. The torque can be calculated at
M = 68,0 (0,65 ⋅ 17,0 − 0 ⋅ 0 ) = 755 Nm
The power can be calculated at P = M ω = 1,98 kW
Formulas in spread sheet, see figure 7.1
14
15
16
17
18
E
=1,5/8
=E14*\$E\$7
=2/3*ATAN(1/\$E\$7/E14)*180/PI()
=E16-\$E\$9
=1/\$E\$8*16*PI()*E14/\$E\$10*(SIN(1/3*ATAN(1/\$E\$7/E14)))^2
Formulas in spread sheet, see figure 7.2
E
22
23
24
25
=Design!E14
=E22*\$E\$6
=Design!E17
=Design!E18*\$E\$6
29
30
31
32
33
24,2317401773297
1,29343334743683
=E29+\$E\$64
=MIN(E30*\$E\$10/(2*PI()*E23);1)
=E23*\$E\$11
E
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
E
0,316464512669469
0,24291710357473
=ATAN((1-E37)/(1+E38)*\$E\$7/(E23*\$E\$11))/PI()*180
=E39-E31
=IF(E40<\$J\$11;\$J\$5+\$J\$6*E40+\$J\$7*E40^2+\$J\$8*E40^3+\$J\$9*E40^4;\$J\$16*SIN(2*E40*PI()/180)+\$J\$17*(COS(E40*PI()/180))^2/SIN(E40*PI()/180))
=IF(E40<=\$J\$11;\$K\$5+\$K\$6*E40+\$K\$7*E40^2+\$K\$8*E40^3+\$K\$9*E40^4;\$J\$14*SIN(E40*PI()/180)^2+\$J\$15*COS(E40*PI()/180)+\$J\$13)
=E41*COS(PI()/180*E39)+E42*SIN(PI()/180*E39)
=E41*SIN(PI()/180*E39)-E42*COS(PI()/180*E39)
=2/PI()*ACOS(EXP(-\$E\$10/2*(\$E\$6-E23)/E23/SIN(E39*PI()/180)))
=4*E45*SIN(E39*PI()/180)^2/E32/E43
=1/(E46+1)
=1/2*(2+E46*(1-2*0,2)-SQRT((E46*(1-2*0,2)+2)^2+4*(E46*0,2^2-1)))
=IF(E37>0,2;E48;E47)
=1/(4*E45*SIN(PI()/180*E39)*COS(PI()/180*E39)/E32/E44-1)
=(E49-E37)/E37*100
=(E50-E38)/E38*100
=\$E\$7*(1-E37)/SIN(PI()/180*E39)
=0,5*\$E\$9*E54^2*E30*E44
=0,5*\$E\$9*E54^2*E30*E43
=\$E\$11*\$E\$10*\$E\$12*E55*E23
=PI()*(((E21+1)*\$E\$13)^2-(E21*\$E\$13)^2)
=E59/\$E\$14*100
=\$E\$10*E56*\$E\$13
E
64
65
66
67
68
69
70
71
72
0
10
88
=SUM(E58:K58)/1000
=E67*1000/E15*100
=E67/E11
=SUM(E61:K61)
=E17
=AVERAGE(E40:K40)
39
Wind turbines
40
The iteration to find a and a’ in row 37 and row 38 (see example on figure 10) can be quite tedious.
The job can be done automatically by the following macro. Run the macro by typing “Ctrl + p”.
Sub Makro1()
'
' Makro1 Makro
' Makro indspillet 13-10-03 af Søren Gundtoft
'
' Genvejstast:Ctrl+p
'
Range("E50:K50").Select
Application.CutCopyMode = False
Selection.Copy
Range("E38").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E49:K49").Select
Application.CutCopyMode = False
Selection.Copy
Range("E37").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E50:K50").Select
Application.CutCopyMode = False
Selection.Copy
Range("E38").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E49:K49").Select
Application.CutCopyMode = False
Selection.Copy
Range("E37").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E50:K50").Select
Application.CutCopyMode = False
Selection.Copy
Range("E38").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E49:K49").Select
Application.CutCopyMode = False
Selection.Copy
Range("E37").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E50:K50").Select
Application.CutCopyMode = False
Selection.Copy
Range("E38").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
Range("E49:K49").Select
Application.CutCopyMode = False
Selection.Copy
Range("E37").Select
Selection.PasteSpecial Paste:=xlValues, Operation:=xlNone, SkipBlanks:= _
False, Transpose:=False
End Sub
Wind turbines
App. C: Formulas, EES-programme
"BEM-model
2008.07.07/SGt"
PROCEDURE C_lift_and_drag(alpha:C_lift;C_drag)
"#####################################################################"
"This procedure calculates the coefficient of lift and drag for the NACA23012 profile"
alpha_s=16
C_dmax=1
"lift and drag at stall"
C_lift_s=1,0318E-01+1,0516E-01*alpha_s+1,0483E-03*alpha_s^2+7,3487E-06*alpha_s^3-6,5827E06*alpha_s^4
C_drag_s=6,0387E-03-3,6282E-04*alpha_s+5,4269E-05*alpha_s^2+6,5341E-06*alpha_s^3-2,8045E07*alpha_s^4
"Some constants"
B1=C_dmax
B2=1/COS(alpha_s)*(C_drag_s-C_dmax*SIN(alpha_s)^2)
A1=B1/2
A2=(C_lift_s-C_dmax*SIN(alpha_s)*COS(alpha_s))*(SIN(alpha_s)/(COS(alpha_s))^2)
IF alpha<alpha_s THEN
C_lift=1,0318E-01+1,0516E-01*alpha+1,0483E-03*alpha^2+7,3487E-06*alpha^3-6,5827E-06*alpha^4
C_drag=6,0387E-03-3,6282E-04*alpha+5,4269E-05*alpha^2+6,5341E-06*alpha^3-2,8045E-07*alpha^4
ELSE
C_lift = A1*sin(2*alpha)+A2*(cos(alpha))^2/sin(alpha)
C_drag = B1*(sin(alpha))^2 + B2*cos(alpha) + C_drag_s
ENDIF
END
PROCEDURE interference_factors(V_0;omega;B;R;ri;c;betai: a_axial;a_tangential;C_x;C_y)
"#####################################################################"
"Initialize"
a_axial=0,3
a_tangential=0,000001
REPEAT
"Reset old iteration values"
a_axial_old=a_axial
a_tangential_old=a_tangential
"Calculate angle of attack and solidity"
phi=arctan(((1-a_axial)/(1+a_tangential))*(V_0/(ri*omega))) {Determine rel. angle of attack, [rad]}
alpha=phi-betai
{Determine actual angle of attack correcting for pitch, [°] }
sigma=(c*B)/(2*pi*ri)
{Determine solidity, [-] }
CALL C_lift_and_drag(alpha:C_l;C_d)
C_y=C_l*COS(phi)+C_d*SIN(phi)
C_x=C_l*SIN(phi)-C_d*COS(phi)
"Perform Prandtl tiploss correction"
g=(B/2)*(R-ri)/(ri*sin(phi))
expmg=exp(-g)
F=(2/pi)*(arccos(expmg)*pi/180)
K=4*F*SIN(phi)^2/sigma/C_y
a_1=1/(4*F*SIN(phi)^2/sigma/C_y+1)
hj=max((K*(1-2*0,2)+2)^2+4*(K*0,2^2-1);0,000001)
a_2=1/2*(2+K*(1-2*0,2)-SQRT(hj))
IF a_axial_old>0,2 THEN
a_axialx=a_2
ELSE
a_axialx=a_1
ENDIF
a_tangentialx=1/(4*F*SIN(phi)*COS(phi)/sigma/C_x-1)
"Underdamping"
underdamping_factor=0,5
a_axial=a_axial_old+underdamping_factor*(a_axialx-a_axial_old)
41
Wind turbines
42
Wind turbines
a_tangential=a_tangential_old+0,5*(a_tangentialx-a_tangential_old)
error_a=sqrt(( (a_axial-a_axial_old)/a_axial_old*100 )^2)
error_t=sqrt(((a_tangential-a_tangential_old)/a_tangential_old*100)^2)
epsilon=0,01
UNTIL((error_a<epsilon) AND (error_t<epsilon))
END "Procedure"
"#####################################################################"
"
MAIN PROGRAM
"
"#####################################################################"
R=5
B=3
V_0=10
n=88/60
rho_a=1,225
beta_p=0
n_elements=8 "[-]
"[m]
"[-]
"[m/s] Wind speed"
"[rps] Rotor speed"
"[kg/m3] Density of air"
"[°]
Pitch angle"
Number of ring elements"
omega = 2*pi*n
"[s^-1]
Angular speed"
"The model devides the blande in 8 ring element, with no air foil in the inner element"
DUPLICATE i=1;7
r_div_R[i]=(i+0,5)/8
END
"chord and pitch angle"
c_div_R[1]=0,2889 : beta_0[1]=23,92
c_div_R[2]=0,2249 : beta_0[2]=14,98
c_div_R[3]=0,1787 : beta_0[3]=9,96
c_div_R[4]=0,1474 : beta_0[4]=6,76
c_div_R[5]=0,1252 : beta_0[5]=4,56
c_div_R[6]=0,1088 : beta_0[6]=2,96
c_div_R[7]=0,0961 : beta_0[7]=1,74
DUPLICATE i=1;7
beta[i]= beta_0[i]+beta_p
END
"BEM"
DUPLICATE i=1;7
r[i]=r_div_R[i]*R
c[i]=c_div_R[i]*R
CALL interference_factors(V_0;omega;B;R;r[i];c[i];beta[i]: a_axial[i];a_tangential[i];C_x[i];C_y[i])
w_rel[i]=sqrt((omega*r[i]*(1+a_tangential[i]))^2+((1-a_axial[i])*V_0)^2)
F_x[i]=0,5*rho_a*w_rel[i]^2*c[i]*C_x[i]
F_y[i]=0,5*rho_a*w_rel[i]^2*c[i]*C_y[i]
P[i]=R/8*B*omega*F_x[i]*r[i]
Fa[i]=B*R/8*F_y[i]
END
P_rotor=SUM(P[i];i=1;7)/1000
F_rotor=SUM(Fa[i];i=1;7)
P_r_max=16/27*pi*R^2*1/2*rho_a*V_0^3/1000
eta_r=P_rotor/P_r_max*100
App. D: Formulas, EES-programme, to solve the integral in (6.4)
"Efficiency of a wind turbine rotor"
"08.07.2008/SGt"
FUNCTION fCp_schmitz(X)
"This function calculates the Cp-wake factor
according to Schmitz's theory
09.07.2008/SG"
i_max = 1000
dRR = 1/i_max
i = -1
x_sum = 0
REPEAT
i = i+1
RR = i/i_max
IF (RR=0) THEN
phi = 90
ELSE
phi = arctan(1/(X*RR))
ENDIF
xx = 4*X*RR^2*(sin(2/3*phi))^3/(sin(phi)^2)*dRR
x_sum = x_sum + xx
UNTIL (i=i_max)
fCp_schmitz = x_sum
END
"========= Main ========="
GR = 100
B=3
X=7
"[-] - Glide ratio = C_L/C_D"