Fluid Mechanics of Vertical Axis Turbines: Simulations and Model Development

(1)

UNIVERSITATISACTA UPSALIENSIS

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

Fluid Mechanics of Vertical Axis Turbines

Simulations and Model Development

ANDERS GOUDE

ISSN 1651-6214 ISBN 978-91-554-8539-9

(2)

Dissertation presented at Uppsala University to be publicly examined in Polhemssalen, Ångströmslaboratoriet, Lägerhyddsvägen 1, Uppsala, Friday, December 14, 2012 at 13:15 for the degree of Doctor of Philosophy. The examination will be conducted in English.

Abstract

Goude, A. 2012. Fluid Mechanics of Vertical Axis Turbines: Simulations and Model Development. Acta Universitatis Upsaliensis. Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 998. 111 pp. Uppsala.

ISBN 978-91-554-8539-9.

Two computationally fast fluid mechanical models for vertical axis turbines are the streamtube and the vortex model. The streamtube model is the fastest, allowing three-dimensional modeling of the turbine, but lacks a proper time-dependent description of the flow through the turbine. The vortex model used is two-dimensional, but gives a more complete time-dependent description of the flow. Effects of a velocity profile and the inclusion of struts have been investigated with the streamtube model. Simulations with an inhomogeneous velocity profile predict that the power coefficient of a vertical axis turbine is relatively insensitive to the velocity profile. For the struts, structural mechanic loads have been computed and the calculations show that if turbines are designed for high flow velocities, additional struts are required, reducing the efficiency for lower flow velocities.Turbines in channels and turbine arrays have been studied with the vortex model. The channel study shows that smaller channels give higher power coefficients and convergence is obtained in fewer time steps. Simulations on a turbine array were performed on five turbines in a row and in a zigzag configuration, where better performance is predicted for the row configuration. The row configuration was extended to ten turbines and it has been shown that the turbine spacing needs to be increased if the misalignment in flow direction is large.A control system for the turbine with only the rotational velocity as input has been studied using the vortex model coupled with an electrical model. According to simulations, this system can obtain power coefficients close to the theoretical peak values. This control system study has been extended to a turbine farm. Individual control of each turbine has been compared to a less costly control system where all turbines are connected to a mutual DC bus through passive rectifiers. The individual control performs best for aerodynamically independent turbines, but for aerodynamically coupled turbines, the results show that a mutual DC bus can be a viable option.Finally, an implementation of the fast multipole method has been made on a graphics processing unit (GPU) and the performance gain from this platform is demonstrated.

Keywords: Wind power, Marine current power, Vertical axis turbine, Wind farm, Channel flow, Simulations, Vortex model, Streamtube model, Control system, Graphics processing unit, CUDA, Fast multipole method

Anders Goude, Uppsala University, Department of Engineering Sciences, Electricity, Box 534, SE-751 21 Uppsala, Sweden.

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

(3)

To my family

(4)

(5)

List of papers

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

I Goude, A., Lundin, S., Leijon, M., “A parameter study of the influence of struts on the performance of a vertical-axis marine current turbine”, In “Proceedings of the 8th European wave and tidal energy conference, EWTEC2009”, Uppsala, Sweden, pp. 477–483, September 2009.

II Goude, A., Lalander, E., Leijon, M., “Influence of a varying vertical velocity profile on turbine efficiency for a Vertical Axis Marine Current Turbine”, In “Proceedings of the 28th International Conference on Offshore Mechanics and Arctic Engineering, OMAE 2009”, Honolulu, USA, May 2009.

III Grabbe, M., Yuen K., Goude, A., Lalander, E., Leijon, M., “Design of an experimental setup for hydro-kinetic energy conversion”,

International Journal on Hydropower & Dams, 15(5), pp. 112–116, 2009.

IV Goude, A., Ågren, O., “Simulations of a vertical axis turbine in a channel”, Submitted to Renewable Energy, October 2012.

V Goude, A., Ågren, O., “Numerical simulation of a farm of vertical axis marine current turbines”, In “Proceedings of the 29th International Conference on Offshore Mechanics and Arctic Engineering, OMAE 2010”, Shanghai, China, June 2010.

VI Dyachuk, E., Goude, A., Lalander, E., Bernhoff, H., “Influence of incoming flow direction on spacing between vertical axis marine current turbines placed in a row”, In “Proceedings of the 31th International Conference on Offshore Mechanics and Arctic Engineering, OMAE 2012”, Rio de Janeiro, Brazil, July 2012.

VII Goude, A., Bülow, F., “Robust VAWT control system evaluation by coupled aerodynamic and electrical simulation”, Submitted to Renewable Energy, September 2012.

VIII Goude, A., Bülow, F., “Aerodynamic and electric evaluation of a VAWT farm control system with passive rectifiers and mutual DC-bus”, Submitted to Renewable Energy, November 2012.

(6)

IX Goude, A., Engblom, S., “Adaptive fast multipole methods on the GPU”, Journal of Supercomputing, DOI 10.1007/s11227-012-0836-0, In Press, October 2012.

Reprints were made with permission from the publishers.

The author has also contributed to the following paper, not included in the thesis:

