Aerodynamics of Vertical Axis Wind Turbines: Development of Simulation Tools and Experiments

UNIVERSITATIS ACTA UPSALIENSIS

UPPSALA 2015

Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1274

Aerodynamics of Vertical Axis Wind Turbines

Development of Simulation Tools and Experiments

EDUARD DYACHUK

ISSN 1651-6214 ISBN 978-91-554-9307-3 urn:nbn:se:uu:diva-260573

Dissertation presented at Uppsala University to be publicly examined in Polhemssalen, Ångströmlaboratoriet, Lägerhyddsvägen 1, Uppsala, Friday, 9 October 2015 at 09:00 for the degree of Doctor of Philosophy. The examination will be conducted in English. Faculty examiner: Mac Gaunaa (Technical University of Denmark, Department of Wind Energy).

Abstract

Dyachuk, E. 2015. Aerodynamics of Vertical Axis Wind Turbines. Development of Simulation Tools and Experiments. Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1274. 86 pp. Uppsala: Acta Universitatis Upsaliensis. ISBN 978-91-554-9307-3.

This thesis combines measurements with the development of simulation tools for vertical axis wind turbines (VAWT). Numerical models of aerodynamic blade forces are developed and validated against experiments. The studies were made on VAWTs which were operated at open sites. Significant progress within the modeling of aerodynamics of VAWTs has been achieved by the development of new simulation tools and by conducting experimental studies.

An existing dynamic stall model was investigated and further modified for the conditions of the VAWT operation. This model was coupled with a streamtube model and assessed against blade force measurements from a VAWT with curved blades, operated by Sandia National Laboratories. The comparison has shown that the accuracy of the streamtube model has been improved compared to its previous versions. The dynamic stall model was further modified by coupling it with a free vortex model. The new model has become less dependent on empirical constants and has shown an improved accuracy.

Unique blade force measurements on a 12 kW VAWT were conducted. The turbine was operated north of Uppsala. Load cells were used to measure the forces on the turbine. A comprehensive analysis of the measurement accuracy has been performed and the major error sources have been identified.

The measured aerodynamic normal force has been presented and analyzed for a wide range of operational conditions including dynamic stall, nominal operation and the region of high flow expansion. The improved vortex model has been validated against the data from the new measurements. The model agrees quite well with the experiments for the regions of nominal operation and high flow expansion. Although it does not reproduce all measurements in great detail, it is suggested that the presented vortex model can be used for preliminary estimations of blade forces due to its high computational speed and reasonable accuracy.

Keywords: wind power, vertical axis turbine, H-rotor, simulations, streamtube model, vortex model, dynamic stall, measurements, blade, force

Eduard Dyachuk, Department of Engineering Sciences, Electricity, Box 534, Uppsala University, SE-75121 Uppsala, Sweden.

© Eduard Dyachuk 2015 ISSN 1651-6214 ISBN 978-91-554-9307-3

urn:nbn:se:uu:diva-260573 (http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-260573)

To my mama & papa, Lidochka,

babushkam & dedushkam

List of papers

This thesis is based on the following papers, which are referred to in the text by their Roman numerals.

I Dyachuk, E., Goude A., and Bernhoff H., “Dynamic Stall Modeling for the Conditions of Vertical Axis Wind Turbines”, AIAA journal, vol. 52, iss. 1, pp. 72 – 81, 2014.

II Dyachuk, E., and Goude A., “Simulating Dynamic Stall Effects for Vertical Axis Wind Turbines Applying a Double Multiple Streamtube Model”, Energies, vol. 8, iss. 2, pp. 1353 – 1372, 2015.

III Dyachuk, E., Goude A., and Bernhoff H., “Simulating Pitching Blade With Free Vortex Model Coupled With Dynamic Stall Model for Conditions of Straight Bladed Vertical Axis Turbine”, Journal of Solar Energy Engineering, vol. 137, iss. 4, pp. 041008, 2015.

IV Rossander M., Dyachuk, E., Apelfröjd S., Trolin K., Goude A., Bernhoff H., and Eriksson S. “Evaluation of a Blade Force

Measurement System for a Vertical Axis Wind Turbine Using Load Cells”, Energies, vol. 8, iss. 6, pp. 5973 – 5996, 2015.

V Dyachuk, E., Rossander M., Goude A., and Bernhoff H.,

“Measurements of the Aerodynamic Normal Forces on a 12-kW Straight-Bladed Vertical Axis Wind Turbine”, Energies, vol. 8, iss. 8, pp. 8482 – 8496, 2015.

VI Dyachuk, E., and Goude A., “Numerical Validation of a Vortex Model Against Experimental Data on a Straight-Bladed Vertical Axis Wind Turbine”, Submitted to Energies, August 2015.

Reprints were made with permission from the publishers.

The following papers are not included in this thesis.

i Dyachuk, E., Goude A., Lalander, E., and Bernhoff H., “Influence of

Incoming Flow Direction on Spacing Between Vertical Axis Marine

Current Turbines Placed in a Row”, In “Proceedings of the 31st In-

ternational Conference on Offshore Mechanics and Arctic Engineering,

OMAE 2012”, Rio de Janeiro, Brazil, July 2012.

ii Dyachuk, E., Goude A., and Bernhoff H., “Simulating Pitching Blade With Free Vortex Model Coupled With Dynamic Stall Model for Condi- tions of Straight Bladed Vertical Axis Turbine”, In “Proceedings of the 33rd International Conference on Offshore Mechanics and Arctic Engi- neering, OMAE 2014”, San Francisco, USA, June 2014.

1

This article was further modified and published as Paper III.

Contents

1 Introduction

13

1.1 Background history of wind energy and different wind turbines 13 1.2 Main features of vertical axis wind turbines

14

1.3 Vertical axis wind turbines research at Uppsala University

16

1.4 Contribution of the thesis

16

1.5 Outline of the thesis

16

2 Theory

17

3 Development of simulation models

20

3.1 Streamtube model

20

3.2 Vortex model

24

3.3 Blade force modeling

27

3.3.1 Unsteady attached flow

27

3.3.2 Condition of dynamic stall

28

3.3.3 Separated flow

28

3.3.4 Application to the turbines

31

4 Measurements of aerodynamic forces

34

4.1 Experimental setup

