• No results found

Wind Turbine Using a Vortex Method

N/A
N/A
Protected

Academic year: 2021

Share "Wind Turbine Using a Vortex Method "

Copied!
42
0
0

Loading.... (view fulltext now)

Full text

(1)

UPTEC

Examensarbete 20 p Oktober 2010

Multi-Body Unsteady Aerodynamics in 2D Applied to a Vertical-Axis

Wind Turbine Using a Vortex Method

David Österberg

(2)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0

Postadress:

Box 536 751 21 Uppsala

Telefon:

018 – 471 30 03

Telefax:

018 – 471 30 00

Hemsida:

http://www.teknat.uu.se/student

Abstract

Multi-Body Unsteady Aerodynamics in 2D Applied to a Vertical-Axis Wind Turbine Using a Vortex Method

David Österberg

Vertical axis wind turbines (VAWT) have many advantages over traditional Horizontal axis wind turbines (HAWT).

One of the more severe problem of VAWTs are the complicated aerodynamic behavior inherent in the concept. In

contrast to HAWTs the blades experience varying angle of attack during its orbital motion. The unsteady flow

leads to unsteady loads, and hence, to increased risk for problems with fatigue.

A tool for aerodynamic analysis of vertical axis wind turbines has been developed.

The model, a Discrete vortex method, relies on conformal maps to simplify the task to finding the flow

around cylinders. After the simplified problem has been solved with Kutta condition, using the Fast Fourier

transform, the solution is transformed back to the original geometry yielding the flow about the turbine.

The program can be used for quick predictions of the aerodynamic blade loads for different turbines allowing the method to be validated by comparing the predictions to experimental data from real

vertical axis wind turbines. The agreement with experiment is good.

ISSN: 1401-5757, F10054

Examinator: Tomas Nyberg

Ämnesgranskare: Hans Bernhoff

Handledare: Paul Deglaire

(3)

Sammandrag

Vertikalaxlade vindkraftverk har m˚anga f¨ordelar ¨over de traditionella horisontalaxlade. Ett stort problem ¨ar den komplicerade str¨omning som uppst˚ar genom s˚adana vindkraftverk. Till skillnad fr˚an horisontalaxlade vindkraftverk ¨ar aerodynamiken per automatik ostadig eftersom bladen upplever varierande anfallsvinkel under sin cirkelbana. Det ostadiga fl¨odet leder till problem med v¨axlande laster vilket kan leda till problem med utmattning.

Ett verktyg f¨or simulering av flera viktiga aerodynamiska effecter hos vertikalaxlade vindkraftverk pre- senteras. Metoden, som tillh¨or klassen diskreta vortexmetoder, anv¨ander conforma avbildningar f¨or att f¨orenkla problemet till fl¨ode runt cylindrar. Efter att det f¨orenklade fl¨odesproblemet l¨osts med Kuttas villkor, m.h.a. Fast Fourier transformen, erh˚alles fl¨odet i den komplicerade geometrin genom en invers conform avbildning.

Programmet har t ex. kunnat anv¨andas f¨or att snabbt f¨oruts¨aga aerodynamiska krafter p˚a olika tur-

biner. Modellen har d¨armed kunnat valideras genom j¨amf¨orelse med experimentella data f¨or vertikalaxlade

vindturbiner. ¨ Overensst¨ammelsen anses vara god.

(4)

Contents

1 Introduction 7

1.1 The Status of Vertical Axis Wind Turbine Technology . . . . 8

1.2 The Goal of This Project . . . . 9

1.3 Previous Experimental Work . . . . 9

1.4 Previous Aerodynamic Prediction Models . . . . 9

1.5 Previous Multibody Potential Flow Solutions . . . . 11

2 Aerodynamic Assumptions 12 2.1 Two-dimensional Approximation . . . . 13

2.2 Inviscid Approximation . . . . 14

2.3 Incompressible Approximation . . . . 14

2.4 Boundary Layer Considerations . . . . 14

3 Flow Around a Moving Rotor 16 3.1 Vortex Dynamics . . . . 16

3.2 Potential Flow . . . . 17

3.3 Boundary Conditions . . . . 18

3.4 Conformal Mapping . . . . 19

3.5 Boundary conditions in Circle plane . . . . 20

3.6 Construction of Complex Flow Potential . . . . 20

3.7 Convection of Vortex Particles . . . . 23

3.8 Routh’s Rule . . . . 24

3.9 Convection in Circle Plane . . . . 25

4 Numerical Validation of the Model 26 4.1 Steady flow . . . . 26

4.2 Blade Loads . . . . 26

4.3 Strickland/Klimas Turbine . . . . 28

4.4 Strickland/Oler Turbine . . . . 28

5 Application to the Marsta 12kW Experimental Wind Turbine 31 6 Conclusions 36 A Method of Images 38 A.1 The Circle Theorem . . . . 38

A.2 Image of a Vortex . . . . 39

(5)

Nomenclature

β Azimuthal angle of the rotor, as measured to a reference blade.

`, ` Tangent vector to ∂ T . Function of ∂ T ω Vorticity ∇ × v [s −1 ]

β ˙ Angular velocity of the rotor. ˙β = 2λV D

∂ T Boundary between the turbine and the fluid [m]

γ Circulation around a vortex [m 2 /s]

γ k Vortex particle strength [m 2 /s]

ˆ

s n Point in Ω s in the visinity of s te,n where vortices are shed

λ, TSR Tip-speed-ratio of the rotor. Defined as the ratio of the undisturbed wind speed and the speed of the fastest rotor blade.

T Stress tensor of the general Navier-Stokes equations v Velocity field, often assumed two-dimensional [m/s]

x k , z k Vortex particle location [m]

x Radius vector, a point in the plane [m]

ˆn, n The unit normal vector to ∂ T . Function of ∂ T ˆn n , n n The unit normal vector to ∂ T n . Function of ∂ T n

Ω Region of the plane containing the fluid [m 2 ] Ω s Area outside N disjoint discs in the plane

∂Ω The boundary of the fluid [m]

φ Velocity potential v = ∇φ. [m 2 /s]

ψ The stream function, i.e., the harmonic conjugate of φ ρ Mass density [kg/m 3 ]

σ Solidity, σ = Nc/R

T Region in the plane occupied by the turbine [m 2 ] T n Region in the plane occupied by a turbine part [m 2 ]

T sn Disc in the Circle plane corresponding to T n in the Physical plane ˆz Unit normal vector to the plane of flow [m]

A Projected frontal area of the turbine

(6)

b n Radius of disk T sn

C P Power Coefficient, measure of aerodynamic efficiency C P =

1

P

2

ρv

3

A

c nk Coefficients of series describing the conformal map f D Diameter of the rotor.

D/Dt Convective derivative, equivalent to ∂/∂t + v · ∇ f Conformal map f : Ω s → Ω

N Number of disjoint parts of the turbine section N 0 Number of terms kept in truncated series describing f p Fluid pressure [N/m 2 ]

R Characteristic core radius [m]

s nk Center of disk T sn

s te,n Point on ∂ T n where the tangential velocity is zero V Undisturbed wind speed