A Yuen, K., Lundin, S., Grabbe, M., Lalander, E., Goude, A., Leijon, M., “The Söderfors Project: Construction of an Experimental Hydroki- netic Power Station”, In “Proceedings of the 9th European wave and tidal energy conference, EWTEC2011”, Southampton, United Kingdom, September 2011.

(7)

Nomenclature

A m² Turbine cross-sectional area

A_∞ m² Asymptotic area of streamtube enclosing turbine A_c m² Cross-sectional area of a channel

A_d m² Area of turbine disc

Ae m² Streamtube area far downstream/center of turbine AR − Aspect ratio of a blade

CD − Drag coefficient

C_Ds − Drag coefficient for strut

CD∞ − Drag coefficient for infinitely long blade C_L − Lift coefficient

CLs − Lift coefficient for strut

CL∞ − Lift coefficient for infinitely long blade CN − Normal force coefficient

CN0 − Normal force coefficient without curvature corrections C_Ns − Normal force coefficient for strut

CP − Power coefficient

CP_e − Power coefficient equivalent for extracted power CPmax − Maximum power coefficient for a given flow velocity

CT − Tangential force coefficient

C_{T s} − Tangential force coefficient for strut

D m Turbine diameter

F − Velocity correction factor

FD N Drag force

FD0 N Drag force at zero angle of attack FL N Lift force

FN N Normal force

F_Nl N/m Normal force per meter FR N Force in radial direction F_T N Tangential force

FT s N Tangential force on strut

Fx N Aerodynamic force from a blade on flow in a streamtube F_xs N Aerodynamic force from a strut on flow in a streamtube

H m Channel height

J kgm² Moment of intertia L m Distance between struts K N Constant to determine lift force

(10)

Ma − Mach number

N − Number of particles (FMM) N_b − Number of blades

Nbox − Number of boxes (FMM)

Nd − Number particles per box (FMM) N_p − Number of panels

Ns − Number of struts N_t − Number of turbines Nv − Number of vortices

NFi − Set of all boxes in near field of box i (FMM) P W Power absorbed by the turbine

Pe W Power extracted from the turbine

P_e_,tot W Total power extracted from the turbines of a farm Ptot W Power available in flow

R m Turbine radius

Rinner m Strut inner attachment point Router m Strut outer attachment point

Re − Reynolds number

Ts Nm Torque from strut V m/s Flow velocity

V0 m/s Flow velocity far upstream of turbine V_∞ m/s Asymptotic flow velocity

V_abs m/s Magnitude of incoming flow velocity at blade position V_b m/s Blade velocity

V_d m/s Flow velocity at turbine disc

Ve m/s Flow velocity far downstream/center of turbine

Vi m/s Vortex velocity

V_r m/s Relative flow velocity for a blade (absolute value) Vrs m/s Relative flow velocity for a strut (absolute value)

Vs m/s Flow velocity at strut position

Vre f m/s Reference flow velocity for estimating angle of attack Vrel m/s Relative flow velocity for blade

V_relz m/s Relative flow velocity for blade in its own reference frame Vs m/s Far downstream velocity of flow passing outside turbine Vs j m/s Flow velocity at strut segment

V_ω m/s Velocity due to vortices W m²/s Complex velocity potential

W_b m²/s Complex velocity potential for blade velocity a − Axial induction factor

a_i − Multipole coefficient i (FMM)

as − Slope of lift coefficient curve in CL/α plot b m Circle radius used for conformal mapping b_i − Local coefficient i (FMM)

(11)

b_y_∞ − Normalized asymptotic streamtube width

bye − Normalized streamtube width at the turbine center bz∞ − Normalized asymptotic streamtube height

bze − Normalized streamtube height at the turbine center

c m Blade chord

c_s m Strut chord

csound m/s Speed of sound

c₀ m Reference chord for struts

g m/s² Gravitational acceleration

h m Turbine height

k₁ kgm² First control system constant k₂ kgm² Second control system constant k₃ kgm² Third control system constant

kd1 − Constant for time estimate of direct evaluation (FMM) k_d2 − Constant for time estimate of direct evaluation (FMM)

l m Blade length in streamtube

p − Number of multipole coefficients (FMM) p0 N/m² Pressure far upstream of turbine

patm N/m² Atmospheric pressure

pd1 N/m² Pressure directly in front of turbine disc pd2 N/m² Pressure directly after turbine disc

pe N/m² Pressure far downstream of turbine

r0 m Box center (FMM)

rs m Radial position on a strut

r m Arbitrary position

ri m Vortex position

s m Position on blade surface in transformed plane

t s Time

tb m Blade thickness

td s Time estimate for direct evaluation (FMM) u − Interference factor

x m Position in the x-direction x₀ m Blade attachment point

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

y_∞ m Asymptotic streamtube position in y-direction y_d m Streamtube position at turbine disc in y-direction y_e m Streamtube position far downstream in y-direction

Δy m Streamtube width

z m Position on blade surface in the blades reference frame z0 m Position on blade surface in the turbines reference frame z_b m Blade position

Δz m Streamtube height

(12)

Γp m³/s Circulation of three dimensional point vortex Γ m²/s Circulation of two dimensional vortex Δ m Cutoff radius used for vortex merging Ω rad/s Turbine rotational velocity

Ωi rad/s Rotational velocity of turbine i