34

4.2 Data treatment

34

4.3 Data treatment

37

4.4 Measurement accuracy

39

4.5 Evaluation of the measurement method

39

5 Results and discussion

42

5.1 Pitching blade

42

5.2 Turbines

51

5.2.1 12 kW VAWT

51

5.2.2 Sandia 17-m VAWT

56

6 Conclusions

65

7 Suggestions for future work

67

7.1 Simulation tools

67

7.2 Experimental work

67

8 Summary of papers

69

9 Acknowledgements

72

10 Sammanfattning

74

11 Аннотация к научной диссертации

77

12 Анотація до наукової дисертації

80

References

83

Nomenclature

A

m

Blade area

A

m

Cross-sectional area of a streamtube

A

m

Streamtube area at the mid-disk of a turbine A

m

Asymptotic cross-sectional area of a streamtube C

– Drag coefficient

C

– Drag coefficient at zero angle of attack C

– Lift coefficient

C

– Normal force coefficient

C

– Slope of normal force coefficient for attached flow C

– Normal force coefficient for trailing edge separation C

– Normal force coefficient during vortex convection C

– Tangential force coefficient

C

– Tangential force coefficient for static separation point C

– Tangential force coefficient for separated flow C

– Scaling constant for tangential force coefficient C

– Scaling factor for low pitch rate

C

– Scaling factor for low angles of attack C

– Static tangential force coefficient D

– Deficiency function for separation point

D

– Deficiency function for geometrical angle of attack E

– Constant for tangential force coefficient

Ω rad/s Rotational speed of a turbine F

N Measured force by load cell 0 F

N Measured force by load cell 1 F

N Measured force by load cell 2 F

N Measured force by load cell 3

F

N Drag force

F

N Centrifugal force

F

N Lift force

F

N Normal force

∆F

N Maximum measured error of normal force

F

N/m Normal force per unit span

∆F

N Maximum measured error of the normal force shape

F

N Tangential force

∆F

N Maximum measured error of tangential force F

N/m Tangential force per unit span

F

N Aerodynamic force on a blade in x-direction L

m Horizontal distance between the load cells L

m Distance from the load cells to the blade L

m Distance between load cells and center of mass N

– Number of turbine blades

P W Absorbed power by a turbine Q Nm Absorbed torque by a turbine

R m Turbine radius

S

, S

– Coefficients for separation point

T

C Air temperature

T

– Empirical time constant for dynamic separation point T

– Empirical time constant for delay in pressure response

V m/s Flow velocity

V

m/s Blade velocity

V

m/s Velocity at the turbine disk V

m/s Velocity at the turbine mid-disk V

m/s Relative flow velocity

V

m/s Velocity due to vortices V

m/s Asymptotic wind velocity

X ,Y, Z – Deficiency functions for unsteady attached flow c m Chord length of a blade

f – Static flow separation point

f

– First-order delayed flow separation point f

– Dynamic flow separation point

g m/s

Gravitational acceleration

h %RH Relative humidity

k – Reduced frequency of a pitching blade

∆l m Length of a blade segment in a streamtube m kg Mass of the blade with support arms p

Pa Pressure directly in front of turbine disk p

Pa Pressure directly behind turbine disk p

Pa Asymptotic pressure

r m Arbitrary position

r

– Critical reduced pitch rate

r

m Vortex position

r

– Lower scaling limit of reduced pitch rate r

– Reduced pitch rate

r

– Upper scaling limit of reduced pitch rate

∆s – Non-dimensional time-step

t s Time

∆t s Time-step

u – Velocity induction factor ˆ

x – unit vector in x-direction

x

– Normalized blade attachment point

∆z m Height of a blade segment

Γ m

/s Circulation of two-dimensional vortex Γ

m

/s Total two-dimensional circulation of a blade

α – Angle of attack

α

– Delayed angle of attack due to pressure delay α ˙ rad/s Blade pitch rate

α

– Angle of attack for a breakpoint of flow separation α

– Effective angle of attack

α

– Critical angle of attack α

– Critical stall onset angle

α

– Lower scaling limit of angle of attack α

– Static stall onset angle

α

– Upper scaling limit of angle of attack β – Angle between a blade and vertical axis

γ – Angle between incoming flow and a streamtube

δ – Blade pitch angle

ε m Cutoff radius of Gaussian vortex kernel

η – Efficiency factor for tangential force coefficient

θ – Blade azimuth angle

θ ˆ – Unit vector in the tangential direction

∆θ – Angular width of a streamtube

λ – Tip speed ratio

ν m

/s Kinematic viscosity ρ kg/m

Air density

ϕ – Angle of relative flow velocity φ m

/s Velocity potential

ω 1/s Vorticity

Abbreviations

BEM Blade element momentum model CFD Computational fluid dynamics

DS Dynamic stall model coupled with streamtube model FEM Finite element method

FVM Finite volume method HAWT Horizontal axis wind turbine LES Large eddy simulation

NACA National Advisory Committee for Aeronautics NASA National Aeronautics and Space Administration RANS Reynolds-averaged Navier Stokes equations TSR Tip speed ratio

VAWT Vertical axis wind turbine

VDS Dynamic stall model coupled with vortex model

1. Introduction

This thesis comprises studies on the aerodynamics of vertical axis wind tur- bines. The focus is on unsteady aerodynamic forces which are usually consid- ered as the major challenge for the development of vertical axis wind turbines.

Within the presented work, simulation tools have been developed and experi- mental work has been conducted.

1.1 Background history of wind energy and different wind turbines

The utilization of wind energy dates back to ancient times when people started to use sails to propel boats and ships. Along with transport, the use of wind energy was adopted in agriculture to grind grain and pump water. The earliest documented wind mills were used by the Persians more than two thousand years ago [1]. Those were drag-driven vertical axis mills with sails made of reeds or wood. Later, knowledge about wind mills was brought from Persia and the Middle East to Europe most likely by the Crusades [2]. In contrast to earlier vertical axis design, European wind mills had a horizontal axis of rotation. The first wind turbine generating electricity was built in 1887 by Scottish scientist and engineer Professor James Blyth [3, 4].

