Modeling

Vortex Lattice Modelling of Winglets
on Wind Turbine Blades
Mads Døssing
Risø-R-1621(EN)
Risø National Laboratory

Technical University of Denmark
Roskilde, Denmark
August 2007
Author: Mads Døssing Risø-R-1621(EN)
Mek-FM-EP 2007-04
Title: Vortex Lattice Modelling of Winglets on Wind Turbine
August 2007
Blades
Departments: Wind Energy Department - Risø & Department of
Mechanical Engineering - DTU
3rd version. 18 oct. 2007
Abstract: ISSN 0106-2840

The power production of wind turbines can be increased by ISBN 978-87-550-3633-8
the use of winglets without increasing the swept area. This
makes them suitable for sites with restrictions in rotor
diameter and in wind farms. The present project aims at
understanding how winglets influences the flow and the
aerodynamic forces on wind turbine blades. A free wake
vortex lattice code and a fast design algorithm for a horizontal
axis wind turbine under steady conditions has been
developed. 2 winglet designs are threated in detail.
Information Service Department

Risø National Laboratory
Technical University of Denmark
P.O.Box 49
DK-4000 Roskilde
Denmark
Telephone +45 46774004
bibl@risoe.dk
Fax +45 46774013
www.risoe.dk
Summary
The power production of wind turbines can be increased by the use of winglets
without increasing the swept area. This makes them suitable for sites with
restrictions in rotor diameter and in wind farms.
The present project aims at understanding how winglets influences the flow and
the aerodynamic forces on wind turbine blades.
A free wake vortex lattice code and a fast design algorithm for a horizontal axis
wind turbine under steady conditions has been developed.
2 winglet designs are threated in detail.

ii
Resumé
Vindmøllers energi-produktion kan øges ved brug af en winglet, uden at det

bestrøgne areal øges. Dette gør dem egnede til brug på anlæg, hvor der er
restriktioner på rotorens diameter.
Formålet med nærværende projekt er at opnå en forståelse for, hvorledes en

winglet påvirker strømningnen og krafterne på en vindmølle vinge.
En vortex-lattice kode og en hurtig design algoritme er blevet udviklet for en

horisontal mølle under tidsuafhængige forhold.
2 winglet design er behandlet i detaljer.

iv
Preface
This thesis was prepared at the Wind Energy Department at Risø in collabora-
tion with the Department of Mechanical Engineering, the Technical University
of Denmark in partial fulfillment of the requirements for acquiring the M.Sc.
degree in engineering.
The thesis deals with different aspects of mathematical and numerical modeling
of wind turbine aerodynamics. The main focus is on development of methods
for prediction the performance of turbine blades with winglets.
I thank my supervisors Mac Gaunaa and Robert Flemming Mikkelsen for their
aid and support.
Kongens Lyngby, August 2007
Mads Døssing
vi
Contents
Summary i
Resumé iii
Preface v
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Winglet Description . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Overview of the Methods Applied . . . . . . . . . . . . . . . . . . 3
1.4 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Theory 7
2.1 Definition of Coordinate Systems . . . . . . . . . . . . . . . . . . 7
2.2 Definition of Blade and Winglet Geometry . . . . . . . . . . . . . 8
2.3 Dimensionless Variables . . . . . . . . . . . . . . . . . . . . . . . 9

viii CONTENTS
2.4 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.5 The General Solution and Boundary Conditions . . . . . . . . . . 10
2.6 The Kutta Condition . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.7 Singularity Elements . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.8 Kutta-Joukowsky and Helmholz Theorems . . . . . . . . . . . . . 15
2.9 Wake Strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.10 Wake Shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.11 The Unsteady Bernoulli Equation . . . . . . . . . . . . . . . . . . 16
3 Numerical Lifting Line Model 19
3.1 Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Classic Analytical Lifting Line Results . . . . . . . . . . . . . . . 22
3.3 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.4 Effect of Winglet Height and Curve Radius . . . . . . . . . . . . 24
3.5 Optimum Circulation distribution on Wings . . . . . . . . . . . . 25
3.6 Generalization of Wing Results to Turbines . . . . . . . . . . . . 27
3.7 Force and Viscous Drag Calculation Using Lifting Line Data . . 27
3.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
4 A Freewake, Vortex Lattice and Panel Method 31
4.1 Freewake Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Evaluation of Forces . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.3 Calculation of Wing Geometries for Vortex Lattice Simulations . 39

CONTENTS ix
4.4 Grid Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.5 Reference Blade Data . . . . . . . . . . . . . . . . . . . . . . . . 41
4.6 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5 Blade and Winglet Design 49
5.1 Defining Bound Circulation . . . . . . . . . . . . . . . . . . . . . 50
5.2 Calculating Velocities . . . . . . . . . . . . . . . . . . . . . . . . 52
5.3 Summary of Method . . . . . . . . . . . . . . . . . . . . . . . . . 56
5.4 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.5 General Design Results . . . . . . . . . . . . . . . . . . . . . . . . 60
5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6 Final Designs 63
6.1 Conservative Winglet Design . . . . . . . . . . . . . . . . . . . . 63
6.2 Optimum Winglet Design . . . . . . . . . . . . . . . . . . . . . . 68
6.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
7 Conclusion 73
A Appendix 75
A.1 Constant Strength Source Panel (3 Nodes) . . . . . . . . . . . . . 75

x CONTENTS
Part 1
Introduction
The present report describes the mathematical methods and the results obtained
regarding winglets mounted on windturbines. The main emphasis is on the
aerodynamics and the associated forces, thus neglecting any elastic behavior of
the blades. Efficient mathematical models for simulation of winglets on turbines
were not available at the beginning of this project and the original idea was
to make Vortex Lattice simulations in order to study the physical properties.
Severe shortcomings of the method soon appeared and this project evolved into
a study of numerical methods and the general properties of winglets. Eventually
a modified Vortex Lattice method was used for simulation and a nearwake design
method was developed which is considered more accurate.
1.1 Motivation
Winglets increases the efficiency of wind turbine rotors by increasing the power
production for a given blade length, making them attractive on sites with design
restrictions on the rotor radius. On sites without restrictions the easiest way to
increase the power is by using a longer traditional flat blade, but in wind farms
this may have a negative impact on the total production due to the increased
sweept area. Therefore, even if there are no design restrictions, winglets may still
be attractive. Since a blade with a winglet is more expensive than a flat blade
2 Introduction
and the manufacturing process is not well established, a clear understanding of

the aeroelastic properties is needed in order for the industry to verify whether
a production is profitable or not.
1.2 Winglet Description
Winglets are essentially extensions of the main wing. The cant, sweep and height
describes the overall shape of a straight winglet, and the local shape is described
by the chord, twist angle and 2D profile. The bend has a curve radius which
is also important. The winglets threated in this work has zero sweep and cant
angle, leaving the curve radius of the bend, the height and the twist and chord
distribution as geometrical parameters. The chord and twist distribution are
very important for the winglet efficiency and determining them is the aim of the
design process. The winglet height can be considered a design parameter. The
efficiency increases with height until the viscous drag outbalances the benefits
of the winglet and usually there is an optimum winglet height around 3-4% of
the blade length.
On airplanes the winglet works by decreasing the induced drag on the main
wing. The induced drag is the component of lift in the direction opposite to the
movement and despite its name it is therefore not associated with any viscous
effects but instead to the change in inflow angle. On airplanes the inflow is
angled down thereby tilting the force vector back. This change in the local
inflow is caused by the complicated 3 dimensional nature of the wings where
pressure variations redistributes the flow. On airplanes the change in flow is a
velocity component, the downwash, perpendicular to the wing in the direction
from the suction to the pressure side. One very visible effect is the tip vortex
created by the pressure difference on the pressure and suction side of the wing
which accelerates the air around the tip in a direction perpendicular to the wing.
Traditionally the winglet has been seen as a device for decreasing this tip vortex.
The principle is that by extending the wing the tip vortex is moved further away
from the main blade and the strength of it is decreased, hence decreasing the
downwash on the center parts of the wing where the aerodynamic forces are
large. The wing can also be extended simply by increasing the span and this is
actually more efficient aerodynamically since the force on the extra span adds
to the lift. The downwash depends not only on the shape of the wing planform,
the force distribution over the blade and winglet is also important.
The physics near wind turbine tips is essentially the same; a tip vortex is trailing
off and generates a downwash on the blade. This tilts the local force vector in
such a way that a force component, acting against the blade direction of rotation,
1.3 Overview of the Methods Applied 3
is created. On the main blade the direction of this unwanted downwash is

against the incoming wind and this confuses matters because a well designed
wind turbine already has a significant induction in this direction. I.e. 1/3 of
the incoming wind under optimum conditions according to the theoretical Betz
limit. However, even if the induction is less than desired the downwash will
not contribute positively because a relatively small winglet creates only local
effects, and will not alter the wake far downstream upon which the momentum
theory and the Betz limit is based. This postulate, that the winglet acts locally
is very important and will be threated in detail under the development of the
design method. On the winglet itself the direction of the unwanted downwash is
radial outward, and more generally the unwanted velocity component is termed
normalwash. It is possible to generate a negative normalwash on the winglet,
in which case a force component is created in the direction of rotation but
the viscous drag may outbalance this. If the induced drag can be reduced it
can have a relatively large effect on the torque generated by the blade because
the reduction is mainly on the outer part. This is especially true if the blade
is heavily loaded near the tip. It is therefore concluded that a reduction of
induced drag is a primary concern. When the winglet is mounted the forces
near the tip becomes larger and this in itself increases the power production
and the thrust. Therefore, the winglet affects not only the induced drag but
also the magnitude of all other aerodynamic forces near the tip. A design should
therefore be evaluated on the power (Cp ) and thrust coefficients (CT ) and not
just the induced drag. An increase in Cp can be caused mainly by heavier
loading near the tip, and not as much by a decrease in induced drag. This is
associated with a large increase in CT , and since the majority of thrust forces are
concentrated near the tip the bending moments also increases substantially. In
the winglet sense, it is better to design for a large reduction in induced drag with
a small associated increase in thrust forces. The forces on the winglet itself does
not contribute to the thrust. Besides from a positive or negative contribution
to the total torque generated by the rotor the winglet induces a radial force
and a twisting moment about the axis of the main blade. These forces are not
threated in this work but are important for the aeroelastic behavior of the blade.
Especially the twisting moment can be expected to be relatively large when
compared to the existing twisting moments found in present turbine blades.
1.3 Overview of the Methods Applied
Because the winglet works by altering the flow and pressure properties it can
be modeled using inviscous theory. It is important that the blade is not stalled,
but in practice windturbines are often stalled near the root. Inviscid theory can
therefore only to limited extend be applied to accurately predict the performance
4 Introduction
of windturbines. The influence of winglets are confined to the outer parts of the
blade where the assumption of inviscous flow is valid and inviscous theory can
be used to study the local variations. The tip speed on turbines is limited to
70 m/s which is approximately 21% of the speed of sound which is well within
the range of incompressible flow (¡ 30%). Locally, near the profile, the velocities
can be much larger and the flow can not be considered incompressible. On the
profile suction peak the local velocity can be as high as twice the incoming flow
velocity. Despite of this the flow will be considered incompressible everywhere
and potential theory will be used.
BEM and actuator disc methods can not predict variations in flow properties
due to winglets and this leaves computational fluid dynamics (CFD) as the only
established alternative. It is desirable to have a fast method for parameter
studies and design validation and this is not the case with CFD. The available
potential methods are roughly; 1) lifing line methods 2) Vortex Lattice methods
3) Panel methods. All 3 methods has been studied and all has been found to have
weaknesses. The Vortex Lattice method is relatively slow and the induced drag
forces are not easily evaluated. The panel method can accurately calculate the
forces but long computational times makes it unusable. Both methods requires
a wake model, and a Free Wake model has been implemented. It is difficult and
time consuming to model the wake accurately and this adds to the uncertainty of
the potential flow methods. The Vortex Lattice method was eventually chosen
over the panel code. The lifting line method is simple and fast, and the accuracy
is relatively good. It is different from the Vortex Lattice method in that the
chord and twist distribution is not part of the solution. Instead the velocity and
bound circulation distributions are found and an equivalent full 3D geometry
can then be defined, something which is relatively simple as long as there is
no substantial spanwise flow. From the velocity and circulation distribution it
is also possible to evaluate all important forces. It is therefore clear that the
velocity and circulation distribution describes the physics along the blade and
the goal of the design process is to optimize these for a given geometry. Note that
in viscous flow other non linear properties are important as well. Because the
velocity and circulation distributions are so important the Vortex Lattice results
are also converted to lifting line data, i.e. chordwise variations are neglected and
the data is distributed along a line describing the span. This makes it easier
to study different designs since data distributed over a Vortex Lattice surface
are hard to interpret. It is also worth noticing that the forces evaluated along a
lifting line includes the leading edge suction, something which requires special
treatment in Vortex Lattice codes. To sum up, the final method consists of a
Vortex Lattice computation of distributed velocities and circulation using a Free
Wake model. Followed by transfer of data to a lifting line representation where
forces are evaluated and results can compared to design input.
Despite a large effort the Vortex Lattice results were never really satisfactory,
1.4 Outline of the Thesis 5
something which is mainly due to the uncertainties of the wake model. Because
of this and based on the experience gained a simpler model was made which
is based on a known velocity and circulation distribution over a traditional flat
blade. The task is then to extend the blade and circulation distribution in order
to simulate a winglet and then predict the changes in velocities due to this.
Under certain conditions it is only necessary to recalculate the influence of the
wake immediately behind the blade and one quarter revolution downstream.
This makes the algorithm fast and, because the wake is prescribed, the uncer-
tainties of the free wake model is avoided. The circulation distribution is based
on results for optimal designed wings, which is allowed as long as the changes
in velocities along the blade is predicted correct. The method is used to design
the winglets and predict the performance. Vortex Lattice simulations are used
for validation.
The methods has been implemented in Matlab and Fortran. Practical consid-
erations and source code has been completely left out of this report.
1.4 Outline of the Thesis
The contents of this thesis are presented in six parts.
Part 2 presents the applied potential flow theory. Coordinate systems and di-
mensionless variables are also introduced here.
Part 3 describes the lifting line method. The method is applied to wings in
steady forward flight for which the optimum design of winglets can be found.
The results are generalized to wind turbines.
Part 4 describes the free wake vortex lattice and panel method. The formulation
applies to panel methods for which the Vortex Lattice method is a special case.
Part 5 describes the design process. Some important general issues regarding
the choice of design parameters are discussed.
Part 6 discusses the predicted performance of 2 selected final winglet designs.
Part 7 contains conclusions.