Ω1 rad/s First control system rotational velocity constant Ω2 rad/s Second control system rotational velocity constant Ω3 rad/s Third control system rotational velocity constant

α − Angle of attack

αb − Corrected angle of attack α^s − Angle of attack for strut

β − Direction of incoming wind δ − Blade pitch angle

ε m Cutoff radius of Gaussian vortex kernel η − Angle of blade relative to the vertical axis η^s − Angle of strut relative to the horizontal plane

θ − Blade azimuthal position shifted 90 degrees θ^b − Blade azimuthal position

λ − Tip speed ratio

λe − Equilibrium tip speed ratio

λ^max − The tip speed ratio that gives highest power coefficient ν m²/s Kinematic viscosity

ρ kg/m³ Density of fluid σ N/m² Stress

ϕ − Angle of relative wind φ m²/s Velocity potential ω 1/s Vorticity

(13)

Abbreviations

CPU Central processing unit

CUDA Compute unified device architecture FEM Finite element method

FMM Fast multipole method FVM Finite volume method

GPU Graphics processing unit L2L Local to local translation L2P Local to particle evaluation M2L Multipole to local translation M2M Multipole to multipole translation

NACA National Advisory Committee for Aeronautics P2M Particle to multipole initialization

P2P Particle to particle interaction (direct evaluation) RANS Reynolds-averaged Navier-Stokes equations

SIMD Single instruction multiple data SSE Streaming SIMD Extensions

(14)

1. Introduction

A turbine is used to convert the energy from a moving fluid into rotational motion, which in turn can drive an electric generator. The best suited turbine for this energy conversion depends on the characteristics of the flow. One case, which is present in gas turbines and traditional hydro power plants, is flow con- strained by walls, such as flow within a pipe. Here, the turbine cross-sectional area usually covers the entire flow and the flow has a large pressure difference that drives the turbine. A second case is the free flow, where no confining walls are present. This is the case for wind power and often a reasonable approxi- mation for tidal power in the ocean. In free flow, the fluid can pass around the turbine and the available energy is the kinetic energy in the flow. This thesis mainly treat the free flow case, but also a hybrid case where there are confining walls, but the turbine does not cover the whole cross-sectional area of the flow, is studied in this thesis. This situation occurs in a river, where the river cross-sectional shape area usually prevents a turbine from covering the entire cross-section.

1.1 Different turbine types

The typical turbine design for wind power is a horizontal axis turbine, where the rotational axis of the turbine is parallel to the flow direction [1, chapter 1].

In this thesis, however, the vertical axis turbine will be investigated. Here, the rotational axis is perpendicular to the flow direction. This kind of turbine is sometimes called “cross-flow turbine”, as the turbine in principle also can be tilted 90 degrees to have a horizontal axis while still having its rotational axis perpendicular to the flow. The traditional name “vertical axis turbine” will be used here, even for situations where the rotational axis is tilted, since this is the most commonly used name.

There are two different types of vertical axis turbines. The first type is based on the drag force and is often called the Savonius rotor after the Finnish inventor Sigurd Johannes Savonius, despite that Savonius only patented an improvement of older designs [2]. This improvement is neither implemented on all present drag-based turbines. Drag-based devices rely on variation of the drag coefficient with respect to the orientation of the object. To create a rea- sonably efficient drag-based turbine, the drag coefficient should be high in one direction and low in the opposite direction, which gives a torque on the turbine. Drag-based devices achieve lower power coefficients than the lift-based

(15)

Figure 1.1. Different types of vertical axis turbines.

devices described in the section below [3, chapters 2, 7]. Another drawback is that the amount of construction material in drag devices is quite high (as can be seen in figure 1.1). This cost is inhibiting the construction of large turbines, as material usage is proportional to the volume, i.e. the cube of the characteristic length of the turbine, while the power absorption is proportional to the cross-sectional area, i.e. the square of the characteristic length.

The second type is the lift-based turbine, which was originally invented by the French engineer George Jean Marie Darrieus [4] in the 1920’s (approximately one year after Savonius patented his design). The patent application of Darrieus covers both the curved blade turbine and the H-rotor (see figure 1.1), as well as turbines with varying pitch angle and ducted turbines. It is suggested in the patent that the designs work both for wind and tidal energy. The aim of the curved blade design is to reduce the bending stresses in the blades due to centrifugal forces. The North American company Flowind commercialized in the 1980s the Darrieus turbine with the curved blade design [5, chapter 1].

During that time, the curved blade turbine was also studied by Sandia National Laboratories, which is the main reason why much of the published work on Darrieus turbines is on the curved blade design. This thesis will instead focus on the straight blade H-rotor design, which currently is in development at Up- psala University. With recent progress for light materials, composites can be used in the turbine construction, which reduces centrifugal forces due to the

(16)

lighter structure. This makes the H-rotor design more feasible. The straight blade design has the advantages that straight blades are easier to manufacture and by attaching the blades with struts, it is possible to place the upper bear- ing much closer to the turbine center, reducing the bending moment on the axis. In addition, the constant radius of the straight blade design gives a larger cross-sectional area. Disadvantages compared to the curved blade design are the addition of extra struts and the higher bending moments due to centrifugal forces.