v s Velocity in Circle plane, corresponding to v in the Physical plane HAWT Horizontal Axis Wind Turbine

VAWT Vertical Axis Wind Turbine

(7)

Chapter 1

Introduction

A wind turbine is a machine that converts the kinetic energy in the air into mechanical energy [1]. If the mechanical energy gained results in a rotational movement around a vertical axis the turbine is called a Ver- tical axis wind turbine or a VAWT (see figure 1.1). There have been many, more or less successful attempts at harnessing the power of the wind using the vertical axis concept. Today, however, the overwhelming majority of commercial turbines are of the competing Horizontal Axis Wind Turbine (HAWT) type.

(a) Danish HAWT (b) Giromill VAWT (c) Darrieus VAWT

Figure 1.1: The three main types of wind turbines

However, the vertical axis principle often simplifies the turbines. In theory that could lead to better economy. The advantages most commonly mentioned are 1) Generator can be placed on the ground where mainainence is simpler and where weight and aerodynamic properties are irrelevant. 2) The turbine is insensitive to the wind direction. Hence eliminating the need for a yawing mechanism. 3) If the blades are straight as in figure 1.1(b) the blades experience the roughly the same type of flow throughout their span.

Hence the blades can have a uniform blade profile, and 4) The primary fatigue loads are aerodynamic loads instead of gravitational loads. This is a strong argument for VAWTs for large scale machines.

In one sense, the price paid for structural simplicity is aerodynamic complexity. No simple theory exists that can adequately account for the complicated interaction of the rotor and it’s wake or between the rotor and the tower. For straight bladed VAWTs of the type displayed in figure 1.1(b) one also have to investigate the interaction between the rotor, the wake and the wing-tip vortices. These problems complicate the design process of VAWTs in comparison with HAWTs.

The power produced by a wind turbine is absorbed from the kinetic energy in the wind. It is; thus, proportional to the projected frontal area A of the turbine. Because the terminology is adopted from HAWT aerodynamics this is sometimes called the swept area of the turbine. It is obvious that a wind turbine cannot capture all kinetic energy of the wind. If so the air would come to a standstill behind the turbine. Hence there exists an upper limit for the aerodynamic efficiency. Often the so called Betz’ limit of (59.3%) is mentioned as this maximum. The technical name for aerodynamic efficiency of a wind turbine is the power coefficient defined as

C P = P

1 ρV 3 A (1.1)

(8)

where P is the absorbed mechanical power, ρ is the density of air (1.225 kg/m 3 , V is the wind speed and A is the projected frontal area of the turbine.

Modern HAWTs has benefited much from research and experience; they have aerodynamic efficiencies close to 50% where as VAWTs never exceeded 40%. Even so, there is no decisive argument showing that VAWTs must have less efficiency than HAWTs. On the contrary, the VAWTs sweep the area twice, opening a slight possibility of reaching C P s exceeding the Betz limit. On the other hand, there are sections of the lap where the blades do not produce a positive torque. Furthermore the blades are obliged to move through the turbulent wake on the downwind part of a revolution.

Efficiency aside, equally important is to build turbines that can function properly for the intended life- time of the device. Specifically this means designing turbines that show less fluctuating loads on the blades and shaft. For example the curved blades of Darrieus vertical axis wind turbines were known to fail from fatigue as early as two years after construction. The emergence of modern materials has somewhat re- lieved this situation. Never-the-less decreasing variance in the loads leads to improved reliability and better economy.

When designing VAWTs it is hence important to correctly understand the structure of the wake and how it corresponds to the blade loads.

(a) Giromill (b) Ropatec (c) Savonius

Figure 1.2: A few types of VAWT

1.1 The Status of Vertical Axis Wind Turbine Technology

During the last thirty years both the horizontal-axis and the vertical-axis turbines have undergone extensive development. It should be noted that vertical axis machines have received much attention in the academic community. The horizontal axis concept, on the other hand, has been successful in attracting commercial research and development.

The largest research effort to date was undertaken by Sandia National Laboratories in the USA [2].

Sandia committed themselves to the curved bladed Darrieus VAWT concept (depicted in figure 1.1(c)) and routinely built and studied these turbines for over fifteen years. In parallel the National Research Council of Canada were investigating a similar VAWT concept. The effort of Sandia and NRC resulted in a boom of interest in VAWTs in North America. Large commercial turbines were produced in by companies such as FloWind, Vawtpower, Indal Tech., Lavalin and Adecon Inc.

When the VAWT trend reached Europe most companies opted for straight bladed giromill (depicted in figure 1.1(b)). Foremost by British company VAWT Ltd. and German company Heidelberg. Although these companies discontinued their VAWT projects there has been recent development in connection to the growing market for wind turbines due to the current climate concerns. Notably a new niche for small scale VAWTs in city environment has emerged. For example, British company Quiet Revolution and dutch company Turby are prominent in this market. Another interesting new type of VAWT is the extreme high solidity turbines introduced by Italian Ropatec that has a cross section as depicted in figure 1.2(b). These companies are emphasizing two properties of the VAWT. 1) The VAWTs are able to absorb power even in the turbulent winds encountered in city landscapes. 2) The low tip-speed of the blades makes the VAWT silent compared to HAWT.

Indeed this diploma work is a part of a current research project at Uppsala University that is trying

to reinvent the VAWT concept by emphasizing simplicity, reliability and especially longterm economic

viability.

(9)

1.2 The Goal of This Project

In this diploma work the goal is to to develop an aerodynamic model capable of predicting the behavior of the turbine wake and the instantaneous aerodynamic forces acting on the turbine parts. Furthermore, it is desired that the model is flexible and fast enough to be used actively in the design process of VAWTs.

1.3 Previous Experimental Work

Given the goals of this project the most relevant studies are those that can be used to validate the predictions of the numerical model. Hence the most relevant experimental studies are based on setups able to:

1. Measure the pressure distribution over the blades.

2. Measure the force acting on the blades.

3. Measure the torque acting on the shaft.

4. Measure the power produced.

5. Visualize the structure of the wake.

Experiments of type 1-4 provides information of the same phenomena with decreasing level of detail (i.e., the turbine performance). Experiments of type 5 provides a general insight into the aerodynamic behavior and can be used when engineering wind farms.

Since the 1970s numerous vertical axis turbines have been built and tested in wind tunnels. However, there are large challenges in producing reliable and detailed performance data. To emulate free streaming conditions in a wind tunnel it is necessary that the model is kept small. In turn, this require that the angular frequency of the model must be kept very high resulting in unwanted vibrations that simultaneously distorts the flow and gives noise to the measurements. Further, it can be assumed that many of the results are kept as trade secrets of the companies who performed them.

The first wind tunnel measurements on Darrieus turbines were performed by the National Research Council (NRC) of Canada in 1972. NRC then proceeded to measure the performance of several turbines [3].

The experimental data that are most often cited were produced as a part of Sandia’s extensive wind project that over several years routinely built and tested Darrieus wind turbines. In particular, Strickland [4]

performed the first (and possibly the only) published flow visualization of the developed wake.