6 Introduction
Part 2
Theory
This part introduces coordinate systems, dimensionless variables and the applied
theory for later reference.
2.1 Definition of Coordinate Systems
Two coordinate systems are defined as seen in Figure 2.1.

Z
z
y
θ Y
Ω
X, x
Figure 2.1: The initial and the rotating coordinate system.
The initial coordinate system corresponds to an earth reference frame. The

8 Theory
coordinates are denoted with capital letters.
R = {X, Y, Z}T
The body coordinate system is rotating about the X-axis of the initial system
with the same rotational speed as the turbine rotor which is stationary in this
frame of reference. The coordinates are denoted with small letters.
r = {x, y, z}T
The rotational speed is Ω and the motion is therefore completely described by

the rotation vector Ω defined in the initial coordinate system.
Ω = {Ω, 0, 0}T
2.2 Definition of Blade and Winglet Geometry
In the body fixed coordinate system the blade is positioned as shown in Figure
2.2. The radius is R. The winglet is extending downstream in the positive
x-direction. The contour of the winglet is described by its height, h, and the
curve radius of the bend ,Rc .
z
Q∞ y
y
h
Rc
x
x
Figure 2.2: Definition of geometry. The winglet contour is shown to the right.
The twist, twistz , and the chord, c, is described along the blade. Notice that a
special definition of the twist is used, which is seen in figure 2.3.
2.3 Dimensionless Variables 9
twistz
y
x twistz
Figure 2.3: Definition of twist angle, twistz .
2.3 Dimensionless Variables
Results will be presented in dimensionless form. Under the assumptions of linear

theory and flow similarity the rules for scaling are given below. Coordinates are
scaled by the turbine radius:
y
Spanwise coordinate
R
x
Downstream coordinate
R
s
Surface parameter along span
R
Velocities are scaled by the wind speed:
u
x-velocity component
Q∞
v
y-velocity component
Q∞
w
z-velocity component
Q∞
Tip speed ratio, thrust and power coefficients. Note that section Cp and CT ,
usefull for describing the distribution of loads, will be defined in section 3.7.
ΩR
λ=
Q∞
Mx Ω
Cp =
1/2ρQ3∞ A
Fx
CT =
1/2ρQ2∞ A
10 Theory
The bound circulation can be scaled by the following relation which holds for
wind turbines with similar flow. I.e. equal Cp , λ and geometrical contour.
Γ
= Constant
Q∞ R
2.4 Governing Equations
The following theory applies to an inviscid, irrotational and incompressible fluid,

where a velocity potential, Φ can be defined
∇Φ = Q
The potential obeys Laplaces equation in both stationary and time dependent
flows.
∇2 Φ = 0 (2.1)
2.5 The General Solution and Boundary Condi-

tions
Equation (2.1) is a linear partial differential equations. Its general solution and
boundary conditions are described below and forms the basis of all the applied
mathematical techniques.
Boundary Conditions
For a body immersed in a fluid the boundaries are the body surface and infinity,
and in general both Dirichlet and Neumann boundary conditions can be defined.
The first condition is used in potential based panel methods where the potential
is specified directly on the boundaries and the potential is solved for. The
second is used in methods where a condition on the flow itself is specified on the
boundary. The latter method is used in this work. The first boundary condition
is tangential flow at the surface which can also be stated as zero normal flow
BC1: (Q − v) · n|surf ace = (∇Φ − v) · n|surf ace = 0 (2.2)

2.5 The General Solution and Boundary Conditions 11
n is an outward pointing normal vector and v is the structural velocity. In the

case of a rotating body this is
v =Ω×R
The second boundary condition is that the disturbance of the body must vanish
at infinity
lim lim
BC2: R→∞ Q =R→∞ ∇Φ = Q∞
Notice that in lifting line methods the boundary conditions may be defined
differently. For example by defining the downwash.
The general solution
The velocity potential can be written as the sum of a constant farfield potential
(Φ∞ ) plus a perturbation potential (Φp )
Φ = Φ∞ + Φp
Φ∞ = Q∞ x
The general solution of the perturbation potential is an integral of the potentials

from doublets and sources, the strength of which are µ and σ respectively,
distributed over the body surface and the wake
Z Z
1 1 1 1
Φ(X, Y, Z) = µn · ∇ dS − σ∇ dS (2.3)
4π body+wake r 4π body r
Notice that the wake is assumed to consist of doublets only. For a derivation of
this equation see e.g. Katz and Plotkin[11]. Doublets and sources are singular
elements inducing potential fields which individually are solutions of (2.1). They
are described in section 2.7.
2.5.1 The general solution on discrete form
Using the knowledge that the perturbation potential can be constructed by

distributing sources and doublets it is also possible to construct the velocity at
a given position as the sum of perturbation velocities induced by doublets and
sources (Qµ and Qσ ) plus Q∞
X X X
Q(X, Y, Z) = Qµ + Qµ + Qσ + Q∞
body wake body
12 Theory
The velocity induced by a singularity element is proportional to the strength of

the element. This means that the equation can be written on the form
N
X NW
X N
X
Q(X, Y, Z) = Qk µk + Ql µl + Qk σk + Q∞
k=1 l=1 k=1
Where indices k and l represents the singularity elements of the body and wake.
N and NW are the number of structure and wake elements respectively. Qk
and Ql is the induced velocity at point (X,Y ,Z) due to unit strength singularity
elements. Inserting this into the Neumann boundary condition (2.2) yields
N
X NW
X N
X
Qk · n µk + Ql · n µl + Qk · n σk = v · n − Q∞ · n (2.4)
k=1 l=1 k=1
This equation has to hold on the body boundary. The distribution of wake
singularities is also important.
2.5.1.1 Remarks
Equation (2.4) forms the basis of the numerical panel method. In the form
given v can include any motion of the body relative to the initial coordinate
system, but in the following only rotation will be included. The implementation
of the equation is discussed in detail in part 4. The vortex lattice method is
a special case of equation (2.4). In this method the sources are omitted and a
flat geometry is chosen. This means that once equation 2.4 is implemented in a
panel code the vortex lattice problem can also be solved.
2.6 The Kutta Condition
The general solution (2.3) is not unique and a condition has to be imposed.
This is the Kutta condition which requires that the flow leaves the trailing edge
smoothly. This is physically reasonable for most flows when the relative velocity
in the flapwise direction is small compared to the velocity in the chordwise
direction. It can be shown that the Kutta condition can be implemented by
setting the vorticity at the trailing edge equal to zero
µT.E. = 0
2.7 Singularity Elements 13
2.7 Singularity Elements
The singularity elements used in this work are described in this section. The
potential fields induced by the elements are all individual solutions of (2.1). For
a more thorough presentation see e.g. Katz and Plotkin[11].
2.7.1 Constant Strength Vortex
A straight vortex with constant strength Γ is the basis of many numerical so-
lutions. The velocity induced at r0 (see Figure 2.4) is given by the Biot-Savart
law
Z
Γ dl × (r0 − r)
q= 3 (2.5)
4π |r0 − r|
Γ dl × (r0 − r)
⇔ dq = (2.6)
4π |r0 − r|3
Where r describes the vortex segment of length dl. In numerical computations
the start and end points of a vortex is given (r1 , r2 ), and a formula for the
velocity induced at a given point can be found in Katz and Plotkin[11].

Γ r1 × r2 r1 r2
q= r0 −
4π |r1 × r2 |2 r1 r2
Alternatively the induced velocity can be found from the view angles β1 and
β2. They are the angles between the vortex segment and the evaluation point
evaluated at r1 and r2 .
Γ
q= (cos β1 − cos β2 ) (2.7)
4πd
Where d is the perpendicular distance from the evaluation point to the vortex
segment.
z β2
r0 r2
Γ
β1
r1
x y
Figure 2.4: A straight vortex of constant strength Γ.

14 Theory
2.7.2 Constant Strength Doublet
The doublet element is equivalent to a distribution of vortices. If the strength

is constant the vortices cancels everywhere except at the edges, and the element
is equivalent to a vortex ring. The induced velocity can be therefore be found
by combining straight vortices of the type descriped in section 2.7.1. Figure 2.5
illustrates this. Even though the element is depicted as flat it can be twisted in
space and it can have any number of edges. Notice that the doublet panel has
an associated vector. The induced flow is being driven through the panel in the
direction of the vector and back around it.
z
µ
4 3 4 3
y µ
µ µ
=
x
1 1 µ
2 2
Figure 2.5: A constant strength doublet element and its vortex ring equivalent.
2.7.3 Constant Strength Source Panel
The potential induced at (x, y, z) by a constant strength source panel of strength

σ positioned in the x-y plane is
Z
−σ dS
Φ(x, y, z) = p (2.8)
4π S (x − x0 )2 + (y − y0 )2 + z 2
The solution in closed form for a 3 node panel (Figure 2.6) is given in appendix
A.1. The flow induced by the panel is everywhere directed away from it.
y
1 3
σ
Figure 2.6: A 3 node sourcepanel.

2.8 Kutta-Joukowsky and Helmholz Theorems 15
2.8 Kutta-Joukowsky and Helmholz Theorems
The Kutta-Joukowsky theorem states that the force on a constant strength di-
rected vortex element subjected to a relative velocity of Q is
F = ρQ × Γ (2.9)
The Helmholz theorem states that the strength of a vortex line is constant along
its length, and it must form a closed loop or extend to infinity.
2.9 Wake Strength
The wake consists of doublet panels which are equivalent to trailing and shed
vortex lines. The trailing lines are parallel to the flow and the shed lines are
parallel to the trailing edge before they are convected downstream.
In the panel and vortex lattice methods the Kutta condition is imposed and
this indirectly determines the strength of the trailing vortices. In order to set
the vorticity at the trailing edge equal to zero it is necessary to attach doublet
panels which cancels the vorticity and this dictates the strengths of the panels.
In the following, since constant strength doublet panels are used and there is no
time dependency the shed vortices corresponding to the panel edges will can-
cel themselves everywhere, and the wake effectively consists of vortices trailing
downstream. This is the same as attaching a horseshoe vortex to the trailing
edge, with the same strength as the trailing edge vorticity.
In the lifting line method the shed vortices cancels by the same argumentation as
above. The strength of the trailing vortices are determined using the Helmholz
theorem. Since the circulation can not end abruptly any change in circulation
along the lifting line must be convected downstream in the wake. I.e.
dΓ
ΓW (s) = (2.10)
ds
Notice that the sign will depend on the defined positive direction of the trailing
wake vortices.
16 Theory
2.10 Wake Shape
The wake is considered to be force free. The Kutta-Joukowski theorem states

that if the force on a wake vortex line is zero then
Q × ΓW = 0
I.e. that the velocity is parallel to the circulation vector everywhere in the wake
Q||ΓW
2.11 The Unsteady Bernoulli Equation
If the flow is time dependent the unsteady Bernoulli Equation 2.11 must be used
for pressure calculations in a panel code.

p q2 ∂Φ p q2 ∂Φ
gz + + + = gz + + + (2.11)
ρ 2 ∂t ρ 2 ∂t ∞
In the following the effect of gravity will be neglected. The equation will be used
in the inertial (X,Y,Z) system where the flow far away is homogenious so that
Φ|∞ =const and the time dependent term is zero. The equation is rewritten as
p − p∞ 2 ∂Φ
1 = q∞ − q2 − 2
2ρ
∂t
The kinematic velocity is defined as the undisturbed velocity relative to the

blade
Qkin = Q∞ − Ω × r (2.12)
The pressure coefficient is defined using the kinematic velocity as reference. I.e.
qref (r) = ||Q∞ − Ω × r||
And the pressure coefficient is