The main aerodynamic advantage of vertical axis turbines, compared to standard horizontal axis turbines, is the independence of flow direction, re- moving the need for a yaw mechanism. For water flow, an additional advantage is that the cross-sectional area can be more flexibly chosen as both height and diameter can be varied (and the diameter can vary with height). This can be useful in shallow water where a turbine with a large width and small height can cover a larger area than a horizontal axis turbine, as the cross-sectional area of a horizontal axis turbine is circular. Disadvantages of the vertical axis turbines are the lower power coefficients and that the turbines are typically not self-starting.

The vertical axis turbine can have its generator on the ground, which in the wind power case simplifies maintenance, tower construction and makes the weight of the generator less important. This is beneficial for direct driven generators, which typically have large diameters. The use of direct driven generators further reduces the number of moving parts in the system. One major concern for vertical axis turbines is the cyclic blade forces in each revolution, which leads to torque oscillations and material fatigue. For further comparisons between horizontal and vertical axis turbines (and also between curved and straight blade turbines) see e.g. [6].

1.2 Comparison between wind and marine current turbines

Even though wind turbines operate in air (gas) while marine current turbines operate in water (liquid), there are many similarities between the two. Tra- ditionally, water is considered an incompressible fluid and can therefore be modeled with the incompressible Navier-Stokes equations. For air, it is typically expected that compressibility effects can be neglected for Mach numbers Ma within the range

Ma= V_rel

csound < 0.3.

Here, Vrel is the relative flow velocity (measured in the blades rest frame) and c_soundis the speed of sound in the fluid [7, chapter 9]. Note that the major con- tribution to Vreloriginates from the blades’ own motion for lift-based turbines.

(17)

The speed of the blades in a wind turbine is typically too low for the Mach number to be above 0.3 and wind turbines can therefore also be modeled with incompressible aerodynamics. Although both wind and marine turbines can be studied with the incompressible Navier-Stokes equations, there are still some characteristic differences. One difference is that for marine current turbines, there is both a sea bed and a free surface that bounds the flow. Another difference is the risk for cavitation at too high flow velocities. Cavitation would modify the flow characteristics and can cause damage to the turbine [8].

The energy absorbed by a turbine is proportional to the fluid density and the cube of the flow velocity. As the density of water is 800 times higher than the density of air, comparatively low fluid velocities are adequate for marine current power generation. For equal cross-sectional area, a wind speed of 10 m/s has the same incoming kinetic power as a water flow speed of 1.1 m/s.

However, as the forces only are proportional to the square of the flow velocity, the marine current turbine experiences approximately 9.3 times higher fluid mechanical forces than a wind turbine at the same conditions (assuming that the turbines are identical) and rotates 9.3 times slower. The increased forces for marine current turbines require both stronger blades and support structure.

Many of the experimental vertical axis turbines for water have used relatively large blades and thereby low optimal tip speed ratios [9–12].

One important parameter for the effectiveness of the turbine is the Reynolds number, which for a blade is defined as

Re=cV ν ,

where V is the flow velocity, c is the blade chord andν is the kinematic viscosity. A higher Reynolds number usually decreases the drag losses and increases the stall angle, which is beneficial for vertical axis turbines. For 20 ^◦C, the kinematic viscosities are 15.1μm²/s for air and 1.00μm²/s for water [13, ap- pendix A]. Under the conditions of equal power extraction mentioned above (i.e. 9.3 times higher flow velocity for the wind turbine), this would give a 63 % higher Reynolds number for the marine current turbine, which is within the same order of magnitude as the wind turbine.

1.3 Vertical axis turbine research at Uppsala University

At the Division of Electricity at Uppsala University, three vertical axis wind turbines have been built. The first turbine had a cross-sectional area of 6 m² and was later followed by a turbine with the cross-sectional area 30 m²and the rated power 12 kW [14–16]. This larger turbine is used for most of the experi- ments. A 10 kW turbine for telecom applications has also been built [17]. Fur- ther, a 200 kW turbine has been constructed by the spin-off company Vertical Wind AB [18]. Additionally, a marine current turbine (described in paper III) is scheduled to be deployed by the end of 2012.

(18)

Several simulation tools for turbine simulations have previously been developed at the division. A two-dimensional inviscid vortex model based on conformal mappings for the blades has been created by Deglaire et al. [19]. In the turbine implementation, each blade is solved independently [20], allowing for coupling to an elastic method developed by Bouquerel et al., see paper IV in [21]. A multibody version for simulating turbines has been developed by Österberg et al., see [22] and paper III in [21]. Two streamtube models have also been implemented by Deglaire and Bouquerel.

1.4 Extended studies within this thesis

This thesis focuses on the fluid mechanical modeling of the vertical axis turbine, and two different simulation tools have been developed. The first simulation tool uses the streamtube model and the development of this tool started from the basic streamtube model implemented by Bouquerel, which is based on the model of Paraschivoiu [3]. All additional modeling and code development have been developed within this thesis.

The second simulation tool uses a vortex model, and this tool has been developed from scratch within this work. This model is based on empirical data for lift and drag coefficients instead of the conformal mapping method by Deglaire, which is based on inviscid theory. The computational speed is crucial for the developed vortex model and large efforts have been put into this. The existing implementation of the fast multipole method by Stefan Eng- blom [23] has been significantly improved and ported to a GPU (paper IX)