The development of wind turbines continued in the 1900s in Europe and the USA. A Finnish inventor Sigurd Johannes Savonius in the early 1920s devel- oped a vertical axis wind turbine (VAWT) which was driven mainly by drag force. Two US patents on this turbine were published, one in 1929 [5] and another one in 1930 [6]. The invention was mentioned as the Savonius rotor in the patents. It was suggested that the turbine could be driven either by wind or flowing water. Another type of VAWT was invented by a French aeronau- tical engineer Georges Jean Marie Darrieus and patented in 1931 [7]. The patent emphasizes simplicity of the turbine and covers turbines with curved and straight blades. It is also suggested in the patent that the Darrieus turbine can be utilized in tidal currents or rivers with low fall. The Darrieus turbine is a lift-based machine (torque is mainly generated by lift force) which is more efficient than the Savonius turbine and requires less material to build it (see figure 1.1).

Due to the cheap electricity from fossil fuels before the 1973 oil crisis, only

small scale wind turbines were used in remote locations which were not con-

nected to an electrical grid. However, the oil crisis has led to the growing

interest in electricity production from renewable energy sources on a larger scale. Wind energy in the USA has received financial support from the gov- ernment. The research was focused on the development of new hardware, resource analysis and cost reduction techniques. Along with the development of horizontal axis wind turbines (HAWT) by the National Aeronautics and Space Administration (NASA), Sandia National Laboratories from the 1970s until the early 1990s focused the research on Darrieus turbines with curved blades. Several companies have commercialized VAWTs with curved blades.

One of these companies was FloWind, which deployed over 500 turbines in California.

In the UK, development of the VAWT technology was mostly applied to turbines with straight blades, (the H-rotor design, see figure 1.1), originally proposed and developed in the late 1960s by British engineer and scientist Peter Musgrove [8]. A mechanism of feathering the blades at strong winds was later included in order to reduce the aerodynamic torque. The turbine with feathering blades was later referred to as the Musgrove rotor. It was commercialized by several companies; among them was VAWT Ltd, which in the 1980s built several VAWTs with rated power up to 500 kW.

However, the knowledge about blade fatigue at that time was limited, and both the American and the British turbines suffered from the blade failures.

In particular, the blades of the Musgrove rotor turbines in UK broke due to manufacture error. The curved blades of the FloWind turbines made of ex- truded aluminium had poor fatigue properties, which resulted in several blade failures. Those blade failures led to a common perception that VAWTs are more prone to fatigue than HAWTs. The British VAWT Ltd was closed in the early 1990s, and the American FloWind was closed in the mid-1990s. Since then, the interest in VAWTs has decreased within the wind energy commu- nity. A more detailed history of the development of VAWTs can be found in Refs. [9–12].

1.2 Main features of vertical axis wind turbines

The absolute majority of the currently deployed wind turbines are HAWTs.

However, the VAWT design has several principal advantages over the HAWT.

VAWTs are omni-directional, i.e. they work with winds from any direction

and the yawing system is excluded. This is a considerable advantage, since

a large portion of the failures in HAWT occur within the yawing mecha-

nism [13–15]. Another advantage is that the generator of VAWTs can be

placed at the ground level, which simplifies installation and maintenance. Ad-

ditionally, this minimizes concerns about the size and the weight of the gener-

ator, which is favourable for the installation of heavy direct driven generators

with permanent magnets [16]. A concern for VAWTs is the cyclic blade forces

resulting in varying torque, which is inherent during operation. Although this

Savonius rotor Darrieus rotor H-rotor

Figure 1.1. Different types of vertical axis wind turbines

problem was handled by Sandia National Laboratories and FloWind by adding compliance to the shafts of their VAWTs [9], the varying torque is still consid- ered the major concern of VAWTs.

The accessibility of generators together with the low mass center of VAWT is specifically advantageous for offshore applications including floating tur- bines. Several studies on VAWTs for offshore use have been published [17–

19]. Sandia National Laboratories has summarized properties of the Darrieus turbine from the perspective of the offshore deployment of VAWTs [9].

This thesis focuses on VAWTs with straight-bladed H-rotors, although ex-

perimental data from a Sandia turbine with curved blades have been used (see

chapter 5). The straight blades of the H-rotor are easier and cheater to man-

ufacture than the curved blades of the Darrieus rotor. The bending moments

on the shaft of the H-rotor are much smaller than on the shaft of the Darrieus

turbine with curved blades, since the upper bearing can be placed close to

the turbine hub, reducing the bending moments in the axis. Instead, it is the

tower in the H-rotor design that takes most of the bending loads from the tur-

bine bearings. Hence, the shaft of the H-rotor is easier to manufacture than the

shaft of the curve-bladed turbines with a small tower. Additionally, the H-rotor

design offers larger cross-sectional area due to the constant radius of the tur-

bine. However, the blades of the H-rotor need support arms, which introduce

additional losses which affect the turbine performance. Another drawback of

the H-rotor, compared to the design with curved blades, is the higher bending

moments in the blades due to centrifugal forces.

1.3 Vertical axis wind turbines research at Uppsala University

The wind power research at the Division of Electricity at Uppsala University has been conducted since 2002. Three H-rotor VAWTs have been built north of Uppsala by the division: a small one with rated power of 1.5 kW, one turbine for telecom applications rated at 10 kW [20] and a 12 kW turbine, which has been used for the most of experiments, see Papers IV, V and Refs. [21–

23]. A large VAWT with rated power of 200 kW has been constructed by the spinoff company Vertical Wind AB in Falkenberg and some published data are available for that turbine [24–26]. The research related to these turbines has resulted in five doctoral theses [27–31]. Within this research, simulation tools for aerodynamics of VAWTs were developed by Anders Goude [30] and Paul Deglaire [28].

1.4 Contribution of the thesis

The focus of this thesis is on the aerodynamic blade forces of the VAWTs.

The work comprises both development of simulation tools and force measure- ments. A model of unsteady forces on a pitching blade for the conditions of VAWTs has been developed and assessed against existing experimental data.

The model is commonly known as the dynamic stall model and it is published in Paper I. The dynamic stall model has been included into a simulation model of the turbine and compared with existing measurements on a VAWT with curved blades. This study is published in Paper II. A further development of the dynamic stall model is published in Paper III, where the model is coupled with a free vortex model for a pitching blade.