p − p∞ 1 2 ∂Φ
1 2 = 2 q∞ − q2 − 2 = Cp
∂t
2 ρ ||Q∞ − Ω × r|| ||Q∞ − Ω × r||
This definition makes the pressure coefficients comparable over the span of the
wing even though the pressures may vary significantly. If Ω=0 the expression
reduces to the standard definition of the pressure coefficient.
2.11 The Unsteady Bernoulli Equation 17
The term ∂Φ∂t can be evaluated in different ways. The method described here is
based on a turbine rotating steadily in a plane perpendicular to the onset flow
which is homogenious. I.e. no wind shear, coning angle, tilt angle and no change
in time of pitch, rotational speed, windspeed etc. Except that no wind shear
is allowed this approximates the working conditions of a wind turbine. Under
these conditions the flowfield seen from the rotating body coordinate system is
steady in time and therefore is
∂Φ
=0
∂t body
The relation between the time derivative in the inertial system and in the body
coordinate systems is derived in [11] page 372

∂ ∂ ∂ ∂ ∂
= −Ω × r · , , +
∂t inertial ∂x ∂y ∂z ∂t body

∂Φ ∂Φ ∂Φ ∂Φ ∂Φ
⇒ = −Ω × r · , , +
∂t inertial ∂x ∂y ∂z ∂t body
The space derivatives are in the body-fixed coordinate system and are therefore
the velocity components in the instantaneous directions of the coordinate vectors
of the rotating system
∂Φ
= −Ω × r · (u, v, w)|body
∂t inertial
The (u, v, w)|body term is not the velocity relative to the rotating system, which
is different from the velocity in the inertial system, but is simply the velocity
evaluated in the inertial system divided into components in the instantaneous
directions of the rotating system.
18 Theory
Part 3
Numerical Lifting Line Model
The lifting line model is a classical model that has been used for calculating
properties of airplanes, propellers and also, to limited extent, windturbines.
The basic idea is to represent all properties on a line along the span, meaning
that only the overall shape of the wing is represented geometrically. The force
system on the blade is described by the velocities relative to the blade and the
bound circulation.
Basic properties of winglets are deduced but the principles of the method will
also be used in part 4 and part 5.
3.1 Description
A numerical lifting line method capable of solving flows over wings of arbitrary
shape is presented here. The wing is required to have a large aspect ratio so
that the chordwise distribution of singularities on the wing can be neglected.
The aspect ratio, AR, is defined as the ratio of the span to the planform area.
The mathematical solution is build from horseshoe vortex elements extending to

infinity. The horseshoe vortex consists of three straight vortices with constant
20 Numerical Lifting Line Model
and equal strength. Two of the vortices extends to infinity and the horseshoe
vortex is therefore in accordance with the Helmholz theorem. In numerical
applications it is enough to let the tail extend a distance far downstream. Figure
3.2.a shows the configuration. The two tails on the horseshoe vortices represents
the wake and the centerpiece represents the bound circulation on the wing.
Because the wake is parallel to the incoming flow it is force free. The bound
circulation is perpendicular to the incoming flow and is therefore subject to a
lift force. The downwash is induced by the wake and is perpendicular to both
the wing and the incoming flow. This component is responsible for the induced
drag.
The wing lifting line is discretized into collocation points and vortex points. The
collocation points are the points of evaluation of velocities. The vortex points
are the start and end points of the bound circulation vortex elements, and the
positions from which the wake vortex lines are attached. The collocation points
and vortex points are indicated with black and red dots respectfully, as seen in
figure 3.1.
Figure 3.1: A wing and wake system represented by horseshoe vortices.
The collocation and vortex points are distributed according to a cosine spacing
scheme. It will be shown that this ensures a fast convergence of the solution.
The formula used to discretize a coordinate, s, in the range s ∈ [a, b] is
cos(xt1 ) − cos(xi )
si = a + · (b − a) (3.1)
cos(xt1 ) − cos(π − xt2 )
Where x is a linearly varying variable in the range x ∈ [0 + xt1 , π − xt2 ]. xt1 and
xt2 are truncation values that can be used to alter the shape of the distribution.
They must be chosen so that x are increasing.
Flow properties are evaluated in the nc collocation points. Since there are nc
unknowns, i.e. the strengths of the horseshoe vortices, nc equations must be
defined in the collocation points. The imposed condition is Munks theorem 2
(see Cone[2]). It states that the induced drag will be a minimum when the
component of induced velocity in the direction normal to the wing (i.e. the
3.1 Description 21
normalwash wN ), is proportional to the cosine of the angle of inclination, γ

w0
wN = Qi · n = cos(γ) (3.2)
2
Where Qi is the wake induced velocity and n is a normal vector. Figure 3.2
illustrates this. The blue arrows indicates the induced velocity components
normal to the wing plane according to Munks theorem 2. The magnitude of
w0 is proportional to the wing loading and is not important in this context
because only geometrical properties, which are independent of the loading, will
be studied.
Figure 3.2: The defined optimum normal-downwash distribution
Since the induced velocity of a vortex line is proportional to the strength off it
the total wake induced velocity in a given point can be written as
NW
X
QW = Ql Γl
l
Where Ql is the velocity induced by wake vortex l with unit strength. Com-
bining this with (3.2) yields an equation that must hold in every collocation
point.
NW
X w0
Ql · n Γl = cos(γ(s))
2
l
Writing this equation in all collocation points forms a linear system of equations
with the strengths as the unknowns. Ql can be evaluated using the formulas
given in section 2.7.1.
The method will be validated against analytical results presented in section 3.2.
Remarks
The forces evaluated using the described method includes the force perpendicu-
lar to the local wing plane due to the incoming flow, and the induced drag due to
velocity components induced by the wake. There is a third velocity component,
the self induced velocity, which is induced by the lifting line onto itself. This
force can not simply be evaluated by the Kutta Joukowsky theorem, and will
be omitted in the following.
3.2 Classic Analytical Lifting Line Results
The simple geometry of airplane wings allows for an analytical treatment if the
wake can be represented in a simple form as well. The plane, straight wing with
a wake trailing downstream is a classical case for which analytical results exists.
An elliptic lift distribution is obtained if the planform itself is elliptic. I.e. if

the chord is elliptic
s 2 s 2
y y
c(y) = C0 1− ⇔ Γ(y) = Γ0 1−
b/2 b/2
Γ0 is the circulation at the center of the wing. The downwash over an elliptic
wing is constant and equal to
Γ0
wi = − (3.3)
2bQ∞
It can be shown that a wing with an elliptic lift distribution has the lowest
induced drag of all flat wings. This is expressed in the following formula which
also defines the efficiency factor k
CL2
CD = , elliptic wing: k=1 (3.4)
πARk
For all plane wings k is less or equal to unity. This formula holds for non planar
wings as well but k can then take values large than unity. The equation is
derived in Cone[2]. Note that the aspect ratio is based on the projected area
and span of the wing. k is a dimensionless parameter describing the effect of
the geometry on the induced drag. It does not include any effect of viscous drag
forces.
3.3 Validation 23
3.3 Validation
A flat wing of span b will be used as a test case to validate the numerical lifting
line code. Q∞ =1 m/s. Since the wing is flat γ=0 everywhere and the optimal
normal downwash according to (3.2) is constant and in the negative z-direction.
This corresponds to an elliptic wing with an elliptic circulation distribution.
The magnitude of the downwash is selected as W0 /2=1 m/s. According to (3.3)
the maximum circulation is max(Γ0 )=4 m2 /s. Figure 3.3.b shows the calculated
elliptical circulation distribution. The relative error of the maximum value is
-0.0168% of the exact value of Γ0 =4 m2 /s. Figure 3.3.a shows the calculated
efficiency factor, k, against an increasing number of collocation points, nc . The
exact analytical result is k=1 and the relative error of the numerical result is
-0.0252% for nc = 80.
nc=80, max(ΓB)=3.99933 m2/s

1
5
0.999 4
0.998 3
ΓB [m2/s]
k [−]
2
0.997
1
0.996
0
0.995 −1
20 30 40 50 60 70 80 −1 0 1
nc [−] 2y/b [−]
(a) Convergence with increasing spanwise (b) The circulation distribution with
discretization. nc =80.
Figure 3.3: Validation results.
It is concluded that the numerical lifting line is accurate and that a relatively
coarse discretization is sufficient. In the test example above a cosine spacing was
used with success. Calculations has been made on a wing with rounded winglets
attached. The winglet height is h=0.2 b and the curve radius is R=h. Figure
3.4.a compares the computed k-values using the linear and cosine discretization
respectively. The values converges very fast to k=1.2527 when the cosine spacing
is used. Further tests shows that 2000 collocation points are needed for the
results of the linear spacing to converge to k=1.2530.
1.3
Linear spacing
Cosine spacing
1.29
1.28
k [−]
1.27
1.26
1.25
20 40 60 80 100 120 140 160 180 200
n [−]
c
Figure 3.4: Convergence of k for a wing with rounded tips using cosine and
linear spacing.
3.4 Effect of Winglet Height and Curve Radius
A parameter study is made for a wing with vertical winglets attached to the
suction side. The corner is rounded with varying curve radius. This shape is
believed to represent important geometrical properties. The following quantities
has been varied in the study
R
h The ratio of curve radius to winglet height. A value of 1 is equivalent to
a completely rounded tip. Rh ∈ [0.25, 0.95].
b
h The ratio of wingspan to winglet height. Large values means small winglets.
b
h ∈ [5, 75].
Figure 3.5 shows the results. The value of k is highest for relatively large
winglets. The effect of curve radius is not as significant but is important for
small values of hb . In general, the smaller curve radius the better, but in a
viscous flow there is an increased risk of separation in a corner and the best
compromise is to round the corner.
3.5 Optimum Circulation distribution on Wings 25
1.4
k [−] 1.3
1.2
1.1
1
0.8 20
0.6 40
0.4 60
b/h [−]
R/h [−]
Figure 3.5: Results of parameter study.
3.5 Optimum Circulation distribution on Wings
Calculations has been made on a wing with 5% winglets (h/b=0.05) and varying
curve radius. The bound and wake circulation strength will be studied. The
bound circulation is the solution from the numerical lifting line code. The wake
strength is found as
dΓb
ΓW = −
ds
Where s is a surface parameter along the span. Note that the wake strength
found in this way is not the strength of the discretized horseshoe vortices but is
instead the distributed strength or strength per length.
Figures 3.6, 3.7 and 3.8 shows the results and the contours of the wing. The
wake strength, ΓW , increases locally in the bend, and especially for small curve
radius the increase is large. This is caused by the bound circulation distribution,
ΓB , having a larger derivative locally in the bend. In all cases the wake strength
has a maximum at the winglet tip and this indicates that the bound circulation
distribution should also have a large derivative near the winglet tip. Except for
local changes in the bend it appears that the optimum bound circulation can be
found by simply extending the elliptic circulation distribution over the winglet.
Winglets with a sharp bend is clearly more effective than rounded winglets. In
this case k=1.056 for Rc /h=0% and k=1.043 for Rc /h=80%.
2
k=1.0562 max(Γb)=4.115m /s
0.8
Γb/max(Γb) [−]
0.6 −Γ /max(Γ )/10 [−]

W b
0.4
0.2
0.8 0.9 1
2s/b [−]
Figure 3.6: Circulation distribution over wing with 5% winglet and 0% curve
radius.
2
k=1.0516 max(Γb)=4.106m /s
0.8
Γb/max(Γb) [−]
0.6 −Γ /max(Γ )/10 [−]

W b
0.4
0.2
0.8 0.9 1
2s/b [−]
radius.
2
k=1.043 max(Γb)=4.089m /s
0.8
Γb/max(Γb) [−]
0.6 −Γ /max(Γ )/10 [−]

W b
0.4
0.2
0.8 0.9 1
2s/b [−]
radius.
3.6 Generalization of Wing Results to Turbines 27
3.6 Generalization of Wing Results to Turbines
Some of the results found in this part can be generalized to turbines. Based on
the discussion and results of part 5 the winglet acts locally and has no effect
on the wake except for immediately behind the trailing edge. It is therefore
expected that physical properties of turbine winglets resembles those of wing
mounted winglets.
The circulation distribution over well designed turbine blades differs substan-
tially from the elliptic distribution over wings. Turbine blades tends to have
a constant circulation over the center part but near the tip the shape approxi-
mates that of wings. Therefore it can be expected that the circulation on well
designed winglets should be distributed somewhat similar to what was found
in section 3.5. This is important from a design point of view and will be used
in part 5 to define the circulation on the outer part of the blade. Note that a
theorem similar to Munks theorem 2, to this authors knowledge, does not exist
for turbines and an optimum distribution of circulation can not be calculated
directly.
Winglets are found to be most effective if the bend is sharp with a small rounding
radius. This is believed to be true for turbines as well.
3.7 Force and Viscous Drag Calculation Using

Lifting Line Data
In the previous sections only the geometrical properties relatated to the induced
drag was investigated. The viscous drag is several times larger than the induced
drag and will be included in the calculations henceforth. A method for cal-
culation of the viscous force on a lifting line is described below. The inviscous
force on a bound vortex of finite length is calculated using the Kutta Joukowsky
theorem (2.9). To get the section force, f (s), divide the force by the length of
the vortex 1
ρQ × Γ
f (s) =
||Γ| |
The total force and moment is found by integration over the lifting line
Z
Fj = f (s)ds (3.5)
1 The bound vorticity is assumed to be discretized into constant strength vectors, Γ

Z
Mj = r(s) × f (s)ds (3.6)
The 3D numerical, inviscid methods are capable of calculating induced drag,