Several studies have been carried out with the two simulation models. The streamtube model has been used to study losses due to struts (paper I), the effects of a velocity profile (Paper II) and to design a turbine for deployment in a river (paper III). The more computationally demanding vortex model has been used to study turbines in channels (paper IV) and turbine arrays (paper V and VI). The vortex model is also coupled to an electrical model to study control systems for a single turbine (paper VII) and extended simulations analyze control systems for a turbine farm (paper VIII).

1.5 Outline of the thesis

After the introduction, theory for vertical axis turbines is presented in chapter 2. This is followed by an introduction to control systems in chapter 3.

The theory and implementation for the simulation models are then presented in chapter 4, which also includes the GPU implementation of the fast multipole method. The results from the simulations are given in chapter 5, where the first part evaluates the accuracy of the simulation models and the second part summarizes the results from the articles. The thesis ends by conclusions,

(19)

suggestions for future work, summary of papers, errata for papers and ac- knowledgments.

(20)

2. Theory for vertical axis turbines

Given a cross-sectional area A perpendicular to a homogenous flow of a fluid, the kinetic power that passes through this area is given by

Ptot=1

2ρAV³, (2.1)

whereρ is the density and V is the flow velocity. If the flow is not confined by any surrounding boundaries, the kinetic power is the available power for a wind/current turbine. The efficiency (i.e. outgoing power divided by incoming power) would be one possible measure of how good the energy conversion is. However, adding a turbine will change the velocity and force parts of the flow to pass outside the turbine area and thereby change the kinetic energy that passes through this area A. Moreover, some kinetic energy is left in the flow and can possibly be used later. Therefore, turbine performance is usually measured with the power coefficient instead, which is defined as

CP= P

12ρAV_∞³, (2.2)

where P is the power absorbed by the turbine and V_∞ is the asymptotic upstream flow velocity. With this expression, the absorbed power is compared to the power that would have passed through the cross-sectional area, if the turbine would be absent, instead of compared to the power that actually passes through the area. Since this expression is normalized against an expression that does not change with the turbine characteristics, it is a better measure than efficiency. Improving the power coefficient will give higher power absorption, which is not always the case with efficiency.

2.1 Basic theory and the Betz limit

One of the most basic approximations of a turbine is the one used in the traditional Betz theory [24], where the turbine is approximated as a single flat disc with a constant pressure drop over the whole turbine surface. All flow passing through the disc is encapsulated in a streamtube that starts far ahead of the turbine and ends far behind. By making the assumption that the pressure at both ends of the streamtube is the atmospheric pressure patm, and by using the

(21)

p_atm patm

V_∞ A_∞

V_d

p_d1 p_d2

V_e Ae

A_d

Figure 2.1. The streamtube used in the Betz limit derivation. The dashed line shows the control volume used for momentum conservation.

Bernoulli equation before and after the turbine p_atm+1

2ρV_∞² = pd1+1

2ρV_d², (2.3)

pd2+1

2ρV_d² = patm+1

2 ρV_e², (2.4)

combined with continuity

A_∞V_∞= AdVd= AeVe (2.5) and momentum conservation for the control volume in figure 2.1 (marked with the dashed line)

ρAeV_e²−ρA∞V_∞²= Ad(pd2− pd1), (2.6) it is possible to show that the velocity at the turbine disc Vd is equal to

V_d =V_∞+Ve

2 . (2.7)

Considering that the Betz theory assumes no losses, the power absorbed by the turbine is given as the difference between incoming and outgoing power in the fluid

P= 1

2ρA_∞V_∞³−1

2 ρAeV_e³. (2.8)

If the axial induction factor a, defined as

Vd= (1 − a)V∞, (2.9)

is combined with expression with equations (2.5) and (2.8), it can be shown that the power will reach its maximum value for a= 1/3, and the optimal power is given by

P=16 27·1

2ρA_dV_∞³ (2.10)

where 16/27 is the traditional Betz limit, limiting the power coefficient to approximately 59.3 %.

(22)

p₀ p₀ V0

V0

V₀

A₀ Ac

V_d p_d1 p_d2

Ve

A_e A_d

p_e Vs

Figure 2.2. Illustration of a streamtube confining the flow that passes through the turbine disc for a channel of cross-sectional area Ac.

2.2 Extension to include channels

The Betz theory assumes no outer boundaries in the system. For a turbine operating in water, it is more common that boundaries are present. One example is a river, where the flow is limited by the width and depth. The outer walls prevent the flow from expanding, pushing more flow through the turbine. This occurs in a traditional hydro power plant, where the entire flow is forced to pass through the turbine, which results in much higher power absorption than the Betz limit [25].

To analyze this case analytically, assume that the flow upstream of the tur- bine has constant velocity V0 and pressure p0 (see figure 2.2). Note that in this case, the pressure upstream and the pressure downstream are not equal.

Instead, there will be a drop in pressure, which, for open channel flow, would correspond to a drop in the surface level. Due to the continuity of the pressure, the pressure inside the streamtube and outside has to be the same downstream (pe). The cross-sectional area of the channel is Acand the cross-sectional area of the turbine is A_d. In this case, the Bernoulli equation gives

p0+1

2 ρV₀² = pd1+1

2ρV_d², (2.11)

p_d2+1

2 ρV_d² = pe+1