Novel measurements of the aerodynamic blade forces on the straight-bladed 12 kW VAWT have been conducted and the evaluation of experimental method is published in Paper IV. The measurements of the aerodynamic normal forces are published in Paper V. A modified vortex model coupled with the dynamic stall model is validated against the new measurements on the 12 kW VAWT.

The validation is documented in Paper VI.

1.5 Outline of the thesis

The introduction of this thesis is followed by the background theory of the

aerodynamics of VAWTs. The overview of simulation models is presented in

chapter 3. This is followed by chapter 4, where the experimental method of

the force evaluation is described. The assessment of the simulations against

the experimental results is presented in chapter 5 along with the discussions

regarding applicability of the models. At the end of the thesis, the conclusions

and suggestions for future work are presented.

2. Theory

The power of wind turbine can be expressed in terms of torque as

P = QΩ, (2.1)

where Q is the torque and Ω is the turbine rotational speed. The focus of this work is on lift-based turbines. The torque of these turbines is generated mainly by lift force, and drag force contributes to losses. The structural loads on the turbine blades are often given by the normal force F

which is the resultant of the aerodynamic force in radial direction, figure 2.1. The tangential force F

is usually used to express the turbine torque during one revolution

Q = N

hF

iR, (2.2)

where N

is the number of blades, R is the turbine radius and hF

i is the aver- age tangential force during one revolution. If the turbine blades are supported by struts (e.g. H-rotor) the tangential force comprises the contribution of the blades and support arms. Both the turbine geometry and operational condi- tions influence the tangential force. To illustrate this, assume that the blade is located at the azimuth angle θ (see figure 2.1). The velocity of the blade is

~ V

= ΩR ˆ θ , (2.3)

where ˆ θ is the unit vector in the tangential direction which is positive in the counter-clockwise direction. The vector of relative wind flow velocity ~ V

at the blade consists of the vector of incident wind flow at the turbine disk ~ V

and the velocity vector at the blade due to its own rotation −~ V

(negative sign indicates that the direction of the flow is opposite to the direction of the blade)

~ V

= ~ V

−~V

. (2.4) Due to the extracted energy from the flow, the magnitude of the wind velocity at the turbine disk V

is generally lower than the asymptotic velocity V

. In Cartesian coordinates, when the asymptotic velocity is aligned with x-axis and no flow expansion is assumed (i.e. ~ V

= V

x ˆ and ~ V

= V

x), the magnitude of ˆ the relative flow velocity is

|~V

| = V

s

 ΩR V

+ sin θ



+ (cos θ )

, (2.5)

V

φ α δ

F

F

F

θ

θ = 0

F

→ V

→

V

→

V

→

Figure 2.1. Illustration of the velocity vectors and forces acting at the blade of a VAWT. Counter-clockwise direction is defined as positive for the angles. Hence, the angles ϕ and α have negative direction for the position of the blade shown in this figure. The normal force F

is positive when pointing outwards, i.e. F

is negative in the figure.

and the angle of the relative wind is

ϕ = arctan cos θ

ΩR V_d

+ sin θ

!

. (2.6)

The estimation of the relative wind angle according to equation (2.6) is only valid assuming the blade as a point. The blade performs rotational mo- tion, indicating that there are flow curvature effects which change the effective angle ϕ. To account for the curvature effects, the relative wind angle is fur- ther modified. The assumption of the potential flow around a flat plate and the Joukowski transformation together with the Kutta condition are used as by Goude [30]. The final expression of the relative wind angle including the flow curvature effects is

ϕ = arctan cos θ

ΩR V_d

+ sin θ

!

− Ωx

c V

− Ωc

4V

, (2.7)

where x

is the normalized blade attachment point which varies from −0.5 to 0.5 with x

= 0 corresponding to the blade attachment at the middle of the chord length.

The angle of attack is the sum of the relative wind angle and the blade pitch angle

α = ϕ + δ . (2.8)

One important operational parameter for wind turbines is a tip speed ratio (TSR) λ , which is defined as the ratio between the speed of a blade tip and the asymptotic flow velocity

λ = ΩR

V

. (2.9)

The term

d

in equation (2.7) should not be interpreted as the TSR, although

ΩR

V_d

is proportional to λ since V

∝ V

. From equation (2.7) it follows that the relative wind angle ϕ increases with decreased TSR, and when the pitch angle δ is fixed

α ∝ 1

λ . (2.10)

To relate this to the turbine torque, the expression for the tangential force F

can be used

F

= F

sin ϕ − F

cos ϕ. (2.11) The lift force F

is orthogonal to the vector ~ V

and the drag force F

is parallel to ~ V

. The lift and drag forces can be estimated as

F

= 1

2 ρ A

V

C

, (2.12)

F

= 1

2 ρ A

V

C

, (2.13)

where ρ is the flow density, A

is the blade area, C

and C

are the lift and drag coefficients, which mainly are the functions of the angle of attack, aspect ratio, the Reynolds number and the airfoil profile. For airfoils, C

increases with increased α and C

is relatively constant. This holds until the point of stall (associated with the flow separation), when C

starts to decrease and C

increases. Thus, the maximum tangential force is obtained when the angles

of attack are close to the stall region. Consequently, the turbine should be de-

signed to operate at TSRs corresponding to the nearly stall limit in order to

obtain the highest power for the given wind conditions. Please note that dur-

ing the operation of VAWT, the angle of attack and the relative wind velocity

change constantly, which causes the phenomenon dynamic stall. During dy-

namic stall, the coefficients C

and C

will have a delay compared to the static

flow case. This should be taken into account when operating VAWT close to

the stall limit. The details regarding the dynamic stall event are presented in

section 3.3.

3. Development of simulation models

The simulation models for aerodynamics of VAWTs can be divided into three groups. The first group includes the models which solve the Navier-Stokes equations on a grid of the entire simulated volume. Due to its very high com- putational complexity, the Navier-Stokes equations are usually combined with turbulence models such as the Reynolds-averaged Navier-Stokes equations (RANS) [32–35] and large eddy simulations (LES) [36–38]. The most com- mon methods for this group of models are the finite element method (FEM) and the finite volume method (FVM) which solve the flow for a confined re- gion. There are commonly available FEM and FVM models, however due to their high computational cost, these models are not used in the presented work.