and this component is directly included in the calculated force. The missing
drag components are the pressure, or form drag, and the friction drag. Both of
them can be estimated using 2D experimental data. Both drag forms are due
to viscous effects, and the total drag is termed viscous drag. Since the effect of
the 3 dimensional wings are effectively to change the local flow properties, care
has to be taken so that the correct reference velocities are used when accessing
2D data.
A method based on lifting line data and Lift-Drag polars are presented in the
following.
The viscous drag is by definition in the direction of the onset flow. In experi-
ments the onset flow is often uniform, but here the local flow velocity Q must
be used. The direction is
Q(s)
t(s) =
||Q(s)||
Q must include all velocities. Since the theory is not linear it is in general not
possible to decompose it, the exception being when a constant lift to drag ratio
is assumed. The velocities induced by the wing itself represents a numerical
problem and can not be evaluated in the lifting line method. The direction
of these velocities are in most cases in the direction of t, and does not affect
the solution (i.e. the Γ-values), but they will have an effect on the reference
velocities and should therefore be estimated and included in the viscous drag
model if possible. As an approximation, if the winglets are small, they can be
neglected and will be in this work. The section force, f (s), is used to calculate
the section lift coefficient, Cl . Notice that this local force is rotated by the local
wing angle, γ, and the induced angle of attack. f (s) is assumed to be known
along the lifting line. The reference velocity in the definition of Cl is the local
flow velocity.
qref = ||Q(s)||
The section lift coefficient is

f (s)
Cl (s) = 2 c(s)
1/2ρqref
Using 2D table look up the section drag coefficient is found as a function of Cl
Cd (s) = Cd (Cl (s))

3.7 Force and Viscous Drag Calculation Using Lifting Line Data 29
The viscous section drag, d(s), is then determined from its definition
d(s)
Cd (s) = 2 c(s)
1/2ρqref
Finally the total friction force and moment can be integrated

Z
Fv = d(s)t(s)ds (3.7)
Z
Mv = r(s) × [d(s)t(s)] ds (3.8)
For turbines the total power and thrust coefficients are based on the moment
about the x-axis and the x-component of the forces
(Mj,x + Mv,x ) Ω
Cp = (3.9)
1/2ρQ3∞ A
Fj,x + Fv,x
CT = (3.10)
1/2ρQ2∞ A
Where A is the swept area. A blade section contributes with the moment, Mx′ ,
about the x-axis and the section power coefficient, Cp′ , is defined as
Mx′ Ω
Cp′ = (3.11)
1/2ρQ3∞ A
Likewise the section thrust coefficient is defined as

Fx′ Ω
CT′ = (3.12)
1/2ρQ2∞ A
The section power and thrust coefficient are dimensionless torque and force
values which will be used for presentation of results.
Validation
The method is validated against work of Gaunaa et al.[5]. They present the
Cp -value at a tip speed ratio of λ=8 calculated using a free wake lifting line
method (FWLL) and CFD on the corresponding 3 dimensional geometry. The
induced velocities and the circulation distribution along the blade are used as
input and can be found in section 4.5.
A constant lift to drag ratio of 110 is used (i.e. Cd = Cl /110). In table 3.1
the reported results are compared to equation (3.9) and (3.10). The Cp and
CT value is 0.4% and 0.7% less than the reported CFD result. Based on this
relative good agreement the method is considered validated. The error on the
thrust is to large to be acceptable but since the major concern in the following
is the power production it will be accepted. There is a possibility that the
disagreement is because of the differences in flow properties predicted by the
various models.
Cp CT
CFD, Gaunaa et al.[5] 0.515 0.872
FWLL, Gaunaa et al.[5] 0.514 0.868
(3.9) using FWLL data, (3.10) 0.513 0.866
Table 3.1: Validation of (3.9) and (3.10).
3.8 Conclusions
In this part preliminary results for winglets are found. The lifting line represen-
tation of data will be used in the following parts for calculation of forces and as
the basis of the design method. The bound circulation distribution over optimal
loaded wings was found and will be used in the design algorithm.
Part 4
A Freewake, Vortex Lattice
and Panel Method
The lifting line method presented in part 3 can be used very successfully for high
aspect ratio wings but lacks the ability to accurately model wings of low aspect
ratio. Even though modern turbine blades are very slender, they do have signifi-
cant chord lengths close to the hub, where the lifting line representation can not
be considered accurate. The vortex lattice method and the panel method both
includes the effect of the chordwise distribution of properties and the methods
are presented in the following.
It turns out that both the vortex lattice and the panel method has severe short-
comings when simulating the forces on rotating turbines. Mainly, the degree of
accuracy when evaluating forces, and the fact that a turbine is almost certain to
stall over a smaller or larger part of the blade near the hub, making the inviscid
flow assumption invalid, has a pronounced effect on the accuracy. The accu-
racy can be improved by using higher order methods, but this is hardly worth
the effort since there is still a number of issues regarding the validity of using
potential flow for simulation of windturbines. Today it is possible to make a
full CFD calculation in 24 hours and this makes it a more appealing alternative
to use simple fast models for initial design and CFD calculations for validation
and final designs.
The method presented is based on constant strength doublet and source panels
32 A Freewake, Vortex Lattice and Panel Method
and is termed a first order method. Higher order methods uses varying strength
panels. The implemented code can be used as both a panel and a vortex lattice
method.
4.0.1 Description
Equation (2.4) repeated below states the condition of zero normal flow on the
body surface and forms the basis of a numerical solution
N
X NW
X N
X
Qk · n µk + Ql · n µl + Qk · n σk = v · n − Q∞ · n (4.1)
k=1 l=1 k=1
To implement it it is nescessery to define the blade and wake singularities. First

consider the N blade source and doublet elements. The blade is described by
nodes that are formed into panels consisting of 4 nodes each. The panels can be
twisted in space which means that any geometry can be represented, as opposed
to using flat panels only. A collocation point is placed at the center of each panel
and it has an associated normal vector pointing away from the structure. The
constant strength doublet panel (section 2.7) is defined based on the 4 corner
points of a panel. This is illustrated in Figure 4.1 along with the collocation
point and normal vector.
3
2 C 4
1
Figure 4.1: A constant strength doublet panel.
A 3 node source-panel has been implemented in order to be able to describe

the twisted source panel and the panel is split into 4 panels with a common
center node. The center node is placed a small distance below the collocation
point. This is important since if the collocation point had been below the panel
it would appear in the equations that the panel is creating suction instead of
blowing. The constant strength doublet panel is equivalent to a vortex ring,
and therefore is the exact placement of the collocation point relative to it not
important, as long as it is approximately at the center.
33
N collocation points is now defined on the body surface which consists of N

panels. Since each panel has associated both a µ and a σ value the number
of unknowns are twice the number of available linear equations. Therefore,
the source strengths can be defined and only the doublet strengths are solved
for. In the Neumann formulation of the problem used here, the strengths can
be choosen freely. However they should be choosen wisely in order to avoid
numerical problems. The source strengths are set equal to the negative of the
magnitude of the normal flow due to freestream and structural velocities
σ = −n · (Q∞ − v)
Directly above a panel the velocity induced by it, is
1
Qi = σn
2
1
⇒ Qi = − (n · (Q∞ − v)) · n
2
Therefore, Q is exactly half the normal flow due to (Q∞ −v) and in the opposite
direction. This choice is made because it makes the solution numerically stable.
Evaluating the terms related to body singularities in the N collocation yields a
linear system of equation which can be put on matrix form
 
:
PN PN  = Aµ µ + Aσ σ (4.2)
k=1 Qk · n µk + k=1 Qk · n σk
:
Where µ and σ are vectors of the panel doublet and source strengths.
The wake is modeled as a distribution of constant strength doublet panels. Since

only steady state solutions will be considered (i.e. no change in time of the bound
circulation) the doublet panels trailing of the same spanwise position will have
the same strength and the circulation at the edges they have in common cancels
everywhere. The wake is therefore equivalent to horseshoe vortices attached
to the trailing edge. The strengths of the horseshoe vortices are defined by
the Kutta condition which states that the circulation at the trailing edge must
cancel. In order to do so the wake horseshoe vortex must have a strength that
depends on the trailing edge doublet strengths
µl = µk1 − µk2
This is illustrated in Figure 4.2 which shows two adjoining panels at the trailing
edge and an attached horseshoe wake element extending downstream.
µk1
µk2 TE
µl =µk1 -µk2
Figure 4.2: A horseshoe vortex attached to the trailing edge.
In (4.1) the influence of the wake is included in the term

NW
X
Ql · nµl
l=1
The contribution from a horseshoe vortex to the sum is

NW
X
Ql · nµl (µk ) = · · · + Ql · n(µk1 − µk2 ) + . . .
l=1
And therefore will the wake influence effectively depend on the strength of the
blade doublet panels at the trailing edge. If a flat vortex lattice geometry
is used the wake strength will only depend on the strength of one panel. In
practice the wake consists of a large number of vortices and only the ones directly
attached to the trailing edge depends on the trailing edge doublet strengths.
The downstream vortices have an absolute fixed strength. This change in wake
strength downstream is caused by the initial developing flow of a time marching
scheme. In the developing phase the horseshoe representation is not valid, but
this part of the wake is quickly convected away from the blade and is eventually
cut away. The total influence of the wake is split into the NT E dependent
vortices at the trailing edge and the remaining constant strength vortices far
downstream
NW
X N
XTE NW
X
Ql · nµl = Ql · nµl (µk ) + Ql · nµl
l=1 l=1 l=1+NT E
Evaluating this term in the N collocation points yields a linear system

 
:
PNT E PNW
l=1 Ql · nµl (µk ) + l=1+NT E Ql · nµl = AW µ + AW 2 µW (4.3)
 
:
Where µW is a vector with the strengths of the downstream independent wake

vortices. The shape of the wake has a significant impact on the solution and
4.1 Freewake Model 35
care must be taken in order to determine it. This is especially important when
calculating the flow past wind turbines, but for wings in steady forward fligth
good results can be obtained simply by letting the wake trail downstream. The
calculation of the wake shape behind wind turbines is discussed in section 4.1.
Combining (4.1), (4.2) and (4.3) yields

(AW + Aµ )µ + AW 2 µW + Aσ σ = a
Where a is the vector of v · n-Q∞ · n evaluated in the collocation points. This
is rearranged as
(AW + Aµ )µ = a − AW 2 µW − Aσ σ
Aµ = b
This system can now be solved for the unknown doublet distribution.
Notice that Aµ and Aσ are constant when there is no deformation of the struc-
ture (rotation is allowed). If Ω and Q∞ is also constant then so is a. When
stepping in time it is therefore only necessery to update AW , AW 2 and µW
before solving the equations. A flowchart for the algorithm is seen in Figure 4.3.
Generate turbine geometry
Calculate Aµ , Aσ and a
Calculate AW and AW 2
Repeat
Solve for µ
Update structure (i.e. rotate turbine)
Update wake strength, µW , and shape
Finished
Evaluate forces
Figure 4.3: Flowchart for solving (4.1).
4.1 Freewake Model
The free wake is generated by letting vortices trail of the blade as the turbine
rotates, and then convect the vortices downstream according to an Euler time
marching scheme. For each time step the position of the vortex nodes are
updated as
∆X = Q∆t
Where Q is the total velocity, but not including velocities induced by the vortex
itself. This is the basic idea but it does not work very well without some mod-
ifications. The problems are related to stability, excessive computational time,
and that the shape of a fully developed wake is needed. The total simulation
time is set according to 8 revolutions of the rotor ttotal .
To enhance speed a wake roll up procedure has been implemented. At time tR1
the procedure begins and at time tR2 the wake consists of just a tip and a root
vortex. The root vortex has the strength of all the trailing vortices from the
root until a position just before the center of the blade. The tip vortex has the
strength of the trailing vortices from the center and out.
To enhance stability a number of improvements has been included. If a wake

node is past a specified x-position, xP , it is only moved in the x-direction
according to the u-velocity component. Thereby it is assumed that the wake has
reached its maximum expansion. The initial wake shape is not correct and at
simulation time tC all vortices older than a specified age, tC,age , is removed. A
prescribed wake of limited size is attached to the oldest vortices that otherwise
would have a free end. Without it the wake becomes unstable. The strength is
the same as that of the vortex to which it is attached. It has a helical shape
according to an assumed downstream velocity Qpre . Its downstream length is
Qpre tpre . A viscous core model has been implemented. The main effect of this
is to eliminate the very large velocities close to the wake vortices since this
makes the wake very unstable. The viscous core model is physical correct since
it models the decrease in velocity near the center of a wake vortex, and therefore
the development of the wake becomes reasonable. From a mathematical point
of view the viscous core model is not correct, and therefore care has been taken
to ensure that it has no effect near the rotorplane, since changes in induced
velocities directly affects the solution. The implemented viscous core model is
described in detail in Sant[13] and the equations are given here. The Biot-Savart
relation (2.6) is modified by a scaling factor, Kv
Γ dl × (r0 − r)
dq = Kv
4π |r0 − r|3
4.2 Evaluation of Forces 37
Kv is calculated from physical properties and the vortex age, t
h2
Kv =
(r4 + h4 )1/2
pc
rc = 4αδv ν(t + Sc ) , t >= -Sc
rc = 0 , t <= -Sc
Where h is the perpendicular distance from the vortex to the evaluation point.
α is a constant equal to 1.25643. ν is the kinematic viscosity and Sc is an age
offset constant. δv is a turbulent viscosity coefficient.
Table 4.1 summaries the parameters used in the simulation of a 30 meter turbine
at a tip speed ratio of λ=8. The resulting wake is seen in Figure 4.4.
tR1 =6 s tR2 =10 s Xp =8 m ttotal =18.85 s