Several authors have measured VAWT instantaneous blade forces. Tow tank studies done at Texas Tech University presented measurements of blade forces on a straight bladed giromill both during operating conditions (tip-speed-ratio 5) and during dynamic stall condition (tip-speed-ratio 2.5) (described in [5]).

The Texas Tech study also measured the velocity profile of the wake. Blade force measurements on a straight bladed VAWT operating in the dynamic stall regime were also produced by Vittecoq et al. [6]

British company VAWT Ltd produced several full scale H-rotors equipped with systems for instanta- neous blade force measurements [7]. These results are interesting because the size of the turbine reduces the uncertainty in the force measurements. On the other hand the uncertainty in the flow conditions cannot be controlled as in wind tunnel measurements.

An in-depth review of the experimental work to date can be found in Paraschivoiu[5].

1.4 Previous Aerodynamic Prediction Models

There has been, in the past, several attempts of modeling lift-driven VAWTs each with their own advantages and disadvantages. The different methods can be classified in six groups:

1. Analytical maximum C P predictions 2. Fixed-wake vortex models

3. Streamtube models 1

1

In the HAWT community these models are often called Blade Element Momentum (BEM)Methods

(10)

4. Mesh-based CFD 5. Free-wake vortex models

The main fully analytical C P prediction model for vertical axis wind turbines is the double actuator disc model. It models the flow as one steady streamtube that passes through the turbine. One can think of it as two Betz turbines [8, 9] in tandem with a certain spacing in-between. The model was conceived by Templin in 1974 and was the first aerodynamic model of VAWTs [10]. Newman showed that the maximum C P ,max of such a double actuator-disc system is C P ,max = 16/25 = 64% [11].

A semi-analytical attempt is the so called fixed-wake vortex model developed by Holme[12] and ex- tended by Fanucci et al. [13]. This model assumes an infinite number of wings but with solidity (σ = Nc/R) fixed to a certain value. At each azimuthal position there is a wing with a certain circulation around. This circulation is calculated from the local angle of attack. Hence there is a vortex sheet bound to the turbine.

There is also vortex sheets, which leaves the turbine due to the change in circulation between adjacent wings, and are modeled to be convected downstream at a uniform velocity. The forces on the blades can be calculated from Kutta-Joukowski principle and integrated to obtain the performance of the turbine. In most fixed-wake analysis the maximum predicted C P is about 54% [14].

The most popular class of prediction models are streamtube models. There is a range of models that vary in how detailed the analysis is. The most sophisticated streamtube model is the Double Multiple Streamtube model (DMST) first described by Paraschivoiu [5]. A simpler model is found in Wilson [15].

In brief they model the flow as composed by a grid of linear streamtubes that partition the swept area. For each streamtube static airfoil data is used for calculation of the average blade forces from lift and drag using the relative velocity to calculate angle of attack and blade Reynolds number. This average force is then used to calculate the loss of momentum and; thus, the slowdown of the wind. These models cannot predict the structure of the wake but are on the other hand extremely fast and can with advantage be used for quick back-of-the-envelope calculations of blade forces and turbine performance. However, important effects like dynamic stall is not easily included. Also, because of the shortage of airfoil data for airfoils in flow with curved streamlines it is not obvious how to study this effect.

In terms of thoroughness the streamtube models has its opposite in the method of direct numerical simulation (DNS). In DNS the Navier-Stokes equations are solved numerically with some PDE-solving method such as Finite Elements (FEM) or Finite Volumes (FVM) using a fine enough grid to capture all relevant effects. To avoid simulating all fine scales an approach commonly used is the Reynolds Averaged Navier-Stokes Simulation (RANS) [5]. In theory this approach can be used for investigating all aspects of the turbine aerodynamics. However, the large computational cost associated with such methods limits the use to specific details such as aeroelasticity but even these applications come with a large computational cost. Various authors, notably Ferreira et al has used Large Eddy Simulations (LES) to study the effect of dynamic stall[16, 17]. The drawback is as with DNS that detailed simulations long enough for the wake to develop completely are very expensive. However, a recent article by Lida [18] shows promising development in this direction.

The fifth group of models, which will be built upon herein, are discrete vortex methods (DVM). The lift force of a blade in Darrieus motion is due to the circulation around the blade. However, because the lift and; thus, the circulation is changing during a revolution a continuous line of eddies are shed from each wing in order to conserve the angular momentum of the air. Discrete Vortex Methods discretize this line into separate eddies and track the resulting eddies as they are convected downstream.

The first vortex simulations of Darrieus turbines in inviscid flow was performed by Strickland [4]. It was based on the principle of a lifting line approach using airfoil data sheets to calculate the circulation.

Lifting line approaches has since then been the most popular of the vortex models. Another possibility are panel methods, which can simulate the behavior of completely general airfoils as long as the flow is attached, without the need for experimental data. However, they have computational disadvantages because they require an inversion of a large full matrix and subsequent multiplication at each time-step.

An early model based on a conformal mapping technique was proposed by Wilson [14]. The idea was

to map an airfoil to a circle. The flow could then be found using the method of images. Conformal mapping

techniques share the principal advantages of the panel methods if an appropriate mapping can be found,

at the same time conformal mappings are very efficient in terms of computational cost. An extension of

Wilson’s conformal mapping approach was made by Deglaire [19] where a general manner of determining

a suitable conformal map was found.

(11)

1.5 Previous Multibody Potential Flow Solutions

Milne-Thomson first stated and proved the circle theorem [20] that allows the introduction of a cylinder in any known potential flow. Williams [21] used the circle theorem iteratively to derive a solution to the flow around two circles. Resent work by Crowdy [22] demonstrated an analytical expression for potential flow around multiple cylinders. In this diploma work, however, an iterative method will be used.

Ives [23] used Laurent-series, which were obtained by the discrete Fourier transform, to find the poten- tial flow around a wide class of two dimensional shapes. Deglaire[19] also used discrete Fourier transform when extending Williams work to a flow around an arbitrary number of circles. Deglaire suggest the use of his technique for application to vertical axis wind turbines. Because most VAWTs rotate as solid bodies it is not necessary to remap at each time-step.

Another approach was described in a paper by Zannetti et al. [24]. They used sophisticated mathemat-

ical techniques to map two specific wings conformally to a simpler geometry. However, their approach is

limited to two blades. Also they appear to have no way to try different blade sections in a general manner.

(12)

Chapter 2

Aerodynamic Assumptions

In general air flow is governed by the Navier-Stokes equation.

ρ Dv

Dt = −∇p + ∇ · T + f (2.1)

where D/Dt ≡ ∂/∂t + v · ∇, T is a tensor describing the viscous properties of the fluid; v the velocity; p the pressure and ρ is the air density (1.225 kg/m 3 ).

There is little hope of finding analytical solutions to (2.1) in all but the simplest of geometries. Even for cases where an analytical solution exist they inevitably break down as the speed of the flow increases because the equations destabilize. In the case of a vertical axis wind turbine the complexity of the geometry and motion forces us to resort to numerical methods. However, solving the Navier-Stokes equations makes, even then, a very challenging task because of the high sensitivity of the solutions on the various simulation parameters such as step-length in the time and space discretization. With today’s computer resources, nu- merical solution in full 3D of (2.1) for a VAWT-geometry is unfeasible, even when using the most powerful supercomputers. In applications it is; thus, necessary to make approximations and simplifications and in order to do that, avaliable a priori knowledge about both the flow must be described and taken into account.