2ρV_e², (2.12)

p₀+1

2 ρV₀² = pe+1

2ρV_s², (2.13)

the continuity equation gives

A0V0= AdV_d = AeVe, (2.14) (Ac− A0)V0= (Ac− Ae)Vs (2.15)

(23)

and momentum conservation for a control volume that encloses the entire channel gives

ρA^eV_e²+ρA^sV_s²−ρAV0²= Ad(pd2− pd1) + Ac(p0− pe). (2.16) The velocity at the turbine can be derived from equations (2.11) – (2.16) as

V_d= Ve(Vs+Ve)

Vs+ 2Ve−V0, (2.17)

which in the free flow limit(Vs→ V0) reduces to equation (2.7). The force on the turbine becomes

Fx= Ad(pd1− pd2) = 1 2A_d

V_s²−Ve²

, (2.18)

giving the power as

P= FxVd= 1

2AdVe(Vs+Ve)

V_s²−Ve²

V_s+ 2Ve−V0 . (2.19)

From equation (2.19), it can be found that the highest power absorption is obtained when V_e= V0/3, which actually is the same as for the free flow. The highest power is thus given by

P=16

27· 1

1−^A_A^d_c2·1

2 ρA_dV₀³ (2.20)

and the pressure drop is

pe− p0=4^A_A^d_c

3−^A_A^d_c 9

1−^A_A^d_c2 ρV0². (2.21) From these results, it can be seen that the maximum theoretical power coefficient increases with the factor (1 − Ad/Ac)⁻² for a channel. In the limit Ad→ Ac, the power coefficient diverges, along with the pressure drop in equation (2.21). Considering that an infinite drop in pressure is unfeasible, when the turbine area is almost as large as the channel, the available pressure difference will start limiting the maximum power coefficient, which will prevent infinite energy extraction. This is the case for hydro power turbines, where the power is limited by the difference in water elevation.

The model above assumes that the cross-section of the channel is constant.

An open channel will have a drop in surface level over the turbine and an extension of the model to include this drop is given by Whelan et al. [26].

This correction has not been included in the present work, as the model is used for comparisons with the two-dimensional vortex simulations where no free surface is modeled.

(24)

R

θb

Ω

x y

Figure 2.3. Illustration of a vertical axis turbine.

2.3 Theory of lift-based vertical axis turbines

The H-rotor is a lift-based design, which means that the aerodynamic torque is generated by the lift force of the blades. Therefore, airfoil profiles are typically used for the blades. Airfoils are generally designed to operate at relatively low angles of attack, where the lift force increases approximately linearly with the angle of attack and the drag force remains low. When the angle of attack increases above the stall angle, which for a typical airfoil occurs for angles around 10 – 15 degrees, the lift force is reduced, and the drag starts to increase substantially. This work is focused towards blades with a fixed pitch angle. To keep the angle of attack low without pitching the blades, the blades must move with a high velocity if the wind comes from the side. To illustrate this, assume that a blade is located at angleθb (see figure 2.3). Using complex notation, this gives the blade position

z_b= Reⁱ^θ^b (2.22)

and the velocity of the blade is therefore

V_b= i ˙θbReⁱ^θ^b. (2.23) The rotational velocity ˙θbis commonly denotedΩ. The incoming wind V, if complex, represents wind from any direction. The blade will now see a relative wind of

Vrel= V − iΩReⁱ^θ^b. (2.24)

(25)

Vrel

α ϕ δ

V_b

Figure 2.4. Definitions of angles and velocities. The positive direction for angles is counter-clockwise, henceα ^andϕare negative for the directions of Vband Vrelin the figure.

To obtain the angle of relative wind, rotate V_rel with the angle ie^−iθ^b, aligning the blade motion with the negative real axis,

V_relz= Vie^−iθ^b+ ΩR. (2.25)

The angle of relative wind will be the argument of this complex number. In the special case of V being real, the angle of relative wind is

ϕ = arctan cosθb ΩRV + sinθb

(2.26) and the absolute value of the relative velocity is

|Vrel| = V

ΩR

V + sinθb

2

+ (cosθb)². (2.27) The angle of attack is given by

α = ϕ + δ, (2.28)

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

λ =ΩR

V_∞ (2.29)

where the asymptotic velocity is used. Due to the energy extracted, the velocity at the blade will generally be lower than the asymptotic velocity.

(26)

The turbine torque is calculated from the tangential force, which is given by

F_T = FLsinϕ − FDcosϕ. (2.30) When the angle of relative wind is low, the approximations sinϕ ≈ ϕ and cosϕ ≈ 1 can be applied. Assume that the pitch angle is zero, hence α = ϕ. For symmetric blades, the blade forces can be approximated as FL≈ Kα, where K is a constant, and FD≈ FD0, where FD0is a constant. The tangential force can therefore be estimated as

FT= Kϕ²− FD0, (2.31)

showing that when the angle of relative wind decreases, the drag force becomes dominating. The conclusion is that for high tip speed ratios, drag will give a more significant contribution, reducing the power coefficient. At too low tip speed ratios, the turbine will enter stall, where lift decreases and drag increases, which also should be avoided. For these reasons, the turbine should be designed to operate with a tip speed ratio close to the stall limit in order to obtain the highest possible power coefficient.