tC =10 tC,age =6 Qpre =1/3 Q∞ tpre =4 s
Sc =-2.5 δv =100 ν=1.51e-5
Table 4.1: Freewake simulation parameters used in simulation of 30m rotor.

λ=8. Q∞ =10 m/s .
X (downstream u=w=0)
p
Q
∞
Prescribed wake
is attached here
Wake roll up
Figure 4.4: The developed wake.
4.2 Evaluation of Forces
The method for evaluation of forces depends on whether the panel or vortex
lattice method is used. A third method that can be applied to both of the
above mentioned, is to transfer data to a lifting line and then evaluate the
forces. The methods are described in the following.
4.2.1 Panel Method
In the panel method the pressures on the collocation points can be calculated
using the unsteady Bernoulli equation (2.11). The force on panel i of area Ai is
Fi = −pi · n Ai
The minus sign reflects that a positive pressure acts against the outward pointing
normal vector.
When using the Bernoulli equation it is nescessery to evaluate the velocities

induced by the doublet panels onto themself in a special way. The reason is
that the constant strength doublet panel mainly induces a velocity onto itself
in the normal direction and the tangential velocities on the collocation points
are therefore to small. This has only a small effect on the solution (i.e. the
doublet distribution) since it is based on the zero normal flow condition, but
it has a large effect on the pressure distribution. The solution is to assume a
more realistic distributed doublet distribution instead of the piecewise constant
distribution. The induced tangential velocities can then be evaluated using ana-
lytical formulas for the induced tangential velocities on doublet panels. If l and
m are chordwise and spanwise surface parameters, then the induced velocities
in the respective directions are
∂µ ∂µ
ql = − qm = −
∂l ∂m
4.2.2 Vortex Lattice Method
In the vortex lattice method the forces are obtained by applying the Kutta-
Joukowsky theorem on each bound vortex. The velocity used is the total velocity
Fi = ρQi × Γi
The best result is obtained if the velocities are evaluated in the collocation
points and then interpolated onto the bound vortices using appropriate surface
parameters. When this method is used an important force component, the
leading edge suction, is missing. Since this has a large influence on the power
production of turbines it was decided to use the method described next.
4.2.3 Lifting Line
The forces evaluated on a lifting line includes the leading edge suction, but any
forces due to velocities induced by the structure onto itself can not be evaluated.
4.3 Calculation of Wing Geometries for Vortex Lattice Simulations 39
This is not a serious problem since these forces are relatively small and confined
to the tip area. The major part of this force will be lift, and omitting it gives a
conservative estimate on the power production and the thrust. The evaluation
of forces on the lifting line is described in section 3.7.
To obtain the lifting line circulation distribution the strength of spanwise bound
vortices are summed over the chord at all spanwise positions. The spanwise
bound vortices are the vortices in the vortex lattice grid which are aligned with
the spanwise direction. The lifting line velocity is more subtle. The induced
velocity is induced by both the wake and the bound vortices, but in the lifting
line representation the velocities induced by the bound vortices must be omitted.
Since the wake induced velocities has large variations over the chord, it is not
clear which value to use. A solution is to recalculate the wake induced velocities
on the lifting line by extending the wake so it is attached directly to the lifting
line instead of the trailing edge of the 3 dimensional blade.
Table 4.2 compares calculated results for a flat AR=10, b=10 elliptic flat wing
in homogenious flow. The lift coefficients are very similar but the induced drag
coefficient calculated in the vortex lattice code, as described in section 4.2.2, is to
high. This is because of the missing leading edge suction. The forces evaluated
by transfering the bound circulation and velocities to the lifting line agrees very
well with the analytical results, and this method will be used henceforth.
CL CDi
Vortex Lattice 0.88289 0.15098
Liftin Line (Transfered Γb ) 0.88223 0.024819
Analytical (Based on Γ0 =max(Γb )) 0.88307 0.024822
Table 4.2: Comparison of methods for force evaluation.
4.3 Calculation of Wing Geometries for Vortex

Lattice Simulations
A major problem with wind turbine simulations is to predict the wake influence
in the rotorplane. Because of this it was eventually decided to focus on the
bound circulation distribution and the wake induced velocities in the lifting line
representation and not any chordwise distribution of properties. Therefore it is
not necessary to design a wing that represents the exact shape of the physical
wing, which will often be designed to satisfy a number of practical design issues,
but instead use a geometry which is more suitable from a computational point
of view. As long as the blade generates the correct circulation distribution in
the given velocity field its trailing wake will be correct and therefore also the
influence of it (assuming the wake model is correct).
The blade profiles used in the Vortex Lattice simulations are selected as flat with-
out camber because they are well defined and therefore eliminates one source of
uncertainty. At a given condition a flat blade can be defined in an infinity num-
ber of ways in order to generate a given lift, since both the angle of attack and
the chord can be varied, but it should be designed to generate correct results at
off design conditions.
The angle of attack of any profile is a function of the lift coefficient, and is
therefore dictated by the choice of it. In the linear range the relationship is
Cl = a0 (α + α0 ) , in general (physical blade)

Cl,f = 2π(α) , flat blade (vortex lattice)
Where Cl and Cl,f are the lift coefficients of the physical (real) blade and the
blade used in simulations, respectfully. The flat plate is given an offset angle
equal to α0 , because the ratio between Cl and Cl,f is then constant and equal to
a0 /2π at all inflow angles. This is important because it guarantees the correct
behavior of the flat blade at off-design conditions. In general Cl,f will be larger
than Cl since a0 < 2π, and the chord of the flat plate (cf ) should be adjusted so
that the lift generated at a given inflow angle is the same for both the physical
and the flat blade. I.e.
2π
Cl,f =Cl
a0
L 2π L
=
1/2ρQ2 cf a0 1/2ρQ2 c
a0
⇒ cf = c
2π
cf is the chord of the flat blade. Figure 4.5 shows an example of a physical
profile and the vortex lattice equivalent.
Figure 4.5: A realistic blade and a flat blade suitable for Vortex Lattice calcu-
lations. Both are at defined zero degree angle of attack
4.4 Grid Generation 41
4.4 Grid Generation
The panel grid has Nupper and Nlower panels on the upper and lower side of the
blade. If a grid for vortex lattice simulations are created Nupper is zero. Nspan
is the number of spanwise stations. The importance of the grid layout is well
described in the literature and will not be threated in detail here. The panels are
distributed according to a cosine spacing scheme in the spanwise direction on
both the upper and lower side, in order to group the panels closer together near
the leading edge/trailing edge and the tip/root. Figure 4.6 shows an example
of a grid layout. The grid is very coarse and in practice a finer mesh is used
with the discretization parameters given in table 4.3. This discretization was
selected because the results were reasonable and finer grids would lead to very
large computational times. xij refers to the truncation values defined in (3.1).
Figure 4.6: Example of a coarse vortex lattice grid.
# panels (root/trailing edge) (tip/leading edge)

Vortex lattice Nchord =25 xt1 =0.4 xt2 =0.4
Nspan =55 xt1 =0.05 xt2 =0.025
Table 4.3: Vortex lattice discretization parameters.
4.5 Reference Blade Data
The turbine blade data presented in this section is used for validation and as
reference values for a flat blade. The tip speed ratio is λ=8.
The calculated power and thrust coefficients and the distributed loads has been
verified against CFD calculations (see Ganuaa and Johansen[5]). Figure 4.7
shows the wake induced velocities and the bound circulation along the span. In
the following only the wake induced velocities will be stated because the other
velocity components are known. I.e. the kinematic velocity, (2.12).
0.5
0.125
0 0.1
Γ /(Q R) [−]
0.075
∞ B
−0.5 0.05
u/Q∞ [−]
v/Q [−] 0.025
∞
w/Q∞ [−]
−1 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
y/R [−] y/R [−]
Figure 4.7: The induced velocities and bound circulation.
Figure 4.8 shows the twist and chord distribution according to the Risø-B118
profile at 8 degrees angle of attack.
Twistz,tip= −2.2676o
100 0.16
0.14
80
0.12
60
0.1
Twistz [o]
c/R [−]
40 0.08
0.06
20
0.04
0
0.02
−20 0
0 0.2 0.4 0.6 0.8 1
y/R [−]
Figure 4.8: Twist and chord distribution.
The power and thrust coefficients are given in table 4.4.
Cp CT
0.513 0.866
Table 4.4: Power and thrust coefficient for flat reference blade.
4.6 Validation 43
4.6 Validation
The vortex lattice code is validated against the results for the flat blade pre-
sented in section 4.5 operating at a tip speed ratio of 8. The turbine was scaled
to a radius of 30 meter, and the free wake parameters given in table 4.1 was
used. The velocities and forces are found by representing data on the lifting line
as described in section 3.7.
Figure 4.9 and 4.10 compares the inviscous and the viscous forces to the design
goal (i.e. the input data). The inviscous forces are forces calculated using
the Joukowsky theorem and viscous forces are drag forces. The agreement is
relatively good but the forces are to large at center of the blade. Near the root
the forces has some error but this has very little effect on Cp .
Joukowsky forces
1
Design goal
Vortex Lattice
0.8
Section Cp [−]
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
y/R [−]
Figure 4.9: Section Cp due to the inviscous Joukowsky forces compared to design
values.
Viscous forces
0
−0.02
Section Cp [−]
−0.04
−0.06
−0.08
Design goal
Vortex Lattice
−0.1
0 0.2 0.4 0.6 0.8 1
y/R [−]
Figure 4.10: Section Cp due to Viscous forces compared to design values.
Figure 4.11 compares the bound circulation. The calculated values are to large
at the center, which explains the discrepancy found above.
0.15 Design goal

Calculated
0.1
Γ [m2/s]
B
0.05
0
0 0.2 0.4 0.6 0.8 1
y/R [−]
Figure 4.11: ΓB compared to design values.
In Figure 4.12 the inflow angles are compared. There is a large error near the
root, and over the center of the blade the inflow angle is to high (approximately
1/2 degree).
25
Design goal
Vortex Lattice
20
15
α [o]
10
0
0 0.2 0.4 0.6 0.8 1
y/R [−]
Figure 4.12: The local inflow angle compared to the design value.
The wake induced velocities are shown in figure 4.13. The data agrees reasonable
well over the center of the blade but some discrepancy is found near the tip and
root.
4.6 Validation 45
0.5
−0.5
Vortex Lattice: u/Q∞ [−]
Vortex Lattice: v/Q∞ [−]
Vortex Lattice: w/Q [−]
∞
−1
0 0.2 0.4 0.6 0.8 1
y/R [−]
Figure 4.13: Comparison of calculated induced velocities and design values

(dashed)
A summary of important data is found in table 4.5. Notice that the axial
induction is to small, which explains why the force and angle of attack is to high.
The error is probably in the wake model which does not accurately predict the
axial induction.
Design goal Vortex Lattice

Axial induction, a (y/R=0.5) 0.3314 0.3233 (-2.4%)
Cp 0.513 0.531 (+3.5%)
CT 0.866 0.887 (+2.4%)
Table 4.5: Comparison of important parameters
It is concluded that the vortex lattice code captures the physics of the flow but
with an error of 3.5% on the power. It can therefore not be expected that the
method can accurately predict turbine performance but it may be usefull for
studying variations in flow properties.
4.6.1 Validation of Panel Code
To validate the panel code a turbine rotating with constant rotational speed is
considered. The wake is completely omitted meaning that the spanwise undis-
turbed velocity distribution relative to the blade is known but the results are
not physically correct. The wake model will not be verified since it was in-
cluded in the vortex lattice calculations above. The undisturbed velocity seen
at a spanwise position is then the kinematic velocity
Qkin = Q∞ − Ω × r
Based on this a twist distribution can be calculated which corresponds to a

given angle of attack.
Calculations has been made on a turbine with a chord that varies linearly from
3 m at R=3 m to 1 m at the tip, R=30 m. The symmetric Van De Vooren profile
with 15% thickness and 12.o trailing edge angle is used because the analytical
pressure distribution is known (see Katz and Plotkin[11]). The pressure coeffi-
y
cients, calculated at a spanwise position of R =0.76 using the unsteady Bernoulli
equation is shown in Figure 4.14. The agreement with the analytical result is
good, considering that 3 dimensional effects can not be completely neglected.
1.5
1
−Cp [−]
0.5
−0.5
−1
0 0.2 0.4 0.6 0.8 1
x/c [−]
Figure 4.14: Comparison of analytical and numerically calculated pressure co-