Wind turbines are operating within the earth surface boundary layer. Therefore the wind speed is sig- nificantly dependent of the distance to the ground i.e., wind shear. Further, depending on the roughness of the earth surface the flow will be more or less turbulent. Moreover still, some aspects of the flow around straight bladed giromills are known. A notable effect is the downwash effect at the tip of each blade. The formation of tip-vortices which decrease the lift and increase the drag. We thus must ask: Can we simplify, model or neglect:

1. Height dependent wind speed differences 2. Turbulence in the incoming flow?

3. Tip effects from the rotor blades?

4. Viscous effects such as skin friction drag on the struts and blades?

The average wind speed is increasing with height over ground. Because a wind turbine tower is typically 1-100 meter high, the velocity gradient at these heights indicate that there is a velocity difference between the bottom and the top parts of the rotor. Further, this difference is smaller the further from the ground the turbine is mounted. In the case of a very small turbine with a 10 meter tower and 5 meter blades the difference in undisturbed wind velocity might be as much as 20% between the bottom and top of the parts of the rotor. An important implication of this effect is that the tip-speed-ratio varies throughout the height of the turbine.

Although the average wind speed is indeed increasing with height over ground. It seems like the short time variance in wind speed, however, is roughly constant. This means that the impact of natural turbulence will be less on larger turbines [9]. That is not to say it is without importance.

Moreover, while the turbine operates it is absorbing kinetic energy from the air. Hence the air speed

downwind from the rotor (the wake) must be lower than the undisturbed wind speed. Because the air is free

streaming it can be assumed that a part of the air that was flowing in the projected frontal area of the rotor

(13)

in the end is flowing around the rotor where it is not producing power. See figure 2.1 for an illustration of the general properties of the expanding flow.

Figure 2.1: The behavior of the air streams around an H-rotor type wind turbine

The lift of an airfoil is due to the pressure difference between its two sides. Air flow follows pressure gradients; the tip of each blade thus provides the air opportunity to move from the lower to the upper side hence evening out this pressure somewhat. This results in two things. 1) There is a reduction of lift due to the lost pressure. 2) The air flowing over the tip of the blade results in a vortex line emanating from each blade-tip. These vortex lines induces an increase in drag which is related to the lift force.

Further, it is known that flow around a cylinder exhibits three dimensional eddies at Reynolds numbers greater than 10 4 and fully turbulent (a large amount of eddies in various sizes) at Reynolds numbers greater than 10 5 . It is thus resonable to assume that the flow is three-dimensional for airfoils in Darrieus motion in operating conditions (typically the blade Reynolds number is exceeding Re > 10 5 ).

2.1 Two-dimensional Approximation

There are indeed complicated three-dimensional effects present in the flow. A treatment of these effects are necessary in order not to overestimate the performances of a given turbine.

It is possible to estimate the losses due to the finite span of the blades through the lifting line theory of Prandtl [25]. But there is an important point to make before the standard corrections for finite span are applied to the results. For steady conditions the starting vortex of a wing is too far away to influence the flow. Thus there is no downwash from it felt by the wing and hence no drag due to lift. Thus the drag due to lift is completely dominated by the downwash from the wingtip vortices. Hence, only 3D wings have drag due to lift in steady flight. On the other hand, the influence from the starting vortex dominates the downwash for a wing that is in the process of building up its lift. And the downwash from a 2D starting vortex is stronger than for a finite length 3D starting vortex.

Given that the blade-tip losses can be estimated from knowledge of the performance of infinite span blades it makes sense to study the simplified case of infinite aspect ratio. Furthermore if the effect of the wind shear is to be studied a rudimentary 3D model can be accomplished by running multiple simu- lations for different height intervals. In addition, small three dimensional eddies can be modelled using Kolmogorov’s theory or other turbulence models [26]. Hence it will be assumed that the flow is two- dimensional; neglecting the influence of the three dimensional turbulence effects.

We thus want to find the two-dimensional flow which satisfies the boundary conditions of the undis-

turbed wind velocity far up-stream, and no-through and no-slip conditions at the rotor boundaries.

(14)

Although a three dimensional simulation would be preferable, there has not been time to develop and test such a tool. Also, the computational cost would be substantially higher.

2.2 Inviscid Approximation

Ludwig Prandtl discovered that given a Reynolds number sufficiently large, the flow over a surface can be divided into two distinct regions. There is one region away from the boundary surface where the fluid behaves as an ideal fluid. In the other region, which surrounds the boundary surface, viscosity is important.

This other region is called the boundary layer. Prandtl also pointed out that the thickness of the boundary layer tends to zero as the Reynolds number increases (see [27] for an interesting article explaining the boundary layer concept; its history and significance). For a VAWT the Reynolds number is relatively high (typically exceeding 10 5 ). It is thus reasonable to make the assumption that the boundary-layer thickness can be approximated to zero. This simplifies the problem at hand to solving the Euler equation with no- through boundary condition.

A problem with inviscid flow around a wing is that the flow is not uniquely defined. There can be an arbitrary circulating flow around the wing. In fact, this circulation determines the placement of the stagnation points in the case of steady flow.

It is further assumed that all shear in the flow takes place in the boundary layer or in a thin wake behind the lifting bodies. Thus the flow outside the boundary layer and not in the wake is irrotational.

2.3 Incompressible Approximation

A rule of thumb is that a fluid behaves in an incompressible manner when the Mach number is smaller than 0.3 [28]. This means that air-flow can be considered incompressible for speeds lower than 100 m/s. As it turns out this approximation applies to a majority of VAWTs in operating conditions.

In an incompressible inviscid fluid, disturbances travel at a very high speed (if the fluid is perfectly incompressible the speed of sound is infinite). Thus, for the flow problem to be well posed it is necessary to specify the flow at some initial time t = 0 as well as on the boundaries of the domain at all times.

Under these assumptions the fluid flow is governed by the Euler equation and the Continuity equation Dv

Dt = − 1

ρ ∇p (Euler equation) (2.2)

∇ · v = 0 (Continuity equation) (2.3)

where v is the velocity field; p the pressure field; ρ the density and D/Dt = ∂/∂t + v · ∇ the convective derivative.

2.4 Boundary Layer Considerations

The Euler equations poses no restrictions on the flow regarding the continuity of the velocity field. Such solutions are however unphysical because all real fluids have a certain viscosity although it might be small.

In particular at the sharp trailing edge of a moving wing the flow comes together from the upper and lower side of the wing. If the stagnation point deviates from the trailing edge the resulting sharp velocity shear gives rise to the formation of a vortex inside the boundary layer. Indeed as the rotor starts up a continuous sheet of vorticity is shed from the trailing edge of each wing. Figure 2.2 and figure 2.3 illustrates the effect.

The angular momentum of the air is conserved through time. The reaction of the angular momentum of the vortex sheet must thus be contained in a counter-circulation around each wing. This circulation is indeed giving rise to the lift force. The variation of angle of attack throughout the Darrieus motion implies a corresponding variation in circulation and that a turbine blade leaves a continuous line of vorticity behind its trailing edge. This line is called the wake of the wing.