The second method for simulating aerodynamics of VAWTs is the use of the vorticity equation, which is based on the Navier-Stokes equations. Many simplifications can be applied within this method, and when a model for blade forces is used, the computational time can be reduced substantially. A vortex model has been used within the work.

The third type of the models for the VAWT simulations is based on the mo- mentum conservation principle. These models are usually referred to as blade element momentum models (BEM), when simulating HAWT, and streamtube models, for simulating VAWT. The models assume steady flow, and the full description of the flow through entire turbine is not obtained. However due to their simplicity and high computational speed, the momentum models are used to quickly evaluate the performance of VAWTs.

This chapter presents the summary of the simulation tool development.

More details concerning the development of each model is presented in Pa- pers I, II, III and VI.

3.1 Streamtube model

A double multiple streamtube model is used in this work, and the implementa- tion of Paraschivoiu [39] is applied. The difference between the implemented model and the version by Paraschivoiu is that the positive flow direction is towards right, see figure 3.1. The flow expansion by Read and Sharpe [40]

is further added into the model. The present version of the model has been

tested against the experimental data on a curve-bladed Darrieus turbine by

Sandia National Laboratories, see Paper II.

A

Δθ

V

V

V

θ

Ω

p

p

p

p

Figure 3.1. Illustration of the double multiple streamtube model. The horizontal lines represent streamtubes, and the vertical line in the middle of the turbine represents the transition from the upwind to the downwind side. The height of a streamtube ∆z

is included into the streamtube area A

. Note, that the flow expansion is not included in this figure.

In the double multiple streamtube model, the turbine is divided into two ac- tuator disks, one upwind and one downwind, figure 3.1. It is assumed that the velocity across the actuator disks remains constant and the pressure in the mid- dle of the turbine equals to the asymptotic pressure. With these assumptions and by using the Bernoulli’s equation at the upwind disk

p

+ 1

2 ρV

= p

+ 1

2 ρV

, (3.1) p

+ 1

2 ρV

= p

+ 1

2 ρV

, (3.2)

combined with the momentum conservation for the control volume of a single streamtube

ρ A

V

− ρA

V

= A

p

− p

, (3.3) and with continuity

A

V

= A

V

= A

V

, (3.4) it can be shown that the the velocity at the upwind disk is the average of the asymptotic velocity and the velocity at the middle of the turbine

V

= V

+V

2 . (3.5)

β

Δz

Δl

V

Figure 3.2. Illustration of the angle β for a curve-bladed VAWT.

The asymptotic flow velocity V

, the velocity at the turbine disk V

and the velocity at the middle of the turbine V

are illustrated in figure 3.1 together with the asymptotic pressure p

, the pressure directly in front of the disk p

and directly behind the disk p

. The cross sectional areas of a streamtube in equations (3.3) to (3.5) are denoted as A

, A

and A

.

To calculate the forces on the turbine, the lift and drag coefficients together with the flow velocity have to be estimated. The expressions for the flow velocity and the angle of attack including the curvature effects are given in equations (2.5), (2.7) and (2.8). For the case with the curve-bladed VAWTs (see figure 3.2), V

and α

become

V

= V

s

 ΩR

V

+ sin θ



+ (cos θ

cos β

)

, (3.6)

α

= ϕ

+ δ = arctan





cos θ

cos β

V_{d, j}

+ sin θ



 −

 Ωx

c

V

+ Ωc 4V



cos β

+ δ , (3.7) where β

is the angle between the blade in a streamtube and vertical axis.

The lift and drag coefficients are obtained from the blade force model, which is described in section 3.3. The lift and drag forces are estimated as

F

= 1

2 ρ c∆l

V

C

, (3.8) F

= 1

2 ρ c∆l

V

C

, (3.9) where ∆l

is the blade length in a streamtube

∆l

= ∆z

cos β

(3.10)

and ∆z

is the height of a streamtube, see figure 3.2. As mentioned in chapter 2, it is the normal force F

(giving the structural loads on the blades) and the tangential force F

(corresponding to the turbine torque) which are of interest when designing the turbine. The normal and tangential force coefficients are defined as

C

= C

cos ϕ

+C

sin ϕ

, (3.11) C

= C

sin ϕ

−C

cos ϕ

, (3.12) which can be used to estimate the normal and tangential forces

F

= 1

2 ρ c∆l

V

C

, (3.13) F

= 1

2 ρ c∆l

V

C

. (3.14) Since the turbine has N

blades, and each blade passes a streamtube during a fraction of time, the average force on the turbine blades in the x-direction is

F

 = N

∆θ

2π F

, (3.15)

where F

is

F

= F

cos θ

cos β

− F

sin θ

. (3.16) Combining equations (3.10), (3.13), (3.14) and (3.16) with equation (3.15), the force F

 becomes

F

 = N

ρ c∆z

∆θV

4π



C

cos θ

−C

sin θ

cos β



. (3.17)

The force F

 is equal to the pressure difference about the streamtube area A

, equation (3.3). Thus, combining an expression for A

A

= R∆z

∆θ | cos θ

|, (3.18) with equations (3.3) and (3.5), it is possible to find the force F

 from mo- mentum conservation

F

 = ρR∆z

∆θ | cos θ

|V

(V

−V

) . (3.19) Combining equations (3.17) and (3.19) with the velocity induction factor

u

= V

V

, (3.20)

and equation (3.5), the following expression is obtained 1

u

− 1 = N

cV

8πR| cos θ

|V



C

cos θ

−C

sin θ

cos β



. (3.21)

Equation (3.21) is solved iteratively to obtain the induction factor u

for each streamtube. The velocity at the middle of the turbine is found based on equa- tions (3.5) and (3.20)

V

= 2V

(1 − u

) . (3.22) To solve the downstream disk, the velocity V

is used as an input, and the blade forces are calculated the same way as described above.

After the estimation of velocities at the both sides of the turbine, the flow expansion model as by Read and Sharpe [40] is included. The model recalcu- lates the size of streamtubes by applying the continuity equation (3.4) at the upwind and downwind sides of the turbine. Please note, that the expansion is modeled in the horizontal plane only. The cross-section of a streamtube is shown in figure 3.3. It is assumed that the flow expands symmetrically about the horizontal line which connects θ = 90