y
efficients at spanwise position R =0.76 .
The panel code is accurate but the computational time is in the order of 12-24
hours. This is mainly because a very high number of panels is needed. Figure
4.15 shows lift coefficients against the number of panels on the upper and lower
side for an elliptic flat wing of aspect ratio 10 at 8o angle of attack. The wake
model is included. 100 panels on the upper and lower side is not enough to
ensure converged results. Figure 4.16 shows the distributed Cp values.
4.7 Conclusions 47
0.8
0.78
CL [−]
0.76
0.74
0.72
20 40 60 80 100
Nupper=Nlower [−]
Figure 4.15: Convergence of lift coefficient against chordwise discretization.
Cp
−1
−2
−3
Figure 4.16: Distributed Cp values.
4.7 Conclusions
The panel and vortex lattice codes have been validated. Because of uncertainties
in the free wake model and the use of relatively coarse discretization due to
excessive computational time, the accuracy is limited to within ≈3-5%. This is
relatively accurate but the uncertainty is larger than the expected increase in
Cp , and it is questionable whether the effects of wingles can be predicted using
the methods described in this part.
Part 5
Blade and Winglet Design
This part issues the design of winglets. Vortex lattice simulations of a poor
design are not interesting and the need for a design algorithm is apparent. A
method will be developed which is found to be fast and numerically consistent.
The method is based on input in the form of a known bound circulation and
velocity distribution of a flat blade. The input velocity should be the total veloc-
ity relative to the blade, but not including velocities induced by the blade itself
since a lifting line representation is used. This means that CFD and experimen-
tal results can not be implemented easily, but the blade element momentum
method can be used to generate input data effectively because the induction
factors corresponds to the wake induced velocities on a flat blade. In the fol-
lowing, data for the flat reference blade presented in section 4.5 will be used as
input.
The design method is based on the assumption that the winglet is relatively
small and that only the circulation distribution on the outermost part of the
blade is altered when winglets are attached to the flat reference blade. In prac-
tice the winglet geometry and the bound circulation will be defined and the
velocities on the new geometry will be calculated. The definition of the bound
circulation and the calculation of the new velocity distribution is described in
detail in the following sections.
50 Blade and Winglet Design
5.1 Defining Bound Circulation
The bound circulation used as input is defined over spanwise positions of the
blade in the range y/R ∈ [0, 1]. On the blade with winglet this will be al-
tered and extended onto the winglet using an appropriate scheme. If the design
method described here is to be used it is necessary to keep the bound circu-
lation constant until a position close to the tip (=y0 ). On physical grounds
the circulation-value, Γ0 , and the derivative , Γ′0 , must be continuous at this
point, because the force distribution is expected to vary smoothly. It will also
be required that the circulation approaches zero at the winglet tip. The defini-
tions are shown in Figure 5.1.a . Since the y-coordinate is not suited to describe
non-flat blades the spanwise surface parameter s is used. s0 corresponds to the
y0 /R position and s >1 corresponds to positions on the winglet. Figure 5.1.b
shows an example of a circulation distribution extended from s0 /R=0.91 onto
a 3% winglet. A second parameter st ∈ [0, stip ] has also been defined, which
descripes the part of the blade where the bound circulation is defined
st = s − s0
slope = Γ0’ st
B
Γ
s0
Γ
Γ0
ΓB, input (without winglets)

Γ (with winglets)
B
s 0.5 0.6 0.7 0.8 0.9 1 1.1

s/R [−]
(a) Definitions of circulation parame- (b) Example of ΓB with and without

ters. winglets.
Figure 5.1: Circulation distributions.
A systematic choice of the tip circulation Γtip (st ) is described in the following..
The tip circulation must satisfy the 3 physical conditions.
1 The derivative at st =0 must be continuous.
dΓtip dΓB
=
dst st =0 dst
5.1 Defining Bound Circulation 51
2 The value at st =0 must be Γ0 .
Γtip |st =0 = Γ0
3 The circulation decays to zero at the tip.
Γtip |st =stip = 0
Suppose a function f (st ) has been defined that satisfies these conditions. In
order to be able to shape it it is multiplied by a function g(st )
Γtip (st ) = f (st )g(st ) (5.1)
Since f already fulfills the conditions on the tip circulation, g must satisfy the
following
1 The derivative at s=0 must be zero.

dg
=0
ds s=0
2 The function value at s=0 must be unity.
g|s=0 = 1
3 The function value at st =stip is arbitrary but must be finite.
Once a basis function f has been defined any function g, that satisfies the above
can be used to shape the tip circulation distribution. f and g has been selected
in order to resemble the results of section 3.5 where the optimum circulation
distribution on wings was found. More specifically f and g has been selected as
f (st ) = as2t + bst + c

2
g(st ) = eαst
α is a shaping parameter. The constants are given as

bstip + c dΓB
a=− b= c = Γ0
s2tip dst st =0
An example of the resulting circulation distribution is shown in Figure 5.2. Note

that negative values of α pushes the circulation to the left, and positive values
to the right onto the winglet. I.e. α represents how heavy the winglet is loaded
and larger values corresponds to a heavier loading.
0.1
α = 100
0.08
α=0
[m /s]
0.06
2
tip
Γ 0.04
0.02 α = −150
0
0 0.05 0.1
s [m]
t
Figure 5.2: Examples of tip circulations against α -values.
5.2 Calculating Velocities
The velocities on the blade with winglet are needed in order to calculate the
forces. The velocites are obtained by decomposing the input velocities on the
flat blade into components that can be transfered onto the blade with winglet,
taking into account the altered circulation distribution in the tip region. The
wake is split into two parts, the nearwake and the farwake. The nearwake is
the wake trailing of the outer part of the blade where the bound circulation
is altered and one quarter revolution downstream. Figure 5.3 illustrates this.
The farwake denotes all other wake singularites which includes the wake trailing
of the inner part of the blade. The velocities induced by the farwake can be
considered constant on the blade. To justify this, first consider the part of the
farwake immediately behind the inner part of the blade. The wake strength is
unchanged and assuming only small changes in the wake geometry, the induced
velocities over the blade are also constant. Next consider the rest of the wake
from the one quarter revolution and downstream. When the bound circulation
near the tip is altered then so is the strength of the outermost part of the wake,
but since the wake is a distance away from the blade the induced velocity on
the blade is largely unaffected by this redistribution of wake strength.
5.2 Calculating Velocities 53
Tip with
Blade winglet
Q
∞
Nearwake
Figure 5.3: The nearwake extending one quarter revolution downstream.
The velocities at a given point on the blade is now assumed to consist of a

constant part, QF F , due to farwake induced and freestream velocities, and a
part, QN W , which is induced by the near wake and which depends on the
distribution of bound circulation near the tip. QN W can be calculated based
on the known strength and an assumed geometrical shape. If a helical shape is
assumed the wake trailing of the position x={0, y, 0}T is described by
xN W = {Q∞ (1 − a), y cos tΩ, −y sin tΩ}T , t ∈ [0, π/(2Ω)]
QF F can now be obtained on the main blade (i.e. not on the winglet) by
subtracting QN W from the total velocity
Blade: QF F = Q − QN W
QN W should correspond to the input data of the flat blade. The winglet will
be mounted along a line extending from the tip downstream. The velocity on
this line is unknown except on the tip, Qtip , and can therefore be written as
Q = Qtip + ∆QF F + ∆QN W
Where ∆QF F and ∆QN W denotes changes in QF F and QN W relative to Qtip .
Even though QF F is constant on the blade this is not necessarily the case on
the line where the winglet is mounted. Specifically, since the distance to the
nearest downstream tip-vortices is relatively small when compared to the length
of the winglet, the velocities induced thereby is expected to vary over the length
of the winglet. ∆QF F , is modelled in section 5.2.1, and it is found to have a
v-component. ∆QN W can be calculated numerically as
∆QN W = QN W − QN W |tip
Subtracting QN W from Q yields the farfield induced velocities on the winglet
Winglet: QF F = Qtip + ∆QF F − QN W |tip
The constant velocity component QF F is now known over the blade with winglet.
To obtain the total velocity it is necessary to add the nearwake induced velocity
QN W , but this should now be evaluated based on the defined bound circulation.
5.2.1 Modelling ∆QF F
To model the change in the farfield velocity along the winglet the wake singular-
ities that are relatively close to it must be considered. This is mainly the wake
trailing of the inner part of the blade and the downstream wake directly behind
the blade tip. The former is parallel to the winglet and with constant distance,
meaning that the induced velocities are relatively constant. The latter is ex-
pected to vary significantly over the winglet since it is directed perpendicular to
it and the distance to it varies significantly over the winglet. In the following the
wake downstream is assumed to consists of concentrated tip and root vortices.
The strength of both is half the maximum value of the bound circulation
Γ = max(ΓB )
Instead of simply modelling a large part of the wake numerically it will be shown
that it is only necessary to take into account a small part of the wake. The final
equation is written as a sum consisting of few terms, which can be included in
analytical models. Figure 5.4 illustrates the important downstream tipvortices.
Direction of rotation
x
vi (x)
nmax ∆x Γ
Γ
Γ ∆x
Figure 5.4:
The distance from the tip to the n’th downstream tipvortex is
2/3π
∆xn = Q∞ (1 − a) n
| {z Ω}
∆x
5.2 Calculating Velocities 55
The ratio of turbine circumference to ∆x is

∆xn 1/3(1 − a)
= n
2πR λ
For a tip speed ratio of λ=8 the ratio is 1/28. I.e. the first tip vortex is
1/28 of the turbine circumference downstream. This is so close to the turbine
that it must be expected that the induced velocity varies over the winglet.
Therefore the influence of nmax downstream tipvortices must be included. To
determine nmax a vortex segment of a tipvortex directly downstream of the tip
is considered. The induced velocity on the winglet is in the positive y-direction
and of magnitude
Γ 1
dQi (r) = ds
4π r2
The ratio of induction at the tip and the root of the winglet is (assuming an
axial induction of 1/3)
dQi (∆xn − h) 1
α= = h
dQi (∆xn ) 1− 4Rnπ/(9λ)
Specifying α and solving for n yields nmax . E.g. α=1.10 ⇒ nmax =4, i.e. by
including 4 downstream tipvortices, the induced velocity over the winglet by
any vortex constituting the farwake, will vary less or equal to 10%. Since most
of the tip-vortices are at large distances from the winglet, only a part of it with
a length equal to the distance downstream of the nmax ’th tip vortex, will be
modelled. This is 11% of the rotor circumference and the circular shape can be
approximated as a straight line vortex. The induced velocity of such a straight
tip vortex segment of length ∆x nmax on the winglet can be written using (2.7)
as

Γ −1 2(n − x/∆x)
vi (x, n) = cos tan
2π(∆x n − x) nmax
The velocity is in the positive outward radial direction. The induced velocity
due to the nmax tip vortex segments is
nX
max
vi (x) = vi (x, n)
n=1
Finally, the change in the farfield induced velocity relative to the velocity at the
tip is
∆QF F = {0, vi (x) − vi (0), 0}T
Figure 5.5 shows ∆QF F over a 5% winglet. The velocity increases to 6% of Q∞ ,

which is not substantial but can be expected to have an influence.
0.06
0.04
v/Q [−]
∞
0.02
0
0 0.01 0.02 0.03 0.04
x/R [−]
Figure 5.5: The v-component of ∆QF F .
5.3 Summary of Method
The velocity over the blade with winglet can now be found by calculating the
velocities induced by the nearwake on all blade and winglet stations and adding
it to the farfield velocity field. The forces are evaluated on the lifting line using
the method described in section 3.7. The complete algorithm is described below.
• Select a point close to the tip, y0 , and determine Γ0 and ΓB .
• Define a tip geometry based on winget height, h/R.
• Calculate the farfield velocities on the blade and winglet, QF F .
• Assume a tip circulation distribution.
• Calculate nearwake strength and shape.
• Calculate nearwake induced velocities on blade and winglet, QN W .
• Calculate total velocity on blade and winglet, Q=QN W + QF F .
• Evaluate forces based on the total velocity and defined bound circulation
distribution.
• Define twist and chord according to velocity and bound circulation distri-
bution.
5.4 Validation 57
5.4 Validation
The method can not be fully evaluated because no data for comparison is at
hand. Here it will be shown that the method is consistent and the sensitivity of
important parameters will be discussed.
5.4.1 Convergence with respect to discretization.
nc denotes the number of collocation points on the tip. They are distributed
according to a cosine spacing scheme (equation (3.1)) in the range st ∈[0, stip ]. A
parameter study has been made using data from table 5.1. Figure 5.6 illustrates
how the CP -values converges when plotted against 1/nc . An infinitely large
discretization is equivalent to 1/nc =0 and the assumed exact result is read of
as CP =0.5269. Comparing this with the results it is found that the results are
within 99.90% of the exact value for nc ≥80. nqw is the number of vortices in
óne trailing nearwake vortex. The results depends very little on this number
and a low value is sufficient.
nqw =20 nc ∈[30, 100] xt1 =π/2 xt2 =0

h/R=0.032 Rc /h=0.10 y0 /R=0.9085 Γ0 /(Q∞ R)=0.1068
Cl /Cd = 110 a=1/3 α=83.33
Table 5.1: Parameters used for studying variation of nc .
0.527
0.5265
CP [−]
0.526
0.5255
0.525
0 0.01 0.02 0.03
1/nc [−]
Figure 5.6: Convergence of CP against 1/nc .

5.4.2 Influence of assumed axial induction, a.
The influence of the choice of a is not obvious, since it affects both the inflow
velocity and the shape of the wake. A parameter study using data from table
5.2 shows that the results are largely unaffected of varying its value. Figure 5.7
shows how the Cp values depends on a.
nqw =20 nc =80 xt1 =π/2 xt2 =0