2.4 Angle of attack including flow curvature

The expression (2.26) is only valid for infinitely small symmetric blades. The blade performs a rotational motion, which leads to additional curvature effects, changing the effective angle of attack. To conclude the theory section, a more proper derivation will be performed using a rotating flat plate instead.

By the use of conformal mappings, a circle can be transformed into a flat plate with the Joukowski transformation. The s-plane represents the circle with radius b and the z-plane a flat plate extending between−2b to 2b giving the blade chord c as c= 4b. The blade coordinates z in its own frame of reference is given by

z= s +b²

s . (2.32)

Using the same transformation as Deglaire [20], the z₀plane can be defined as z0=

(z + x0)e^−iδ+ iR

eⁱ^θ = (z + D)eⁱ^(θ−δ), (2.33) where D= x0+ iReⁱ^δ. By assuming that the blade only rotates around the center, the blade velocity is given by

Vb = i ˙θ

(z + x0)e^−iδ+ iR

eⁱ^θ = iΩz0 (2.34) withΩ = ˙θ. Introduce a complex velocity potential W, with complex conju- gate W , such that

dW

dz₀ = V, (2.35)

(27)

hence the potential can be used to calculate the flow velocity V . Assume that the potential is given on the form

W(s) = Vabseⁱ^{(−β+θ−δ)}s+Vabse^{−i(−β+θ−δ)}b² s − iΓ

2πlog(s) +W1(s), (2.36) where V_abseⁱ^β is the flow velocity at the blade position. Now, construct a potential such that

dW_b

dz₀ = Vb, (2.37)

which gives dWb

ds = Vbdz0

ds = −iΩz0dz0

ds = −iΩ

z+ D dz

ds. (2.38)

Note that on the boundary, z= z. Integrated, this is

Wb = −iΩ

1 2z²+ Dz

= −iΩ 1

2

s²+ 2b²+b⁴ s²

+ D

s+b²

s

. (2.39)

The no-penetration boundary condition states that the stream function should be constant (possibly time-dependent) on the boundary. Given that the boundary is moving, the condition becomes

Im[W (s)−Wb(s)] = C. (2.40)

The first part of W(s) already fulfills the condition, while W1(s) remains to be determined, hence

Im[W1(s) −Wb(s)] = C. (2.41) The boundary condition at infinity states

dW1

ds _s

→∞= 0, (2.42)

and on the boundary,

s= b²

s (2.43)

applies. Write equation (2.41) in terms of complex conjugates iC = W1(s)−W1(s)+ iΩ

1 2

s²+ 2b²+b⁴ s²

+ 1

2

s²+ 2b²+b⁴ s²

+ D

s+b²

s

+ D

s+b²

s

= W1(s)−W1(s)+

iΩ

b⁴ s² +b⁴

s² +b² s b+b²

s D+b²

s D+ Db² s + 2b²

. (2.44)

(28)

Note that constants can be excluded from the potential. Therefore, W₁(s) can be identified as

W1(s) = −iΩ b²

s

D+ D +b⁴

s²

= −iΩ

b² s

iR

eⁱ^δ− e^−iδ + 2x0

+b⁴ s²

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

dW ds

_s

=b = 0 ⇒

−iΓ 2π

1

b = −Vabseⁱ^{(−β+θ−δ)}+Vabse^{−i(−β+θ−δ)}− iΩ

iR

eⁱ^δ− e^−iδ

+ 2x0+ 2b

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

− iΓ 2π

1

b = Vre f2isinα ⇒

sinα = −Vabssin(−β + θ − δ) − Ω(−Rsinδ + x⁰+ b)

Vre f , (2.47)

and that for small values of x0

Vre f ≈

(Vabscos(θ − β) + ΩR)²+V_abs² sin²(θ − β). (2.48) Now, redefine x0 in terms of the chord as x0= x0rc, which means that x0r= 0.25 is the quarter chord position and use that 4b = c

sinα = −Vabssin(θ − β − δ) − Ω

−Rsinδ + x^0rc+₄^c V_{re f}

= (Vabscos(θ − β) + ΩR)sinδ −Vabscosδ sin(θ − β)

V_{re f} −

Ω

x_0rc+^c₄ Vre f

= sin

δ + arctan −Vabssin(θ − β) V0cos(θ − β) + ΩR

−Ω

x0rc+₄^c

Vre f . (2.49) Assuming small angles of attack, one can approximate

arcsin(α + β) ≈ α + β (2.50)

(29)

which gives the simplifications

α = δ + arctan −Vabssin(θ − β)

V_abscos(θ − β) + ΩR−Ωx0rc V_{re f} − Ωc

4V_{re f}. (2.51) Note that at the positionθ = 0, the blade is at the position θb=π/2. The substitutionθ = θb−π/2 gives

α = δ + arctan V_abscos(θb−β)

Vabssin(θb−β) + ΩR−Ωx0rc Vre f − Ωc

4Vre f, (2.52) which with β = 0 would correspond to equations (2.26) and (2.28), but in- cludes mounting position x0rand flow curvature. As an example, for a turbine with chord 0.25 m and radius 3 m (i.e. the experimental turbine in Marsta [16]), at very high rotational velocities

Vre f ≈ ΩR