The demand continuous velocity around the trailing edge is called the Kutta condition after the German

mathematician and aerodynamicist Martin Wilhelm Kutta [29]. Herein we will call any point that is required

a priori to be a stagnation point a point with Kutta condition.

(15)

(a) zero circulation (b) zero circulation

(c) Kutta condition (d) Kutta condition

Figure 2.2: Streamlines and Velocity field at the trailing edge of a 2D wing at an angle of attack.

−0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35

0.6 0.65 0.7

Figure 2.3: Simulated vortex sheet with starting vortex. The upper image depicts a complete boundary layer vortex

simulation at Re=8000 (by courtesy of Paul Deglaire) whereas the lower image shows the approximation of an infinitly

thin boundary layer (Re→ ∞) used in the present model.

(16)

Chapter 3

Flow Around a Moving Rotor

A shear layer is is formed from the trailing edge of each wing, as explained in section 2.4. This sheet will be modeled as a discrete chain of distinct eddies (called vortices, sing. vortex). The dynamics of such vortices can be determined from the principal equations of fluid dynamics.

3.1 Vortex Dynamics

In nature a two-dimensional vortex consists of two regions that exhibits different behavior. In the vortex center there is a rotational core that rotates as a solid body. Outside the core the flow is almost irrotational.

The velocity field from a single vortex placed at the origin can; accordingly, be accurately modeled [26] as