h/R=0.032 Rc /h=0.10 y0 /R=0.9085 Γ0 /(Q∞ R)=0.1068
Cl /Cd = 110 a ∈[1/4, 1/2] α=83.33
Table 5.2: Parameters used for studying variation of a.
0.527
0.5265
C [−]
0.526
P
0.5255
0.525
0.2 0.3 0.4 0.5 0.6
a [−]
Figure 5.7: Dependency on variation of a.
5.4.3 Influence of choice of y0 .
The choice of y0 influences many parameters and therefore it is hard to study

the effect of varying it. Selecting a value very close to the tip means that
the circulation that can be distributed onto the winglet is very small, hence
decreasing the potential effect of it. The optimum choice would be to select a
value where the tip effects are not yet dominating the circulation distribution.
In practice this means that Γ0 should have a relatively large value and that
Γ0 should be relatively small, but the freedom of choice is still limited by the
condition that the farwake must be considered unaffected. In practice it is found
that good results can be obtained by setting y0 ≈0.85-0.90.
Figure 5.8 illustrates the problem. For different values of y0 /R the optimal value
of α and therefore Cp is found. The Cp value apparently increases when y0 /R
decreases but the results are uncertain because the farfield of the wake can then
no longer be assumed to be unaffected.
5.4 Validation 59
0.55
0.54
CP [−]
0.53
0.52
0.51
0.8 0.85 0.9 0.95 1
y0/R [−]
Figure 5.8: The optimized Cp values for different choices of y0 /R .
5.4.4 Influence of Curve Radius Rc /h .
In this study the optimal value of α is found for each value of the curve radius
Rc /h. Parameters are given in table 5.3. Figure 5.9 shows how the largest Cp
values are found for small curve radii. This result is consistent with the results
for the optimum loaded wing given in section 3.5.
nqw =20 nc =80 xt1 =π/2 xt2 =0

h/R=0.032 Rc /h ∈[0.05, 0.50] y0 /R=0.90 Γ0 /(Q∞ R)=0.1068
Cl /Cd = 110 a=1/3 α=α(Rc /h)
Table 5.3: Parameters used for studying variation of Rc /h.
0.5264
0.5262
0.526
0.5258
CP [−]
0.5256
0.5254
0.5252
0.525
0 0.1 0.2 0.3 0.4 0.5
Rc/h [−]
Figure 5.9: Influence of varying curve radius, Rc /h .

5.5 General Design Results
In this section the influence of the shaping parameter α and the winglet height
are studied. Based on the parameters in table 5.4 the power and thrust coeffi-
cients has been calculated for 40 values of h/R and α.
nqw =20 nc =80 xt1 =π/2 xt2 =0

h/R ∈[0.01, 0.07] Rc /h=0.10 y0 /R=0.8878 Γ0 /(Q∞ R)=0.1068
Cl /Cd = 110 a=1/3 α ∈[-50, 200]
Table 5.4: Table of parameters.
Figure 5.10 shows a contour plot of the Cp values. The maximum Cp is found
for a winglet height of approximately 3.5% and α=75. If α is increased further
the Cp values decreases fast. This is because the viscous drag associated with
the heavier loading of the winglet starts to dominate. For larger winglets this
happens for smaller values of α.
C [−]
p
200
0.52
150
0.51
0.5
100
α [−]
0.49
50 0.48
0.47
0
0.46
−50 0.45
0.01 0.02 0.03 0.04 0.05 0.06 0.07
h/R [−]
Figure 5.10: Calculated power coefficients.
Figure 5.11 shows only the range of Cp values larger than the values for a flat
blade. Note that only a band of α values yields an increase in Cp , and therefore
that only well designed winglets contributes positively.
5.5 General Design Results 61
C [−]
p
200
0.528
150 0.526
0.524
100
0.522
α [−]
0.52
50
0.518
0 0.516
0.514
−50
0.01 0.02 0.03 0.04 0.05 0.06 0.07
h/R [−]
Figure 5.11: The range of increased Cp relative to flat wing
The thrust coefficients are shown in Figure 5.12. As expected they increases
with increasing α values, because the loading on the winglet and the outer part
of the blade increases.
C [−]
T
200
1.1
150
1.05
100
1
α [−]
50 0.95
0 0.9
0.85
−50
0.01 0.02 0.03 0.04 0.05 0.06 0.07
h/R [−]
Figure 5.12: Calculated thrust coefficients
A large increase in thrust is unwanted and therefore is it not possible to point out
the best design. Figure 5.13 shows a contour plot of the ratios of power to thrust
(Cp /CT ). A large value is desired but this is not possible without selecting a
Cp lower than the maximum possible. On the figure is also shown with black
the contour lines corresponding to Cp values of 0.52 and 0.528 respectfully.
They show that choosing the optimum Cp value (Cp =0.528) corresponds to a
low value of (Cp /CT ), and therefore a large CT . Choosing Cp =0.52 instead is
possible without an associated large increase in thrust. Two sets of h/R and
α values has been selected based on this and these designs will be threated in
detail in part 6. They are denoted optimum and conservative design respectfully.
Table 5.13 summaries their parameters. Note that a third, interesting design is
found for a winglet height of h/R=0.02 and α=50. This winglet is small with
a relatively high Cp value, and an increase in thrust which lies between the
conservative and the optimum design values. Due to time limitations it will not
be threated further.
C /C [−]
p T
200
Optimum design:
C = 0.528
p 0.59
150
Cons. design:
C = 0.52
p 0.58
100
α [−]
50 0.57
0 0.56
−50 0.55
0.02 0.04 0.06
h/R [−]
Figure 5.13: Ratio of Cp to CT and the choices of conservative and optimum

design.
Conservative design Cp =0.521 CT =0.871

α=7.69 h/R= 0.048
Optimum Cp design Cp =0.528 CT =0.892
α=52.56 h/R=0.035
Table 5.5: Design parameters for defined conservative and optimum design.
5.6 Conclusions
Using the described design algorithm important general design conclusions has
been established. Especially how the winglet height and the bound circulation
influences the overall design. Based on this two designs has been defined which
will be studied in part 6.
The design algorithm is fast. One evaluation of a design takes approximately

1 second in Matlab using a fully numerical code which is not optimized. More
efficient coding in C or Fortran can substantially decrease the computational
time, and this makes the model fast enough for implementation in an unsteady
BEM algorithm.
Part 6
Final Designs
A conservative and an optimum winglet design is described in detail. The design

parameters are based on the discussion in section 5.5. Unless otherwise stated
the results are calculated using the design algorithm described in part 5 using
the flat blade reference data of section 4.5 as input.
6.1 Conservative Winglet Design
Design parameters are given in table 6.1.
nqw =20 nc ∈80 xt1 =π/2 xt2 =0

h/R=0.048 Rc /h=0.10 y0 /R=0.8878 Γ0 /(Q∞ R)=0.11083
Cl /Cd = 110 a=1/3 α=7.69 λ=8
Table 6.1: Design parameters for 4.8% conservative winglet design.

64 Final Designs
6.1.1 Overall Design Results
Table 6.2 summaries the geometry and performance parameters. Notice that
the Cp value is increased by 1.60% for a 0.65% increase in thrust relative to a
flat blade.
h/R Rc /h Cp CT
0.048 0.10 0.521 (+ 1.60 %) 0.871 (+ 0.65%)
Table 6.2: Final design results
The calculated wake induced velocities and the bound circulation distribution
over the blade is seen in figure 6.1.a and 6.1.b. Note the large negative v-
velocity component, in the direction from the winglet towards the hub. This is
desirable because it decreases the normalwash on the winglet. The circulation
is distributed smoothly onto the winglet and differs from the flat blade only in
a region very close to the winglet. This shape is physically realistic and because
the changes are confined to a very small region the assumptions in the design
algorithm are expected to be fullfiled.
0.5
0.12
0.1
0
Γ /(Q R) [−]
0.08
∞
0.06
B
−0.5
u/Q [−] 0.04
∞
v/Q [−]
∞ Γ (input, without winglets)
0.02 B
w/Q [−]
∞ Γ (with winglets)
B
−1 0
0.8 0.85 0.9 0.95 1 1.05 0.6 0.8 1
s/R [−] s/R [−]
(a) Wake induced velocities (b) Bound circulation
Figure 6.1: Induced velocities and bound circulation.
The twist and chord distribution on the outer part of the blade and the winglet
is shown in Figure 6.2 and the blade itself is seen in figure 6.3. The blade is
designed using the Risø-B118 profile at 8 degrees angle of attack.
6.1 Conservative Winglet Design 65
4 0.05
2
0.04
0
Twist [ ] −2 0.03
o
c/R [−]
z
−4
−6 0.02
−8
0.01
−10
−12 0
0.4 0.6 0.8 1 1.2
s/R [−]
Figure 6.2: Twist and chord distribution on outer part of blade.
Figure 6.3: The final blade with winglet.
6.1.2 Design Predicted Performance
The section Cp values over the blade are shown in Figure 6.4. The inviscous
(Joukowsky) forces on the winglet contributes with negative values and the
winglet is therefore being dragged by the main blade. The winglet is relatively
slim and is lightly loaded on the outer part. Because it is desirable to have a low
winglet height this design can probably be improved by using a smaller winglet
with heavier loading.
66 Final Designs
1 0
input (without winglets)
0.8 with winglets
−0.02
0.6
Section C [−]
Section C [−]
−0.04
p
p
0.4
0.2 −0.06
0
−0.08
−0.2 input (without winglets)
with winglets
−0.4 −0.1
0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1
s/R [−] s/R [−]
(a) Joukowsky forces. (b) Viscous forces.
Figure 6.4: Section Cp due to Joukowsky (lift) and Viscous forces.
The section thrust coefficient distribution is shown in figure 6.5. There is only
a very small increase in thrust forces near the tip. Notice that the winglet itself
is not subject to thrust forces.
1.5
Section CT [−]
0.5
with winglets
0
0.5 0.6 0.7 0.8 0.9 1
s/R [−]
Figure 6.5: Section thrust coefficients
The normalwash on wind turbines, Qn , is defined in Figure 6.6 as the component

of the total velocity perpendicular to both the blade lifting line and the direction
of rotation. It is defined positive in the direction from the suction to the pressure
side, and is the velocity component responsible for the generation of torque
according to the Joukowsky theorem. Notice that positive values generates
negative torque, and it is therefore desirable that Qn is decreased relative to the
flat blade.
6.1 Conservative Winglet Design 67
Direction of rotation
y
Qn
Figure 6.6: Definition of normalwash on a wind turbine blade.
Figure 6.7 shows the calculated normalwash. On the main blade the normalwash
has been decreased near the tip, which is important because it decreases the
induced drag. The normalwash on the winglet is positive and it yields a negative
contribution to the power production.
0.4
0.2 with winglets
0
Q /Q [−]
−0.2
∞
−0.4
n
−0.6
−0.8
−1
0.5 0.6 0.7 0.8 0.9 1
s/R [−]
Figure 6.7: Normalwash velocity magnitude.
6.1.3 Vortex Lattice Simulation Results
Using the free wake vortex lattice code the power and thrust coefficients are
Cp =0.534 and CT =0.916. This is an increase of 0.6% and 3.8% respectively,
relative to the vortex lattice results for the flat blade given in table 4.5. This is
not in agreement with the results of the design algorithm stated above. Figure
6.8.a compares the calculated bound circulation with the design goal. There is
some agreement but the relative error is large. Figure 6.8.b shows the section
Cp values due to the Joukowsky forces on the blade. The agreement is not
good. One source of error is the relatively coarse grid near the tip of the blade
which can not resolve the distribution of properties accurately. Due to excessive
computational time it is not possible to increase it and it was therefore decided
to make no further vortex lattice simulations.
68 Final Designs
0.15 Joukowsky forces

1.5
Section C [−]
0.1
Γ /(Q R) [−]
p
0.5
∞
0.05
B
−0.5 Design goal

0 Design goal Vortex Lattice
Calculated −1
0.6 0.8 1 0.6 0.8 1
s/R [−] y/R [−]
(a) The bound circulation distribu- (b) The section Cp due to Joukowsky
tion. forces.
Figure 6.8: Vortex lattice simulation results compared to design values.
6.2 Optimum Winglet Design
Design parameters are given in table 6.3.
nqw =20 nc ∈80 xt1 =π/2 xt2 =0

h/R=0.035 Rc /h=0.10 y0 /R=0.8878 Γ0 /(Q∞ R)=0.1108
Cl /Cd = 110 a=1/3 α=52.56 λ=8
Table 6.3: Design parameters for 3.5% optimum winglet design.
6.2.1 Overall Design Results
Table 6.9 summaries the geometry and performance parameters. Notice that
the Cp value is increased by 2.9% but the increase in thrust relative to a flat
blade is now relatively large, 3.1%.
h/R Rc /h Cp CT
0.035 0.10 0.528 (+ 2.9 %) 0.892 (+ 3.1%)
Table 6.4: Design parameters for 3.5% optimum winglet design.
The calculated wake induced velocities and the circulation distribution over the
blade is seen in figure 6.9.a and 6.1.b. The circulation differs substantially from
the conservative design, i.e. the outer part of the blade and the winglet is now
heavily loaded.
6.2 Optimum Winglet Design 69
0.5 0.12
0.1
0
Γ /(Q R) [−]
0.08
−0.5
∞
0.06
B
−1 u/Q [−] 0.04
∞
v/Q [−]
∞ Γ (input, without winglets)
−1.5 0.02 B
w/Q [−]
∞ Γ (with winglets)
B
0
0.8 0.85 0.9 0.95 1 1.05 0.6 0.8 1
s/R [−] s/R [−]
(a) Wake induced velocities (b) Bound circulation
Figure 6.9: Induced velocities and bound circulation.
The twist and chord distribution on the outer part of the blade and the winglet
is shown in Figure 6.10 and the blade itself is seen in figure 6.11. The blade is
designed using the Risø-B118 profile at 8 degrees angle of attack.
6 0.05
2 0.04
0
0.03
Twist [ ]
o
c/R [−]
−2
z
−4
0.02
−6
−8 0.01
−10
−12 0
0.4 0.6 0.8 1 1.2
s/R [−]
Figure 6.10: Twist and chord distribution on outer part of blade.