, the change in angle of attack due to flow curvature is approximately -1.2 degrees. This gives higher angles of attack upstream and lower downstream.

(30)

3. Control strategy for vertical axis turbines

Optimizing the power from a wind/marine current turbine does not only require that the turbine is designed with the highest possible power coefficient.

Another important factor is to make sure that the turbine actually runs at the tip speed ratio associated with the peak power coefficient. Therefore, a control strategy which keeps the turbine tip speed ratio near this optimal value is preferable.

Wind turbines in general are either controlled by pitch or stall regulation, where the most common design today is a horizontal axis turbine with pitch regulation [28]. The advantage with pitch control is that it introduces an additional parameter that can be controlled, allowing for a more flexible control system. Pitch control is mainly used in the region above rated wind speed (see figure 3.1) to keep a smoother power and reduce mechanical loads, as the pitch angle can be changed to reduce the blade forces [1, Chapter 8]. Stall regulation, instead, reduces the tip speed ratio, which increases the angle of attack.

This will eventually cause stall, which reduces the lift force and increases the drag force. This will be most prominent for the tangential force, and thereby the turbine torque, due to the significant increase in the drag force. Pitch control has been used for vertical axis turbines, mainly to improve performance at low tip speed ratios [29,30] where stall is avoided by actively altering the pitch angle to reduce the angle of attack. The angle of attack oscillates between positive and negative values as the blade moves between the upstream and downstream section of the turbine. Hence, reducing the angle of attack with active pitch requires a change in pitch angle during each revolution. An active pitch mechanism would complicate the turbine further. No pitch mechanisms are included in the turbines studied here to reduce the sources of mechanical failure.

Without a pitch mechanism, the remaining parameter to control is the rotational velocity, where the turbine power is controlled by regulating the tip speed ratio and thereby the power coefficient. Even with a pitch mechanism installed, it is common for horizontal axis wind turbines to use a fixed pitch angle in the variable rotational speed region illustrated in figure 3.1 [1, Chap- ter 8]. The following sections will focus on the variable rotational speed region, where the aim is to maximize the extracted power.

3.1 Control of a single turbine

One way to control a turbine is to perform real time flow velocity measurements and adjust the tip speed ratio to optimal values. However, this would,

(31)

0 2 4 6 8 10 12 14 16 0

0.5 1 1.5

Fractionofratedpower

Wind speed (m/s) Cut in

Variable rotational speed region, constant CP

Constant power Constant rotational

speed Rated wind speed

Figure 3.1. Example of the different control strategies for different flow velocities.

The turbine starts operating at wind speed 3 m/s and operates at optimal power coefficient (approximately constant tip speed ratio) up to rated rotational velocity (wind speed 10 m/s), where the rotational velocity is kept constant until rated power is achieved (wind speed 12 m/s). At higher wind speeds, power is kept constant. There is also a cut-out wind speed where the turbine is stopped (not shown in the figure).

Rotational velocity

Power

Speed increases Turbine stops Stable region

Equilibrium PeP P_e Speed

decreases

Speed decreases

Figure 3.2. Illustration of a control system, where the extracted power only depends on the rotational velocity.

(32)

rely on the accuracy of the flow measurements. An alternative approach is to let the extracted power Pe be a function only of the rotational velocity of the turbine. This type of control has been used for horizontal axis wind turbines [31], but due to the low power coefficient at low tip speed ratios for vertical axis turbines, some care has to be taken when transferring this control system to a vertical axis turbine to avoid that the turbine ceases to rotate for low tip speed ratios.

An example of a control system only using the rotational speed as input parameter is illustrated in figure 3.2. The angular acceleration ˙Ω of the turbine is

˙Ω = P− Pe

JΩ , (3.1)

where J is the moment of inertia of the system. Here, extracted power P_e includes both power from the generator and electrical and mechanical losses in the system. If the extracted power Peis less than the turbine power P, the turbine accelerates, while if the extracted power is larger, the rotational velocity decreases. In Figure 3.2, the turbine would stop if the turbine has too low rotational velocity as the turbine power at this low rotational velocity is very limited(P < Pe). This region where the turbine will stop is characterized by a low power coefficient at low tip speed ratios, which may apply to a vertical axis turbine. The existence of such a region depends on the power absorption characteristics of the turbine and for high enough wind speeds, this region typically becomes smaller.

Figure 3.2 shows that there will be an equilibrium where extracted power equals turbine power(P = Pe). If this equilibrium occurs at the peak of the power curve in figure 3.2, maximum energy is extracted. If the extracted power is normalized the same way as turbine power, the extracted power coefficient is

C_P_e = P_e

1

2ρAV_∞³. (3.2)

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

where ˙Ω = 0

CP_e= CP. (3.3)

Equation (3.2) can be written in terms of tip speed ratio and rotational velocity as

CP_e= P_eλ³

12ρA(RΩ)³. (3.4)

Denote λ^max the tip speed ratio with the peak power coefficient CPmax. By choosingλmaxas the desired equilibrium, equations (3.3) and (3.4) give

C_Pmax= P_eλmax³

12ρA(RΩ)³ ⇒ Pe= 1

2ρAC_Pmax

RΩ λmax

3

= k2Ω³, (3.5)

Fluid Mechanics of Vertical Axis Turbines: Simulations and Model Development