v(x) = ( γ

2π ˆz × x if |x| ≤ R

γ 2π

ˆz×x

|x|

2

if |x| > R. (3.1)

where R is the core radius and γ is the circulation around the vortex. This type of vortex is often called a Rankine vortex. In the special case of R = 0 the vortex shall be called a point vortex (in 3D it is called a line vortex). A smoother vortex model is the so called Lamb-Oseen vortex, which is modeled as

v(x) = γ 2π



1 − e −(

|x|R

)

2

 ˆz × x

|x| 2

= γ 2π



1 − e −(

|x|R

)

2



∇log(|ˆz × x|) (3.2)

It is important to note that v(x) is essentially conservative outside the vortex core (i.e when |x|  R). Plots of the velocity profile of the different vortex modeles are displayed in figure 3.1

0 1 2 3 4 5 6

0 0.2 0.4 0.6 0.8 1 1.2

Radial distance to the vortex center

Velocity

Rankine vortex Oseen vortex Point vortex

Figure 3.1: Velocity profiles of the most common vortex models

It is possible to determine rules for the kinematics of vortices from Euler’s equation Dv

Dt = − 1

ρ ∇p (3.3)

(17)

where v is the velocity; p the pressure and ρ is the density.

Taking the curl of (3.3)

D(∇ × v)

Dt = 0

Vorticity is defined as ω = ∇ × v, which gives the vorticity equation of incompressible flow Dω

Dt = 0 (3.4)

Given the vorticity field (i.e., a solution to (3.4)) it is possible to express the velocity field by use of the Biot-Savart relation, in our case for two dimensions [30]

v(x) = ∇φ 0 − 1 2π

ZZ

ω(x 0 ) × (x 0 − x)

|x − x 0 | 2 dΩ (3.5)

where φ 0 is the irrotational velocity potential generated by the boundary conditions.

The trick of vortex particle methods is to temporarily assume that all vorticity is concentrated around a countable set of points {x k } called vortex particles. The vorticity field can then be represented as

ω(x) = ∑

k

γ k δ(x − x k )ˆz (3.6)

where δ(x) is the Dirac distribution. The velocity field can be calculated from (3.5) and (3.6) to be v(x) = ∇φ 0 (x) − ∑

k

γ k

2π ∇log |ˆz × (x − x k )|. (3.7)

To avoid the singularity and obtain a more physical velocity field the point vortices are now replaced with Lamb-Oseen vortices with a small core radius R. Thus, we have

v(x) = ∇φ 0 (x) − ∑

k

1 − e



|x−xk|

R



2

! γ k

2π ∇log |ˆz × (x − x k )| (3.8) Clearly (3.8) and (3.7) coincide almost identically outside the vortex cores.

Now, Helmholtz proved (in a seminal 1858 paper [31]) that a fluid element initially located at the center of a vortex will continue to stay at the center for all time. Another way of saying that is that a vortex move with the same velocity as the fluid does.

This means that the movement of vorticity is governed by a system of ordinary differential equations dx k

dt = ∇φ 0 (x) x

k

− ∑

j6=k

1 − e



|xk−xj|

R



2

 γ j

2π ∇log |ˆz × (x − x j )|

x

k

(3.9)

3.2 Potential Flow

Outside the cores of the vortices it holds

v(x) = ∇φ 0 (x) − ∑

k

γ k

2π ∇log |ˆz × (x − x k )|.

Factoring out the ∇ from the previous expression

v(x) = ∇ φ 0 (x) − ∑

k

γ k

2π log |ˆz × (x − x k )|

!

def ≡ ∇φ

The continuity equation (2.3) for incompressible fluids reads

∇ · v = 0

(18)

This implies that ∇φ(x) satisfies the Laplace equation

∇ · ∇φ = ∆φ = 0 (3.10)

i.e φ is a harmonic function.

It is beneficent to construct an analytic function from φ(x, y) and its harmonic conjugate ψ(x, y). The well known Complex Flow Potential is obtained, which in this case is induced by the boundaries of Ω.

w(z) = w(x + iy) = φ(x, y) + iψ(x, y) (3.11)

ψ is called the Stream function and is of great importance for fluid dynamics [20]. Physically ψ(z) mea- sures the fluid flux through the line connecting z and a reference point z 0 . In the following, the common association between a vector x ∈ R 2 and the complex number z ≡ x + iy ∈ C will be used. Then (3.11) can be written

w(z) = w ext (z) − ∑

k

γ k

2π log(z − z k ) Equations (3.8) and(3.9) can now be written concisely as

v(z) = dw

dz (3.12)

dz k dt = dw

dz z

k

(3.13)

3.3 Boundary Conditions

The potential flow can be described as superposition of four types of flow. The first part is the wind; the second is the headwind due to the rotation of the rotor; the third is the flow due to free vortices in the flow;

and the fourth and last is the action from the boundaries of the turbine.

It is assumed that there exists a complex flow potential that governs the flow as described in section 3.2 w(z) = w (z) + w γ (z) + ˆ w γ (z) + η(z). (3.14) Here w (z) denotes the flow of the undisturbed wind; w γ (z) denotes the flow induced by free vortices; and η(z) denotes the disturbance induced by the turbine. The function ˆ w γ is the complex flow potential of the special vortices that are currently forming in the imagined boundary layers of the blades as described in section 2.4.

The complex potential η(s) in turn, can be written as

η(z) =

N

n=1

η n (z).

where η n (z) denotes the disturbance of each turbine part ∂ T n . Equipped with this notation the boundary conditions can be determined for the parts individually. The far field is only effected by the wind. Hence the boundary conditions become

dw ∞

dz → V dw γ

dz → 0 d w ˆ

dz → 0 dη

dz → 0 as |z| → ∞ (3.15)

The physical reality is that no fluid can penetrate the surface of the rotor. I.e., the no-through boundary condition. This means that the wings of the turbine must be streamlines. When the turbine is in operation the rotor moves in rigid body rotation. Assuming that β denotes the azimuthal angle measured to a reference blade of the turbine then the angular velocity of the rotor is defined by

β = ˙ 2λV

D

where λ is the tip-speed-ratio and D the turbine diameter.

(19)

Accordingly, the turbine wings must be streamlines in the headwind −i ˙βz. This corresponds to a stream- function ˙β|z| 2 /2. The potentials w ∞ (z) and w γ (z) are induced by the wind and the free vortices; they are free to take any value on ∂ T . It is; therefore, necessary that the complex flow potential η(z) satisfy

ℑ(η(z) + w ∞ (z) + w γ (z) + ˆ w γ (z)) + 1 2

β|z| ˙ 2 = C n z ∈ ∂ T n (3.16) where C n is a real constant that may be different on different wings. The notation ∂ T n has been introduced to mean the boundary of blade n. Similarly ∂ T will denote the boundary of the whole rotor i.e., the union of the blades and Ω denotes the fluid region.

In applications it is known that certain points {z te,k } on the boundary are stagnation points (in particular the trailing edges of the blades). These relations will be called Kutta conditions even though they mean a slightly more general thing than what is normally called Kutta condition.

v(z te,k ) = dw

dz (z te,k ) = i ˙βz te,k z te,k ∈ ∂ T (3.17)

3.4 Conformal Mapping

It is a property of certain functions called conformal maps that a solution φ(x, y) to the Laplace equation (3.10) can be transformed to a different geometry, using the conformal map, yet still solve the Laplace’s equation in the new coordinates.

This is of interest if a conformal map between the area outside the rotor and a simpler geometry can be found.

Specifically, if a function φ(x, y) that solves Laplace equation is transformed with a conformal map f (u + iv) = x(u, v) + iy(u, v) to the new function φ 0 (u, v) def ≡ φ(x(u, v), y(u, v)) then

2 φ 0

∂u 2 + ∂ 2 φ 0

∂v 2 = 0

In a recent paper Deglaire proposed a method for mapping the area outside N distinct simply connected subsets of C to the area outside N distinct discs conformally [19].

Figure 3.2: The idea of using conformal mappings to simplify the problem

Deglaire suggests that the conformal mapping can be written as a kind of Laurent-series

z = f (s) = s +

N n=1 ∑

k=1

c n k b k n

(s − s c n ) k ≈ s +

N n=1 ∑

N

0

k=1

c n k b k n

(s − s c n ) k . (3.18) where s c n are the circle centers and b n the circle radii. He then develops an algorithm to determine the complex coefficients {c nk }. This technique provides a straight forward manner to find maps as in figure 3.2 for an extended family of rotors. The task is; thus, changed to the simpler problem of finding φ for the flow around N cylinders.

Accordingly, let z = f (s) be a function that maps conformally, the area outside N circles s : |s − s c n | = b n ,

in what shall be denoted the circle plane, to the area outside N airfoils, denoted the physical plane.

(20)

To obtain the velocity in the physical plane the fact that the complex flow potential is invariant under conformal transformations can be used (see for example [32]). Thus, the velocity carries over to the physical plane by the simple formula

v(s) = v s (s) 1

f 0 (s) . (3.19)

here a note of caution is necessary. Strictly speaking (3.19) is only valid if v s (s) is irrotational at s. Strictly this is not the case when there are vortices in the flow with core radius R > 0. Nonetheless (3.19) will be used with the restriction that the core radius must be kept small.

Because the two functions φ(x, y) and ψ(x, y) both solves Laplace’s equation it is easy to varify that the complex flow potential w also is invariant under a conformal transformation. From equation (3.14):

w(s) = w (s) + w γ (s) + ˆ w γ (s) + η(s). (3.20)

3.5 Boundary conditions in Circle plane

We now want to express the boundary conditions (3.15) and (3.16) in the simpler circle plane. The wind correspond to a boundary condition at infinity. The conformal map f , (3.18), which relate the physical plane to the circle plane, has the property

f 0 (s) → 1, as |s| → ∞

This implies that the boundary condition (3.15) — which describe the behaviour as |s| → ∞ — can be expressed independently of f as

dw

ds → V dw γ

ds → 0 d w ˆ γ

ds → 0 dη

ds → 0 as |s| → ∞ (3.21)

The boundary condition at the turbine surface is slightly more complicated. That is due to the rotation of the rotor. This means that the shape of the wings matter because the distance to the center of rotation defines the boundary condition of a point. The boundary condition (3.16) takes the form

ℑ(η(s) + w ∞ (s) + w γ (s) + ˆ w γ (s)) + 1 2

β| f (s)| ˙ 2 = C n s ∈ ∂ T sn (3.22) In the circle plane the Kutta condition can be written as

v s (s te,k ) = i ˙β f (s te,k ) f 0 (s)

Because the velocity is otherwise bound to follow the circle it is enough to require

v s (s te,k ) − i ˙β f (s te,k ) f 0 (s) ⊥ `(s te,k ) (3.23) where s te,k denotes a point on one of the circles z te,k = f (s te,k ) and `(s) denotes a tangent vector to the circle.

3.6 Construction of Complex Flow Potential

From the boundary condition (3.21) and section 3.1 the following forms of the partial complex potentials of w in equation 3.20 are proposed

w (s) = V s w γ (s) = ∑

k

k

2π log(s − s vk ) w ˆ γ (s) = ∑

n

iˆγ n

2π log(s − ˆ s n ) η(s) =

N n=1 ∑

η n (s)

(21)

the wind velocity V ; the vortices {(γ k , s vk )} and the kutta vortex release positions ˆ s n are all assumed to be known. The goal of this section is to present an algorithm to determine η(s) and {ˆγ n } such that the no-through and the Kutta conditions are satified. In order to describe the algorithm a notation for partial results under iteration step j is used as in the following

w ˆ ( j) (s) γ = iˆγ ( j) n

2π log(s − ˆ s n ) η ˆ ( j) n (s) = − iˆγ ( j) n

2π log(s − ∆ n ( ˆ s n )) η n (s) = ∑

j=1



η ( j) n (s) + ˆ η ( j) n (s) 

where a symbol ∆ n (s) = s−s b

2n

cn

+ s cn for the point symmetric with respect to circle n has been introduced.

The condition (3.22) can be reformulated as ℑ(η n (s)) = −ℑ(w (s) + w γ (s) + ˆ w γ (s) + ∑

m6=n

η m (s)) − 1 2

β| f (s)| ˙ 2 +C n s ∈ ∂ T n

The approach will be to first find the potential as if the blades has no influence on each other (i.e.,

∑ m6=n η m (s) is neglected for the moment). Also, for the moment the Kutta condition is ignored. Hence determine η (0) n (s) such that

ℑ(η (0) n (s)) = −ℑ(w ∞ (s) + w γ (s)) − 1 2

β| f (s)| ˙ 2 +C n s ∈ ∂ T n .

This is possible because all entities are known. Then a first approximation {ˆγ 0 n } can be determined such that d w ˆ (0) γ

ds s

te,n

+ d ˆ η (0) n ds

s

te,n

≡ iˆγ (0) n

"

1

s te,n − ˆs n − 1 s te,n − ∆ n ( ˆ s n )

#

= − dη (0) ds

s

te,n

− dw ds

s

te,n

+ dw γ ds

s

te,n

+ i ˙β f (s te,n ) f 0 (s) s

te,n

The iteration step is to determine η ( j+1) n (s) such that ℑ(η ( j+1) n (s)) = −ℑ( ∑

m6=n



η ( j) m (s) + ˆ η ( j) m (s) 

+ ˆ w ( j) γ (s)) +C n

and {ˆγ n j+1 } such that d w ˆ ( j+1) γ

ds s

te,n

+ d ˆ η ( j+1) n

ds s

te,n

≡ iˆγ ( j+1) n

"

1

s te,n − ˆs n − 1 s te,n − ∆ n ( ˆ s n )

#

= − ∑

m6=n

( j+1) m ds

s

te,n

It shall be assumed that η n (s) can be written as a series

η n (s) =

k=1

a nk b k n (s − s cn ) k

furthermore, is is also assumed that the functions η ( j) n (s) have the same form

η ( j) n (s) =

k=1

a ( j) nk b k n (s − s cn ) k

Under this assumption η is analytic on Ω and the boundary condition (3.21) is automatically fullfilled. To satisfy (3.22) ∂ T sn is parameterized by s cn + b n e . On ∂ T sn it; thus, holds that

η ( j) n (s cn + b n e ) =

a ( j) nk e −ikθ

(22)

This expression can be identified as negative-frequency part of a Fourier series. If the series is truncated after N co terms and then evaluated at 2N co equally spaced points on the circles. Then the coefficients can be obtained by use of the inverse Discrete Fourier Transform ( F −1 )

{a (0) nk } k = F −1 {−iℑ(w ∞ (s np ) + w γ (s np )) − i 2

β| f (s ˙ np )| 2 } p {a ( j+1) nk } k = F −1 {−iℑ( ∑

m6=n



η ( j) m (s np ) + ˆ η ( j) m (s np ) 

+ ˆ w ( j) γ (s np )} p

where s np = s cn + b n e iπp/N

co

In practice, however, if there are vortices very close to circle n the number of terms N co becomes very high. The performance of the method is improved if the close vortices is handled separatelly using the method of images.

For each disc, one can divide w γ (s) in one part coming from vortices close to the disc and another coming from vortices far from the disc. Empirically, the best cut-off distance has been determined to be 0.1b n from the boundary of circle n [19].

w γ (s) = ∑

s

vk

k: ∈N

n

k

2π log(s − s vk ) + ∑

k:

s

vk

∈N /

n

k

2π log(s − s vk ) N n def ≡ {s : |s cn − s| ≤ 1.1b n }

now, η (0) n (s) is written

η (0) n (s cn + b n e ) =

k=1

a (0)

0

nk e −ikθ + ∑

s

vk

k: ∈N

n

k

2π log(s − ∆ n (s vk ))

where the coefficients {a (0)

0

nk } k are calculated by {a (0) nk

0

} k = F −1 {−iℑ(w ∞ (s np ) + ∑

k:

s

vk

∈N /

n

k

2π log(s np − s vk ) + i 2

β| f (s ˙ np )| 2 )} p

Figure 3.3 and 3.4 displays the convergence of the algorithm for the case of steady flow around a Ropatec wind turbine.

Figure 3.3: Streamlines of steady flow around a standing Ropatec rotor in the transformed and physical plane respec-

tively. Note there is no Kutta condition applied to the central body.

(23)

0 2 4 6 8 10 10 −8

10 −6 10 −4 10 −2 10 0

Iterations (M)

Error (max flux trough boundary)

Figure 3.4: The convergence of the algorithm described in section 3.6

3.7 Convection of Vortex Particles

The movement of the vortices are governed by equation (3.13) which reads dz k

dt = dw dz z

k

(3.24) with the current notation this can be written

dz k

dt = V + dη dz z

k

+ ∑

j6=k

1 − e



|zk−z j|

R



2

 iγ k

1

z k − z j (3.25)

which is a system of N v ordinary differential equations. These equations give the velocity of each vortex particle and; accodingly, contain the information needed to advance the vortices. Simply summing (3.24) for each vortex—called Direct Summation or sometimes Na¨ıve Summation—of course yield the velocities.

It is called na¨ıve because the number of basic operations in the computation scales as O (N v 2 ). A respectable work given that it must be computed in every time-step and for perhaps 100000 vortices. More efficient methods was developed in the eighties. These exploit the fact that a cluster of vortices itself behave as a large vortex placed at the center of vorticity (defined analogously to the center-of-mass concept of particle mechanics). A cluster of vortices; thus, moves around this point while at the same time moving somewhat as shown in figure 3.5.

This self-similar property can be efficiently exploited if recursively applied on single vortices, clusters

of vortices and clusters of clusters of vortices and so fourth. During the eighties this approach spawned

a number of efficient O (N v log N v ) divide-and-conquer algorithms that has been applied successfully. The

approaches mainly differ in how they determine the clusters, a common method beeing the Particle-in-Cell

algorithm [30], which relies on superimposion of a grid over the plane. In this method a cell of the grid

is thought of as defining a cluster. Another method was proposed by Barnes and Hut [33] in 1989. In

their method the box containing all vortices is recursively divided until each box only contains one or zero

vortices. In this manner a tree is built where each branch is defining a cluster. The idea is that the branches

far from each other only interact through their centers of vorticity. The most sophisticated algorithm to

date is the Fast Multipole Method (FMM) of Greengard and Rokhlin [34]. FMM does not simply replace

a cluster far away with a replacement vortex at the center of vorticity. The cluster is instead represented

(24)

Figure 3.5: Qualitative picture of the hierarchical motion of vortex clusters. The center-of-mass symbol represents the center of vorticity.

by a truncated Taylor series. FMM has the advantage that the error can be kept arbitrarily small. However, FMM relies on the potential field representation of the velocity. The error from this is however not believed to be large.

In this project the fast multipole method has been used for finding the velocities of the vortices. The author of this report has had access to a very efficient implementation by Stefan Engblom of the Division of Scientific Computing Department at Uppsala University. Further explanation can be found in [35].

Once the velocities of the vortices are known; they can be advanced in accordance with some numerical integration method. In the continuous case ∆t → 0 it is impossible for a vortex to be convected into a solid region. This is however not true when ∆t is discrete. It is therefore necessary to perform a check at each time-step whether a vortex has moved into a solid body. Two actions are possible when a vortex is found in a forbidden zone. 1.) The vortex could be removed from the flow or 2.) the vortex can be moved to the closest boundary. In the case of the former alternative the boundary layer can be said to have absorbed the vortex. The reaction must; accordingly, be a change in circulation around the corresponding wing. In this project the latter option has been chosen for the simplicity of implementation.