and θ = 270

. An iterative process is required to find the size of each streamtube for known V

at the upwind and the downwind side. New azimuthal angle θ

is obtained as

θ

= θ

− γ

, (3.23)

and the angular streamtube width ∆θ is recalculated. The angle γ

represents the deviation of streamlines from the direction of the asymptotic flow due to the difference between the streamtube area at the upwind and the downwind sides. After the size of each streamtube is estimated, the flow velocities and the blade forces are calculated again.

3.2 Vortex model

The motion of the fluid can be described with the continuity equation and the Navier-Stokes equations. The continuity equation

∂ ρ

∂ t + ∇ ·  ρ~ V



= 0, (3.24)

represents the conservation of mass. The flow is assumed to be incompressible due to the low flow velocities at VAWTs, and equation (3.24) reduces to

∇ ·~ V = 0. (3.25)

The Navier-Stokes equations are derived from the Newton’s second law, and in the case of incompressible flow, the equations are

∂~ V

∂ t + 

~ V · ∇ 

~ V = − ∇ p

ρ +~g + ν∇

~ V , (3.26) where ~g is the gravitational acceleration and ν is the kinematic viscosity.

Equation (3.26) represents the system of non-linear partial differential equa-

tions, which can be solved analytically only in rare cases. One method to

V

γ A

V

A

V

θ

θ

θ = 90

θ = 270

Figure 3.3. Illustration of the flow expansion model. The doted line represents the center of a streamtube. Indices up and dw denote the upwind and the downwind side accordingly. Please, note that the counter-clockwise direction is defined positive for γ , i.e. γ is positive for the streamline in this figure.

solve the Navier-Stokes equations numerically is to use the finite element or finite volume models, where the equations are solved on a grid. An alternative method is to use the vorticity equation, which is obtained by taking the curl of equation (3.26)

∂ ~ ω

∂ t +

 ~ V · ∇ 

~ ω = (~ ω · ∇)~ V + ν∇

~ ω , (3.27) where the vorticity is

~ ω = ∇ ×~ V . (3.28)

In the two dimensional (2D) case, the vorticity is aligned with the z-axis and the flow occurs in the x − y plane, i.e. ~ ω = ω ˆ z and ~ V = V

x ˆ + V

y. Thus in the ˆ two dimensional case, equation (3.27) reduces to

∂ ~ ω

∂ t + 

~ V · ∇ 

~ ω = ν ∇

~ ω , (3.29) Since experimental data and a blade force model are used to obtain the lift and the drag forces (i.e. the boundary layer does not have to be resolved), viscosity in equation (3.29) can be neglected. This significantly increases the computational speed of the vortex method.

The flow velocity is obtained from

~ V = ∇φ +~ V

, (3.30)

where ∇φ is the potential flow solution and ~ V

is the contribution from re- leased vortices. Under the assumption that VAWTs are not confined with boundaries and that the blades are approximated as single points, the poten- tial flow solution equals to the asymptotic flow velocity, i.e. ∇φ = V

. The contribution from the vortices is calculated from the Biot-Savart law by using complex variables instead of vectors, which in the 2D case becomes

V

(r) =

Nv

∑

k=1

iΓ

2π 1

(r − r

) 1 − e

|r−rk|²

ε 2

!

, (3.31)

where r is the arbitrary position, r

is the position of vortex k, Γ

is the circu- lation of vortex k (r and r

denote the complex conjugate of r and r

), and ε is the smoothing parameter for the vortex when applying a Gaussian smoother.

The vortices are allowed to drift with the flow velocity when using the La- grangian description. By neglecting the viscosity outside the boundary layers of the blades (due to the high Reynolds numbers), the vortices are propagated according to

dr

dt = V (r

) . (3.32)

The velocities V (r

) at each vortex position are efficiently evaluated using the fast multipole method as described in Ref. [41].

The blade force model (described further in section 3.3) requires the flow velocity and the angle of attack to calculate the forces at the blade. The flow curvature effects have to be taken into account when calculating the flow ve- locity and the angles of attack. The flow curvature can be handled by using linear panels with linear distribution of vorticity to model the blade geometry according to Refs. [30, 42]. By applying the Kutta condition at the trailing edge, the circulation around the blade can be determined. The angle of attack corresponding to the blade circulation can then be calculated as

α

= k

 arcsin

 Γ

π cV



− α



, (3.33)

where k and α

are two constants which are determined by curve-fitting from a set of static attached flow simulations at known angles of attack. The ref- erence flow velocity V in equation (3.33) is calculated with equations (3.30) and (3.31) at the quarter chord position for the blade, under the condition, that the panels for blades are not present during this calculation. This whole pro- cedure is described in detail and denoted as “explicit method” in Paper III.

This method adds the released vortex into the panel equations and uses the assumption that the total circulation of the released vortex and the bound cir- culation is equal to the bound circulation form previous time-step, i.e. the total circulation is conserved.

The blade force model estimates the lift force coefficient C

for the known

angle of attack α

and flow velocity V . This lift force coefficient is used to

calculate a reduced circulation around the blade by using the Kutta-Joukowski lift formula

Γ = 1

2 C

cV. (3.34)

Due to the conservation of circulation, a vortex has to be released at each time- step, with a strength corresponding to the change of circulation between the time-steps.

To ensure the convergence of the model results, the simulations are per- formed for 100 revolutions. This value is chosen from the convergence stud- ies in Ref. [43]. Each revolution has 120 time-steps, i.e. one time-step corre- sponds to ∆θ = 3

. The results from the last revolution are used in the model assessment.

3.3 Blade force modeling

The blade forces in both the streamtube and the vortex models are estimated with lift and drag coefficients. The calculation of the lift and drag coefficients have to account for the continuous change of the angle of attack and the rel- ative wind velocity, which are natural during the operation of VAWTs, see equations (2.5), (2.7) and (2.8). For high wind speeds, where the rotational velocity is limited, the TSR decreases, which causes an increase in the ampli- tude of the angles of attack. Consequently, the blades fall into a stall, causing a drop in the lift force. Since the flow around the blades is not steady, dynamic stall is present, which is associated with the delay in the lift and drag force coefficients compared to the static values. Moreover, the amplitude of the os- cillations of the blade forces changes with dynamic stall, which is related to fatigue problems.