70 Final Designs
Figure 6.11: The final blade with winglet.
6.2.2 Design Predicted Performance
The section Cp values over the blade are shown in Figure 6.4. The winglet is be-
ing dragged by the main blade and the viscous forces are increased substantially
due to the heavier loading of the outer part of the blade and the winglet.
1.5 0
with winglets
−0.02
1
Section C [−]
Section C [−]
−0.04
0.5
p
−0.06
0
−0.08
−0.5
−0.1
with winglets
−1 −0.12
0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1
s/R [−] s/R [−]
(a) Joukowsky forces. (b) Viscous forces.
Figure 6.12: Section Cp due to Joukowsky (lift) and Viscous forces.
The section thrust coefficient distribution is shown in figure 6.13. There is now
a large increase in thrust near the tip, and this indicates that the increase in Cp
is caused mainly by the large forces near the tip.
6.3 Conclusions 71
1.5
Section CT [−]
1
0.5
with winglets
0
0.5 0.6 0.7 0.8 0.9 1
s/R [−]
Figure 6.13: Section thrust coefficients
Figure 6.14 shows the normalwash. The reduction is small compared to the
conservative design, and this confirms that the large increase in Cp is due to the
larger forces near the tip.
0.5
with winglets
0
Q /Q [−]
−0.5
∞ n
−1
−1.5
−2
0.5 0.6 0.7 0.8 0.9 1
s/R [−]
Figure 6.14: Normalwash velocity magnitude.
6.3 Conclusions
The distributed properties calculated using the design algorithm are physically
reasonable and indicates very well how the winglet affects the tip region. The
vortex lattice results are not in agreement with the design algorithm. This can
be caused by the numerous problems already described, or it can be caused by
the poor discretization of the winglet. It is also a possibility that the design al-
gorithm is not correct, and at this point it is desirable to make a CFD validation
before progressing further.
The physics of the two designs confirms the general design conclusions of section
5.5.
72 Final Designs
Part 7
Conclusion
The implemented vortex lattice and panel code were not satisfactory. Among
other issues the problems were related to uncertainties in the free wake model
and excessive computational time. The latter implicated that a relatively coarse
discretization had to be used which further added to the uncertainty of the
results.
The developed design method is fast and consistent. It can predict the change in
flow properties near the tip and the associated forces introduced by the winglet.
The method have not yet been validated, something which the vortex lattice
code was intended for, but the results are physically reasonable and comparable
to other results reported in the literature. The design method is based on the
decomposition of the total velocity into a constant component and a nearwake
induced component. This approach appears promising and can form the basis
for aeroelastic computations.
Winglets was found to affect only the outermost part of a wind turbine blade
(y/R >≈ 0.85) indicating that it affects the flow locally. The induced drag
component is decreased in this region but the thrust forces are increased. A
general design study has shown that it is possible to design for a large increase
in power with a relatively small associated increase in thrust forces relative to a
flat blade of the same radius. I.e. an increase in Cp of 1.65% for an increase in
CT of 0.65%. The power can be increased by as much as 2.9% but the increase
74 Conclusion
in CT is then 3.1%. It was further shown that winglets of a given height can
have both a positive and a negative effect depending on the distributed loading.
I.e. winglets has to be well designed. Winglet with heights in the range of 2-4%
of the rotor radius appears to be most effective.
Appendix A
Appendix
A.1 Constant Strength Source Panel (3 Nodes)
A solution to (2.8) in closed form for both the induced potential and the velocity
of a 4 node constant strength sourcepanel is given in Katz[11]. Often a 3 point
panel is required and the transformation to a 3 point panel is described below.
The transformation is based on equations (10.95) - (10.97) in Katz[11] which
are not stated here. The final set of equations will be stated below. Notice that
the panel is entirely in a x-y plane and therefore a transformation of coordinates
and of the resulting induced velocity is required in practical applications.
76 Appendix
A.1.0.1 Transformation of 4 Node Constant Strenght Source Panel

Into 3 Node Panel
set (x3 , y3 )=(x4 , y4 ). This yields
d34 = 0
d41 = d31
m41 = m31
r3 = r4
e3 = e4
h3 = h4
Notice that m34 is undefined, but all terms including it cancels since

−1 m34 e3 − h3 −1 m34 e4 − h4
tan = tan
zr3 zr4
It remains to find the limit of the following 2 terms as (x3 , y3 ) → (x4 , y4 )

y4 − y3 r3 + r4 − d34
ln (A.1)
d34 r3 + r4 + d34

x3 − x4 r3 + r4 − d34
ln (A.2)
d34 r3 + r4 + d34
Substituting d34 into (A.1) and rewriting it yields

y4 − y3 r3 + r4 − d34
p ln
2
(x4 − x3 ) + (y4 − y3 ) 2 r3 + r4 + d34
s
1 r3 + r4 − d34
= (x −x )2 ln
r3 + r4 + d34
(y4 −y3 )2 + 1 |
4 3
{z }
| {z } →ln(1)=0
∈[1,0+ ]
The argument to the logarithm is 1 in the limit since d34 goes to zero and
2
r3 =r4 is finite. In the limit (x4 −x3 )
(y4 −y3 )2 can take the values [0, +∞]. Therefore the
square root is limited and takes values in the range [1, 0+ ]. Muliplying this with
ln(1)=0 gives the result that the whole term is zero. By the same argumentation
(A.2) is equal to zero.
A.1 Constant Strength Source Panel (3 Nodes) 77
A.1.0.2 Velocity Induced by a 3 Node Constant Strength Source

Panel
The equations describing the velocity induced at point (x, y, z) by a 3 point

source panel of strength σ (Figure 2.6), is stated below. Notice that the panel
must lie in the x-y plane. I.e. the nodes must have the coordinates (x1 , y1 , 0),
(x2 , y2 , 0) and (x3 , y3 , 0).
p
d12 = (x2 − x1 )2 + (y2 − y1)2
p
d23 = (x3 − x2 )2 + (y3 − y2)2
p
d31 = (x1 − x3 )2 + (y1 − y3)2
y2 − y1
m12 =
x2 − x1
y3 − y2
m23 =
x3 − x2
y1 − y3
m31 =
x1 − x3
p
rk = (x − xk )2 + (y − yk )2 + z 2 , k = 1, 2, 3
ek = (x − xk )2 + z 2 , k = 1, 2, 3
hk = (x − xk )(y − yk ) , k = 1, 2, 3
The induced velocity components:

σ y2 − y1 r1 + r2 − d12 y3 − y2 r2 + r3 − d23 y1 − y3 r3 + r1 − d31
u= ln + ln + ln
4π d12 r1 + r2 + d12 d23 r2 + r3 + d23 d31 r3 + r1 + d31

σ x1 − x2 r1 + r2 − d12 x2 − x3 r2 + r3 − d23 x3 − x1 r3 + r1 − d31
v= ln + ln + ln
4π d12 r1 + r2 + d12 d23 r2 + r3 + d23 d31 r3 + r1 + d31

σ m e
12 1 − h 1 m e
12 2 − h 2
w= tan−1 − tan−1 +
4π zr1 zr2

m23 e2 − h2 m23 e3 − h3
tan−1 − tan−1 +
zr2 zr3

m31 e3 − h3 m31 e1 − h1
tan−1 − tan−1
zr3 zr1
78 Appendix
Bibliography
[1] John D. Anderson. Fundamentals of Aerodynamics. 2005.

[2] Jr Clarence D. Cone. The theory of induced lift and minimum induced
drag of nonplanar lifting systems, nasa r-139. 1962.
[3] Eppler. Induced drag and winglets. 1995.
[4] Maughmer et. al. The design and testing of a winglet airfoil for low speed
aircraft. 2001.
[5] Mac Gaunaa and Jeppe Johansen. Can Winglets Increase Cp. 2007.
[6] Stig Øje. Instationære aerodynamiske kræfter på todimensionalt vingepro-
fil. 1980.
[7] Kroo. Drag due to lift. concepts for prediction and reduction. 2000.
[8] Jones & Lasinsky. Effect of winglets on the induced drag of ideal wing
shapes. 1980.
[9] Maughmer. The design of winglets for high-performance sailplanes. 2001.
[10] Jack Moran. Theoretical and Computational Aerodynamics. 1984.
[11] Katz & Plotkin. Low Speed Aerodynamics. 2005.
[12] Jeppe Johansen & Niels N. Sørensen. Aerodynamic investigation of winglets
on wind turbine blades using cfd. 2006.
[13] Gijs Van Kuik Tonio Sant and G. J. W. van Bussel. Estimating the angle
of attack from blade pressure measurements on nrel phase vi rotor using a
free wake vortex model: Axial conditions. 2006.
Risø’s research is aimed at solving concrete
problems in the society.
Research targets are set through continuous

dialogue with business, the political system and
researchers.
The effects of our research are sustainable energy

supply and new technology for the health
sector.
www.risoe.dk

Modeling

Uploaded by

Copyright:

Available Formats

Modeling

Uploaded by

Document Information

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Modeling

Uploaded by

Copyright:

Available Formats

Vortex Lattice Modelling of Winglets

on Wind Turbine Blades

Risø National Laboratory

3rd version. 18 oct. 2007

Abstract: ISSN 0106-2840

Information Service Department

2 winglet designs are threated in detail.

Vindmøllers energi-produktion kan øges ved brug af en winglet, uden at det

Formålet med nærværende projekt er at opnå en forståelse for, hvorledes en

En vortex-lattice kode og en hurtig design algoritme er blevet udviklet for en

2 winglet design er behandlet i detaljer.

Kongens Lyngby, August 2007

1.2 Winglet Description . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Overview of the Methods Applied . . . . . . . . . . . . . . . . . . 3

1.4 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1 Definition of Coordinate Systems . . . . . . . . . . . . . . . . . . 7

2.2 Definition of Blade and Winglet Geometry . . . . . . . . . . . . . 8

2.3 Dimensionless Variables . . . . . . . . . . . . . . . . . . . . . . . 9

2.4 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.5 The General Solution and Boundary Conditions . . . . . . . . . . 10

2.6 The Kutta Condition . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.7 Singularity Elements . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.8 Kutta-Joukowsky and Helmholz Theorems . . . . . . . . . . . . . 15

2.9 Wake Strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.10 Wake Shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.11 The Unsteady Bernoulli Equation . . . . . . . . . . . . . . . . . . 16

3 Numerical Lifting Line Model 19

3.2 Classic Analytical Lifting Line Results . . . . . . . . . . . . . . . 22

3.4 Effect of Winglet Height and Curve Radius . . . . . . . . . . . . 24

3.5 Optimum Circulation distribution on Wings . . . . . . . . . . . . 25

3.6 Generalization of Wing Results to Turbines . . . . . . . . . . . . 27

4 A Freewake, Vortex Lattice and Panel Method 31

4.1 Freewake Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.2 Evaluation of Forces . . . . . . . . . . . . . . . . . . . . . . . . . 37

4.3 Calculation of Wing Geometries for Vortex Lattice Simulations . 39

4.4 Grid Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

4.5 Reference Blade Data . . . . . . . . . . . . . . . . . . . . . . . . 41

5 Blade and Winglet Design 49

5.1 Defining Bound Circulation . . . . . . . . . . . . . . . . . . . . . 50

5.2 Calculating Velocities . . . . . . . . . . . . . . . . . . . . . . . . 52

5.3 Summary of Method . . . . . . . . . . . . . . . . . . . . . . . . . 56

5.5 General Design Results . . . . . . . . . . . . . . . . . . . . . . . . 60

6.1 Conservative Winglet Design . . . . . . . . . . . . . . . . . . . . 63

6.2 Optimum Winglet Design . . . . . . . . . . . . . . . . . . . . . . 68

A.1 Constant Strength Source Panel (3 Nodes) . . . . . . . . . . . . . 75

and the manufacturing process is not well established, a clear understanding of

1.2 Winglet Description

is created. On the main blade the direction of this unwanted downwash is

1.3 Overview of the Methods Applied

1.4 Outline of the Thesis

The contents of this thesis are presented in six parts.

Part 6 discusses the predicted performance of 2 selected final winglet designs.

Part 7 contains conclusions.

2.1 Definition of Coordinate Systems

Two coordinate systems are defined as seen in Figure 2.1.

Figure 2.1: The initial and the rotating coordinate system.

The initial coordinate system corresponds to an earth reference frame. The

coordinates are denoted with capital letters.

The rotational speed is Ω and the motion is therefore completely described by

2.2 Definition of Blade and Winglet Geometry