3.8 Routh’s Rule

The convection of vortex particles is complicated by the use of conformal mappings to transform the do- main. A condition for the relation between velocities in the two planes (3.19) to be valid is that v s is irrotational at the point s. This condition is satisfied everywhere except on the vortex particles themselves.

The velocity of a vortex is determined by the velocity induced at it’s center by all other causes except the vortex itself. Thus, to derive the velocity in the physical plane it is logical to start with the complex flow potential in the physical plane. It is; hence, assumed that a vortex is located at z v = f (s v ). Then the invariance of the complex flow potential, w, can be written as

w(s) ≡ w s

v

(s) + iγ v

2π log(s − s v ) = w z

v

(z) + iγ v

2π log(z − z v ) ≡ w(z)

where w z

v

(z) and w s

v

(s) denotes the complex flow potential induced from all other causes except the vortex with center at z v and s v respectively. From this expression Clements [36] used a Taylor series expansion to show that the velocity of the vortex can be written as

v z

v

= dw z

v

dz

z

v

= dw s

v

ds

s

v

ds dz z

v

+ iγ v 2π 1 2

d f

−1

dz

z

v

d f

−1

dz

z

v

(3.26)

Using the implicit function theorem (3.26) is rewritten as

v z

v

= dF z

v

dz

z