This section presents the dynamic stall model which is based on the model by Leishman-Beddoes [44, 45]. The model is semi-empirical and requires the local Reynolds number, the rotational speed, the time-step, the angle of attack, the flow velocity and the airfoil profile. Additionally, the model uses experi- mental data on static lift and drag coefficients [46]. This model is described in detail in Paper I and the modifications for the conditions of VAWTs are presented in Papers II and IV. The model consists of three parts: unsteady attached flow, stall onset and separated flow.

3.3.1 Unsteady attached flow

The time varying bound vortex during the unsteady attached flow is repre- sented by an effective angle of attack

α

= α

− X

−Y

− Z

, (3.35)

where α is the geometrical angle of attack and X , Y and Z are the deficiency functions, which are empirically derived based on the flow velocity and the pitch rate, and they can be found in Paper I. Indices n and n − 1 in equa- tion (3.35) represent the current and previous time-steps. The delayed angle of attack due to the lag in pressure response is calculated as

α

= α

− D

, (3.36)

where D

is the deficiency function D

= D

exp



− ∆s T



+ (α

− α

) exp



− ∆s 2T



, (3.37)

with an empirically derived time constant T

, which is found in table 3.1. The non-dimensional time-step ∆s in equation (3.37) is calculated as

∆s = 2V ∆t

c . (3.38)

Here and in equation (3.40), the velocity V stands for the relative flow velocity V

from equation (3.6) for the streamtube model and for the relative flow velocity |~ V −~V

| from equations (2.4) and (3.30) for the vortex model.

3.3.2 Condition of dynamic stall

Due to flow reversal within the boundary layer, a leading edge vortex forms at the airfoil surface. The critical angle of attack α

is used to define the condition at which the dynamic stall may begin

α