v

= dF s

v

ds

s

v

− iγ v 2π 1 2

f 00 (s v ) f 0 (s v )

! 1

f 0 (s v ) (3.27)

However, the algorithm for satisfying the boundary conditions is using the circle plane. Therefore it

(25)

would be convenient if the velocity of the vortices in this plane also could be found. We write v s = lim

h→0

f −1 (z v + v z h) − f −1 (z v ) h

= lim

h→0

f −1 (z v ) + v z h d f dz

−1

z

v

+ 1 2 v 2 z h 2 d

2

f

−1

dz

2

z

v

+ O (h 3 ) − f −1 (z v ) h

= v z d f −1 dz

z

v

= v z 1

f 0 (s v ) (3.28)

Finally, by inserting (3.27) in (3.28) the velocity of the vortex particles, expressed in the circle place is obtained

v s = dF s

v

ds

s

v

− iγ v 2π 1 2

f 00 (s v ) f 0 (s v )

! 1

| f 0 (s v )| 2 (3.29)

3.9 Convection in Circle Plane

Equation (3.28) gave the relationship between how a vortex move in the physical plane and how the corre- sponding vortex moves in the circle plane. It is obvious that the two expressions

z v,k (t) = Z t

0

v z (z v,k (t 0 ))dt 0 (3.30)

s v,k (t) = Z t

0

v s (s v,k (t 0 ))dt 0 (3.31)

are equivalent in the sense that

f (s v,k (t)) = z v,k (t) ∀t > 0

However, this is not necessarily true for a numerical integration scheme. By comparing the predicted motion of a single vortex around an airfoil it was discovered that (3.31) gave predictions of higher quality for the same time-step. An illustration of this interesting fact can be viewed in figure (3.6)

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1

−0.6

−0.4

−0.2 0 0.2 0.4 0.6

0.8 Physical plane ∆t=0.01

Circle plane ∆t=0.1 Physical plane ∆t=0.1 Panel method ∆t=0.01

Figure 3.6: Single vortex moving in counter-clockwise path around a Joukowski airfoil using the midpoint integration method and three different ways of calculating the velocity

The experiments visualized in figure 3.6 indicate that (3.31) can be used to calculate the paths of the

vortices instead of (3.30) without disadvantage. In this particular case it is interesting to note that (3.31)

even performs better than (3.30) at the specified time-step. Why this is so, and whether or not this is a

general property remains to investigate.

(26)

Chapter 4

Numerical Validation of the Model

In this chapter the validity of the flow model is examined by comparing computational results with values obtained by experiments and other models.

4.1 Steady flow

The steady flow around a NACA0018 airfoil has been computed, as a first means for validation of the model.

For this simple case it is sufficient to compare the results with those of the well known XFoil code with viscous modeling enabled[37]. The most suitable parameter for comparison is the pressure coefficient that is defined as

C p = p − p

1 2 ρV 2

where p V ∞ and ρ ∞ are the undisturbed pressure, velocity and density respectively. The results are shown in figure 4.1. The agreement is excellent.

−0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5

−1

−0.5 0 0.5 1 1.5 2 2.5

−Cp

chord

Present model Xfoil

Figure 4.1: Streamlines and pressure coefficient over a NACA0018 airfoil. Angle of attack 8 degrees. Re = 6 · 10 5

§11

4.2 Blade Loads

The operational loads on the turbine blades has been computed as a means for the validation of the complete numerical model.

The forces can be calculated numerically from integration of the unsteady Bernoulli equation.

p = −ρ  |v| 2 2 + ∂φ

∂t



+ p 0 (4.1)

(27)

where v denotes the velocity relative to the airfoil and ∂φ

∂t is the time derivative of the velocity potential.

When the chord-radius ratio is small it is believed that the ∂φ

∂t term can be neglected. For example, in the case of a single blade at tip-speed-ratio five and c/R = 0.25 the error of this approximation is small (see figure 4.2). None-the-less this can be improved in a future version of the code.

0 100 200 300 400

−35

−30

−25

−20

−15

−10

−5 0 5 10 15

Azimutal angle (deg) − 1080

Non−dimensional normal force Exact formula

Unsteady Bernoulli Bernoulli relative velocity

(a)

0 100 200 300 400

−0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5

Azimutal angle (deg) − 1080

Non−dimensional tangential force

Exact formula Unsteady Bernoulli Bernoulli relative velocity

(b)

Figure 4.2: Comparison of blade loads calculated with an exact formula; unsteady Bernoulli principle and static Bernoulli principle

We thus calculate the pressure force on a wing by

F n = −ρ Z

∂ T

n

|v| 2 2 ˆn ds

By convention this is then divided in tangential and normal parts using a cylindrical coordinate system.

F N,n = F n · ˆr (4.2)

F T,n = F n · ˆ θ (4.3)

the force F T is providing the torque that produces power. The normal force is important for the study of

the structural loads on the turbine. For ease of comparison the forces are then nondimensionalized by the

References

Related documents

Object A is an example of how designing for effort in everyday products can create space to design for an stimulating environment, both in action and understanding, in an engaging and

In conclusion, it can be understood that there are elements of the State acting both as an investor and as a public authority when a public tender is conducted.

Figure 11: Simulated delivered power by each turbine for a system with three wind turbines with constant wind speeds and a voltage source load on the DC-bus.. The value of the power

In the end of 2009 a new kind of vertical axis wind turbine, a giromill 3 blades type, has been built in Falkenberg, by the Swedish company VerticalWind2. The tower of this

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

The accuracy of the dynamic stall model decreases with increased angle of attack, which is shown in [25,29], where the dynamic stall model was tested against wind tunnel data for

Optimization of cycloidal water turbine and the performance improvement by individual blade control. In Seong Hwang, Yun Han Lee, Seung Jo Kim. Performance investigation of

The goal in this project was to design a 20kW, PM, outer rotor-type generator for a vertical axis wind turbine. A number of generators has been designed and simulated with the