= (

α

|r

| ≥ r

,

α

+ (α

− α

)

0

|r

| < r

, (3.39) where the reduced pitch rate r

is

r

= α ˙

c

2V . (3.40)

Here, ˙ α is the pitch rate, r

is the reduced pitch rate which delimits the quasi- steady stall and the dynamic stall; r

= 0.01 for symmetrical NACA-airfoils, see Paper I. The static stall onset angle α

and the critical stall onset angle α

are listed in table 3.1. The following dynamic stall condition is used

  α

 > α

→ stall. (3.41)

3.3.3 Separated flow

The effects of separated flow are divided into two groups: trailing edge sep-

aration and leading edge vortex convection. The trailing edge separation is

associated with the time delay in the movement of the boundary separation point, and it is obtained via Kirchhoff’s flow approximation

f

=



 



 



1 − 0.4 exp  |

|

S₁



|α

| < α

, 0.02 + 0.58 exp



α₁−

|

|

S₂



|α

| ≥ α

,

(3.42)

where f

is the delayed separation point and α

, S

and S

are constants based on the airfoil profile and the local Reynolds number, found in Paper I. The boundary layer around the blade itself is time dependent, which is superim- posed on the pressure response delay. This additional delay is represented by the dynamic separation point

f

= f

− D

. (3.43) Here, the deficiency function D

is

D

= D

exp



− ∆s T



+ f

− f

exp



− ∆s 2T



, (3.44)

where T

is an empirically derived time constant, table 3.1. The normal force coefficient for the unsteady conditions before the dynamic stall onset is ob- tained as

C

n

= C

α

1 + p f

2

!

. (3.45)

After the stall condition is met, the leading edge vortex convects over the surface of the airfoil towards the trailing edge. During this process, a signifi- cant increase in the normal force is present

C

= B

f

− f
V

, (3.46) where C

is the normal force coefficient during the vortex convection (so- called “vortex lift”), V

and B

are parameters based on the local Reynolds number and the airfoil profile and are found in Paper I. The normal force decreases rapidly when the vortex passes the trailing edge. The total normal coefficient is estimated as

C

= C

+C

. (3.47) Figure 3.4 shows an example of the normal force coefficient of a pitching blade, simulated by the dynamic stall model, and the features described above are present.

The calculation of the tangential force coefficient is based on the Kirch- hoff’s flow relation using the dynamic separation point with the modification by Sheng et al. [47]

C

= ηC

α

p

f

− E



, (3.48)

0 5 10 15 20 25 0.2

0.4 0.6 0.8 1 1.2 1.4 1.6

Angle of attack α [deg]

Normal force coefficient C

Formation of

leading edge vortex, dynamic stall onset

Vortex convection

Vortex passes trailing edge, full stall Start of the flow

reattachment

Figure 3.4. Normal force coefficient for the pitching NACA0021 airfoil of a periodical motion α = 12 + 10 sin ( ˙ αt), ˙ α = 7.20 rad/s, V

= 33 m/s, c = 0.55 m.

where η and E

are the empirical constants listed in table 3.1.

For low angles of attack, the values of C

approach the static values accord- ing to equations (3.45) and (3.47). However, the tangential force coefficient C

does not approach the static value at low angles of attack when calcu- lated with equation (3.48). Therefore, a modification is applied to solve this limitation. The coefficient C

is calculated for the static separation point f

C

= ηC

α

 p f − E



. (3.49)

Depending on the absolute value of the geometrical angle of attack α and the reduced pitch rate r

, the tangential force coefficient C

is estimated as

C

= C

+ C

−C

C

, (3.50) where C

is the tangential force coefficient from the static lift and drag data [46]. The scaling constant C

includes the scaling factor for the angle of attack C

C

=



 

 

0 |α| > α

, 1 |α| < α

,

αu−|α|

αu−αl

α

≤ |α| ≤ α

,

(3.51)

Table 3.1. Empirical constants for the dynamic stall model model

Airfoil T

α

[deg] α

[deg] B

η E

NACA0012 3.90 14.95 18.73 0.75 1 0.25

NACA0015 5.78 14.67 17.81 0.50 1 0.25

NACA0018 6.22 14.68 17.46 0.50 1 0.20

NACA0021 6.30 14.33 17.91 0.50 0.975 0.15

NACA0025 6.95 13.59 17.22 0.50 0.90 0.18

and the scaling factor for the pitch rate C

C

=



 

 

0 |r

| > r

, 1 |r

| < r

,

r_u−|r_n|

ru−r_l

r

≤ |r

| ≤ r

.

(3.52)

The scaling constant C

in equation (3.49) is defined as the product of the scaling factors

C

= C

C

. (3.53) The parameters α

and α

in equation (3.51) together with r

and r

in equa- tion (3.52) represent the upper and the lower scaling limits of the angle of attack and the reduced pitch rate respectively. These limits are used to ensure the smooth transition of the tangential force coefficient to its static values at low angles of attack and low pitch rate. This method is similar to the one pub- lished in Paper IV with the exception that the pitch rate is taken into account in the present implementation. The scaling limits of the angle of attack are α

= α

and α

= 0.5α

and the limits of the reduced pitch rate are r

= 2r

and r

= r

. Please note that these values of the scaling limits have been cal- ibrated for the simulation cases in the current work. Further investigation of the scaling limits is recommended if data for other turbines become available.

3.3.4 Application to the turbines

Due to the circulatory motion of the blades during the turbine operation, the blade vortices may be released faster at the downwind disk. Figure 3.5 shows the vortex shedding structure of a Darrieus turbine at low TSR. This figure is based on the velocity measurements of the straight-bladed Darrieus turbine operating at a TSR of 2.14 in a water channel, obtained by Brochier et al. [48].

It is shown that the dynamic stall vortices are released and swept away when the blade passes Quadrant III; thus, the flow is fully separated. This can be both due to the circulatory motion of the blade and due to the highly turbulent flow, as noted in Ref. [48].

An additional modification to the dynamic stall model is applied to account

for the vortex shedding. This is done by modeling fast release of the trail-

ing and leading edge vortices: the delays in the angle of attack, the leading

0º

90º

180º

4

5 4

5

1 1

2

3

3 2

a

a a' a

a'

a a

b

b b b

b

IV c I

III II

V

Figure 3.5. Flow visualization in the dynamic stall condition at λ = 2.14, taken from Ref. [48]. a, a

, b and c denote vortices.

edge separation point and the vortex lift are set to zero when the blade passes Quadrant III:

Quadrant III → α

= α

, C

= 0 (3.54) When the normal and tangential force coefficients are found, the lift and drag coefficients are estimated as

C

= C

cos ϕ +C

sin ϕ, (3.55)

C

= C

sin ϕ −C

cos ϕ +C

1 −C

, (3.56) where C

is the drag coefficient at the zero angle of attack. Note that C

is already taken into account in C

for low angles of attack and low pitch rate (C

= 1) through C

according to equation (3.50). To assure a smooth transition at increased angles of attack and pitch rate, C

is linearly interpolated to avoid applying the value twice.

The integration of the dynamic stall model is different between the stream- tube and the vortex model. The angles of attack and the flow velocities in the streamtube model are calculated for the whole revolution forming an input matrix to the dynamic stall model. Then, the dynamic stall model runs two loops, which is sufficient to approach the convergence, see Paper I. The lift and drag coefficients from the dynamic stall model are then used within the streamtube model to estimate the blade forces.

Since the vortex model is time-dependent, the dynamic stall model is used

every time-step. Hence, the inputs to the dynamic stall model are from the

current and the previous time-steps. Note, that the changes of flow velocity

due to the vortex contribution and the effective angle of attack are calculated

within the vortex model. Thus, the attached flow part of the dynamic stall

model (equation (3.35)) is omitted when combined with the vortex model.

4. Measurements of aerodynamic forces

This chapter summarizes the force measurements on the 12 kW VAWT with the diameter of 6.48 m which is operated at an open site north of Uppsala, Sweden (N 59

55

32”, E 17

35

12”). The measurement campaign is described in detail in Paper IV and the experimental data on the turbine forces is pub- lished in Paper V.

4.1 Experimental setup

The measurements of the forces were conducted on the VAWT designed and constructed at Uppsala University in 2006. The measurements of the power coefficient on this turbine were conducted in 2009. The results have shown that the power coefficient had the maximum value of 0.29 at the TSR of 3.3 [21]. In 2014, this turbine was renovated and equipped with force sensors, see figure 4.1. Due to the installed force sensors, the radius of the turbine has increased from 3 to 3.2 m. The force sensors are single-axis load cells which measure tension and compression and are installed between the turbine hub and support arms of one blade. Spacers of equal radial distance and weight are installed between the hub and the support arms of other two blades. The load cells assembly is depicted in figure 4.2. The specifications of the turbine with installed load cells are listed in table 4.1.

4.2 Data treatment

Four load cells are used to measure forces on one blade and its support arms.

The sum of the measured forces F

, F

, F

and F

represents the radial force.

Table 4.1. Nominal parameters of the 12 kW VAWT used for the force measurements Number of blades 3

Turbine radius 3.24 m

Hub height 6 m

Blade length 5 m

Blade profile NACA0021

Chord length 0.25 m

Blade pitch angle 2

Figure 4.1. The 12 kW turbine, designed and built at the Division of Electricity at Uppsala University. The turbine is equipped with load cells used for the force mea- surements. The torbine location is N 59

55

32

, E 17

35

12

.

The turbine rotational speed can be kept constant by controlling the electric load of the generator, see Ref. [21] and Paper IV. When the speed of the tur- bine is nearly constant, it can be assumed that the rate of change in angular momentum approaches zero, and centrifugal force on the rotor is constant.

The normal force is the difference between the radial force and the centrifugal force

F

= F

+ F

+ F

+ F

− F

, (4.1) and the centrifugal force F

is estimated as

F

= mΩ

L

, (4.2)

where m is the mass of the blade and support arms and L

is center of mass of the blade and support arms. The tangential force on the blade and support arms is related to the difference between F

+ F

and F

+ F

, and based on the assumption of the constant centrifugal force, it can be estimated as

F

= L

2L

(F

+ F

− F

− F

) . (4.3)

Dimensions L

, L

and L

in equations (4.2) and (4.3) are illustrated in fig-

ure 4.2 and specified in table 4.2.

R

L

center of mass of blade and support arms

Load cells setup F

F

L

(a)

F

F

F

F

L

spacers

(b)

Figure 4.2. The load cells assembly on the 12 kW VAWT used for the force measure-

ments. (a) Notation of the normal and tangential force; (b) Notation of the measured

loads.