Performance analysis of wind turbines at low tip-speed ratio using the Betz-Goldstein model

Energy Conversion and Management 126 (2016) 662–672

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Performance analysis of wind turbines at low tip-speed ratio using the Betz-Goldstein model Jerson R.P. Vaz a,⇑, David H. Wood b a b

Faculty of Mechanical Engineering, Institute of Technology, Federal University of Pará, Av. Augusto Correa, N 1 Belém, PA 66075-900, Brazil Department of Mechanical and Manufacturing Engineering, Schulich School of Engineering, University of Calgary, Calgary T2N 1N4, Alberta, Canada

a r t i c l e

i n f o

Article history: Received 12 June 2016 Received in revised form 23 July 2016 Accepted 15 August 2016

Keywords: Wind turbine performance Swirl Low tip-speed ratio General momentum theory Tip-loss factor

a b s t r a c t Analyzing wind turbine performance at low tip-speed ratio is challenging due to the relatively high level of swirl in the wake. This work presents a new approach to wind turbine analysis including swirl for any tip-speed ratio. The methodology uses the induced velocity field from vortex theory in the general momentum theory, in the form of the turbine thrust and torque equations. Using the constant bound circulation model of Joukowsky, the swirl velocity becomes infinite on the wake centreline even at high tipspeed ratio. Rankine, Vatistas and Delery vortices were used to regularize the Joukowsky model near the centreline. The new formulation prevents the power coefficient from exceeding the Betz-Joukowsky limit. An alternative calculation, based on the varying circulation for Betz-Goldstein optimized rotors is shown to have the best general behavior. Prandtl’s approximation for the tip loss and a recent alternative were employed to account for the effects of a finite number of blades. The Betz-Goldstein model appears to be the only one resistant to vortex breakdown immediately behind the rotor for an infinite number of blades. Furthermore, the dependence of the induced velocity on radius in the Betz-Goldstein model allows the power coefficient to remain below Betz-Joukowsky limit which does not occur for the Joukowsky model at low tip-speed ratio. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction It is well-known that mathematical models for horizontal-axis wind turbines (HAWTs) design have increased in the last century [1,2]. Most models are based on the general momentum theory with some simplifying hypothesis for the complex flow in the wake [3,4]. As the tip-speed ratio X ¼ XR=V 0 decreases, the relative magnitude of the swirl velocity increases, which makes this region of turbine operation particularly challenging to analyze. In the definition of X; X and R are the angular velocity and tip radius of the blades, respectively, and V 0 is the wind speed. There have been a number of analyses of HAWT operating at low X, and least two, Wood [5] and Rio Vaz et al. [6], have suggested that the power coefficient may exceed the BetzJoukowsky limit of 59.3% as X # 0. This behavior is associated with the modeling of the trailing vorticity. The swirl, or circumferential velocity, represents the angular momentum generated by the blades and hence is directly related to the torque and power extracted by the rotor. The commonly used ‘‘free-vortex” assumption for the wake is equivalent to the Joukowsky model of the ⇑ Corresponding author. E-mail address: [email protected] (J.R.P. Vaz). http://dx.doi.org/10.1016/j.enconman.2016.08.030 0196-8904/Ó 2016 Elsevier Ltd. All rights reserved.

blades in which the bound circulation is constant. This causes an infinite velocity along the axis of rotation. Glauert [7] avoided this problem by terminating the integration of the blade element equations for momentum and angular momentum at the radius where the angular velocity of the swirl, x, equals X. Wilson and Lissaman [8] split the wake using X=xmax . As for a Rankine vortex, the swirl velocity near the axis is now proportional to the radius, r, and no infinities occur. The strong influence of swirl on the angular momentum balance for a rotor carries over to the energy and thrust equations. The power coefficient must be zero at zero X and it is considered as axiomatic that the power coefficient, C P asymptotes to the BetzJoukowsky limit as X becomes large. All references to maximizing performance in this paper will be to maximizing C P . Glauert [7] determined its monotonic increase with X to approach the BetzJoukowsky limit for an ideal rotor with an infinite number of blades, B. The simple analysis that leads to the Betz-Joukowsky limit ignores swirl. This can be justified on the grounds that HAWTs at high X produce power by the multiplication of a small torque and a high angular velocity, but the accuracy of this approximation has never been made clear. Wood [5] proposed a more general model for the hub vortex, and suggests that some account of the vortex structure of the wake will be required to resolve fully

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

663

Nomenclature a; a0 ab ; a0b A; A1 b B CP CQ CT dC P dC Q dC T dQ dT dr dr1 F FP p p0 p1 q

axial and tangential induction factors at the disc average axial and tangential induction factors at the rotor area of the disc and the cross section far in the wake axial induction factor in the wake number of blades power coefficient torque coefficient thrust coefficient local power coefficient local torque coefficient local thrust coefficient local torque (N m) local thrust (N) radial increment at the rotor (m) radial increment in the far wake (m) tip-loss factor Prandtl tip-loss factor dimensionless pitch for the helical vortex pressure in the free stream (Pa) pressure in the wake (Pa) dimensionless circulation

the effects of swirl. Based on the analytical solution for the velocities induced by helical vortex filament, Okulov and Sørensen [9] analyzed in detail the original formulation of the Goldstein function [10] for optimum rotor loading. They found that the model can be extended to heavily loaded rotors in full accordance with the general momentum theory at all X. Sørensen and Kuik [11] proposed a refined model based on the momentum including the friction and pressure on the lateral boundary of the annular streamtube intersecting each blade element. In their model, the lateral pressure is proportional to the area expansion and the pressure drop over the rotor at the edge of the disc, Dpjr¼R ðA1  AÞ. However, the term  is not fully understood because it is not known a priori how to determine such a term, and the authors are unaware of reference with further details on . Rio Vaz et al. [6] present an optimization model taking into account the influence of the swirl, based in the formulation proposed by Wilson and Lissaman [8] for the power coefficient. They found that the efficiency of a HAWT is heavily influenced by wake rotation at small X. Mikkelsen et al. [12] used a CFD actuator line model to analyze the wake at low X. It was shown that the excessive swirl near the centreline at low X generates vortex breakdown, causing a recirculation zone in the wake that limits the power yield of the rotor. The appearance of vortex breakdown has a similar effect on the flow behavior as the vortex ring state that usually appears at higher X. Wood’s [13] lifting line analysis of optimal loading at very low X rediscovered Goldstein’s [14] result for the asymptotic behavior of his function. Wood [13] found that the induced velocity is linear in radius for any number of blades. This inviscid analysis avoided the infinite velocities that usually result from the assumption of inviscid flow. Wood [15] considered the limiting case of the Goldstein and Glauert analyses as X # 0 for an ideal rotor with an infinite number of blades, B. He showed that the former is physically untenable as it requires significant energy extraction and a constant induced velocity for a stationary rotor. A modification of Goldstein’s analysis produced a more satisfactory physical model in which the swirl was linear in r near the axis and inversely proportional to r at higher values of the local-speed ratio, X r . The

r r1 rc S u u1 uh uh1 uh;max V0 X w

radial position at the disc (m) radial position in the wake (m) vortex core radius (m) swirl number axial velocity at the disc (m/s) axial velocity in the far wake (m/s) swirl velocity in the near wake (m/s) swirl velocity in the far wake (m/s) maximum swirl velocity in the wake (m/s) freestream velocity (m/s) tip-speed ratio w=V 0 , normalized velocity of the vortex sheet

Greek Symbols exponent in Vatistas vortex model n r=r c q density of air ðkg=m3 Þ / flow angle (rad) x angular speed of the flow in the wake (rad/s) X angular speed of the disc (rad/s) C vortex circulation ðm2 =sÞ

a

analysis, however, was limited to infinite B without wake expansion. The analysis for finite B and wake expansion using an induced velocity field dependent on r are the main topics of this paper. A detailed analysis is developed for the swirl influence on the thrust, torque and power coefficients that is valid for all X. General analytical formulations are developed for the azimuthal velocity distributions of the Joukowsky [16] model combined with the Rankine, Vatistas [17] and Delery [18] models for the vortex core. The effect of finite B on the general theory is accounted for by using two formulations for the tip-loss factor: that of Prandtl [19] and the recent method of Wood et al. [20]. The results show that a swirl velocity dependent on r produces an acceptable physical behavior. The generalized formulations for the thrust, torque and power coefficients, are in good agreement with previous results. The torque coefficient has a maximum value of about 0.5 when B tends to infinity, which agrees with Wood [15]. The remainder of this paper is organized as follows. The next section introduces the general momentum theory, showing the expressions for the performance coefficients using the Joukowsky model with Rankine, Vatistas and Delery vortices. Section 2 describes the general form of the Betz-Goldstein (BG) model. Section 3 shows the results and discussion, where the dependence on radius of the induced velocity field is presented, as well as the integral form of the performance coefficients taking into account the tip-loss factor. Section 4 shows the conclusions of this study.

2. General formulations for the performance coefficients General momentum theory differs from the simple axial momentum theory that leads to the Betz-Joukowsky limit, by accounting for the rotational motion. To achieve this, it is necessary modify the usual actuator disc theory. A flow model that includes the angular momentum, was presented in Joukowsky [16], and applied by Glauert [7] for the study of propellers, and later modified by Wilson and Lissaman [8] for HAWTs. In that case, the induction factor in the far-wake has a non-linear relationship

664

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

with the induction factor at the rotor plane. These induction factors are the induced velocities normalized by the wind speed, V 0 . Wilson and Lissaman [8] provide a more general one-dimensional mathematical model than Glauert’s [7], for the calculation of the thrust, torque and power coefficients. Sørensen and Kuik [11] used the general momentum theory to study the behavior of the classical free vortex wake of Joukowsky [16], and confirm some results found by Wilson and Lissaman [8], as described in Rio Vaz et al. [6]. In this section, the general performance coefficients are presented. Fig. 1 illustrates the rotor and wake and the flow axial velocities. The axial induced velocities u and u1 at the rotor plane and in the wake, respectively, are written as




V 0  v ¼ u  ð1  aÞV 0

V 0  v 1 ¼ u1  ð1  bÞV 0

ð1Þ

In the simple axial momentum theory, the axial induction factor in the wake is twice that at the rotor. At low X, Joukowsky’s [16] analysis of the momentum, energy and mass conservation for the control volume illustrated in Fig. 1, b can be expressed as a function of a as

a¼

" # 2 b b ð1  aÞ 1 2 2 4X ðb  aÞ

ð2Þ

It is important to note that Eq. (2) is a result of Joukowsky’s assumption of constant bound circulation. The continuity equation is

udA ¼ u1 dA1

ð3Þ

Conservation of circulation gives

uh r ¼ uh1 r 1 ¼ C=ð2pÞ

ð4Þ

The thrust on the blades is found using the Kutta-Joukowsky theorem. It matches the pressure drop across the rotor caused by the swirl in the wake. Applying the Kutta-Joukowsky theorem and ignoring the shear stress and lateral pressure at the boundaries of the annular streamtube intersecting B blade elements, the incremental thrust, dT, on the element is

   uh  1 uh 2Xr uh dA dT ¼ qC Xr þ þ dr ¼ qV 20 2 2 V0 V0 V0

ð5Þ

Thus, the incremental form of the thrust coefficient, dC T , becomes

dC T ¼ uh; ð2Xr  þ uh; Þ

According to the general momentum theory with the rotor represented by an actuator disc, C ¼ 2pruh , the dimensionless circulation q is expressed in terms of the swirl velocity, as described by Sørensen and Kuik [11].

q

 r  uh;

ð7Þ

The local thrust coefficient assumes the form:

  q dC T ¼ q 2X þ 2 r

ð8Þ

Thus dC T  1=r  as r  # 0 as observed by Mikkelsen et al. [12]. The thrust coefficient for the rotor, C T , is obtained integrating Eq. (5):

CT ¼ 1 2

Z

T

q

V 20 A

1

¼2 0

uh; ð2Xr þ uh; Þr  dr 

ð9Þ

The torque dQ on the annular element is found using the balance between angular momentum and the blade torque, determined from the Kutta-Joukowsky equation. The local torque can be obtained directly from Euler’s turbine equation Hansen [21], Sørensen and Kuik [11], and is given by

dQ ¼ quuh rdA

ð10Þ

In elemental form:

dQ

dC Q ¼ 1

qV 20 RdA 2

¼ 2u uh; r  ¼ 2ð1  aÞuh; r 

ð11Þ

The torque coefficient for the rotor is found by integrating across the rotor:

CQ ¼ 1 2

Q

¼1

q

V 20 RA

2

Z

1

q

V 20 RA

R

0

quuh r2prdr ¼ 4

Z 0

1

ð1  aÞuh; r2 dr  ð12Þ

Since only steady operation is considered, dC P ¼ XdC Q . In terms of the dimensionless circulation q, inserting Eq. (7) in (11), the local power coefficient becomes that obtained by Sørensen and Kuik [11]:

dC P ¼ 2Xð1  aÞq

ð13Þ

Note that, if the circulation is independent of r, then a and q are constant, and the latter is given by

ð6Þ

where r  ¼ r=R and uh; ¼ uh =V 0 . Throughout this paper, the subscript  indicates normalization by the blade tip radius and wind speed. Eq. (6) shows that the thrust is entirely due to the swirl in the wake at low X. As X increases, 2uh; Xr remains finite while u2h; # 0.

C

2pRV 0

2

q¼

b ð1  aÞ 2Xðb  aÞ

ð14Þ

The turbine power coefficient is

C P ¼ XC Q ¼ 1 Z

2

P

q

V 30 A

¼1 2

1

q

V 30 A

Z 0

R

qXuuh r2prdr

1

¼ 4X 0

ð1  aÞuh; r 2 dr 

ð15Þ

Combining Eqs. (13) and (14), give the result found by Wilson and Lissaman [8] for the power coefficient. This result is obtained if uh =V 0 ¼ q=r  is substituted in Eq. (15) and a is assumed constant. Then the blade element and turbine power coefficients are equal: 2

dC P ¼ C P ¼

b ð1  aÞ2 : ðb  aÞ

ð16Þ

2.1. Joukowsky model with a Rankine vortex

Fig. 1. Simplified illustration of the velocities in the rotor plane and in the wake.

The wake of a Joukowsky rotor contains a free vortex system consisting of helical vortices trailing from the tips of the blades

665

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

and a rectilinear hub vortex [22,23]. Hence, as described by Mikkelsen et al. [12], uh  1=r, which can be avoided by introducing a hub vortex with a cut-off radius, rc , such that uh # 0 at r # 0, corresponding to a finite hub vortex. The Rankine vortex has uh  r for r  < r c; and uh  1=r for r P r c; . Using n ¼ r=r c , the Rankine vortex is defined as

(q r q r

uh; ¼

dC T ¼

2

n ; 06n61

ð17Þ

; 16n

8 h i > < rq n2 2Xr  þ rq n2 ; 0 6 n 6 1   dC T ¼ > : rq 2Xr  þ rq ; 1 6 n  

ð18Þ

2ð1  aÞqn2 ; 0 6 n 6 1 2ð1  aÞq; 1 6 n

ð19Þ

To determine the turbine thrust and torque coefficients, it is necessary to insert Eq. (17) in Eqs. (9) and (12). For the thrust coefficient:

)    Z 1  q 2 q 2 q q r  dr  CT ¼ 2 n 2Xr þ n r  dr  þ 2Xr þ r r r 0 rc; r  i  o 1 n h ¼ q q 1  ln r 4c; þ 2X 2  r2c; 2 ð20Þ r c;

" ðaþ1Þ 4 2Xr þ

 ðaþ1Þ # aþ1 4 qn rc; a þ n4

j



aþ1 a þ n4

ðaþ1Þ 4

ð27Þ

j

aþ1

qða þ 1Þ 4 r c;      9 8 12 a 34 a aþ1 1a 3a > 2a 4 > >  r 2c; a 2 4r c; X rac; 1 þ r4c; a  r3c; a 4 > = > a a 1  3  > > ; :

ð28Þ where the first term in the brackets is zero when a ¼ 1, and

4ð1  aÞjqða þ 1Þ r c; ð3  aÞ

aþ1 4

   34 a 3a  r 3c; a 4 r ac; 1 þ r 4c; a

C Q ¼ 4ð1  aÞ 0

r c;

q 2 2 n r dr  þ r 

 ¼ ð1  aÞq 2  r 2c;

Z

1

r c;

q 2 r dr  r 

Delery’s [18] review of vortex breakdown included an empirical representation of the swirl velocity distribution. In the most useful form for present purposes:

uh; ¼

q 1  exp½1:256n2  r

ð30Þ

Using Eq. (30), the elemental thrust and torque coefficients are obtained from Eqs. (6) and (11) as

dC T ¼






q q 1  exp½1:256n2  2Xr  þ 1  exp½1:256n2  r r

ð31Þ

#

and, for constant a,



dC Q ¼ 2ð1  aÞq 1  exp½1:256n2  ð21Þ

2.2. Joukowsky model with a Vatistas vortex

uh uh;max

ð22Þ

ð24Þ

where j is a parameter to be determined later. Rearranging Eq. (22), and introducing Eq. (24), yields:

uh; ¼

 ðaþ1Þ  ðaþ1Þ uh;max aþ1 4 j aþ1 4 n ¼ qn V0 r c; a þ n4 a þ n4

"

1

! #) 1:256  1 r2c;

2 q2 1  exp½1:256n2  dr  r

ð33Þ

for the thrust and, for the torque,

C Q ¼ 4ð1  aÞq

( " ! #) 1 1:256 1 : þ 0:398089r 2c; exp  2 2 r c;

ð34Þ

2.4. Betz-Goldstein model

ð23Þ

Combining Eqs. (7) and (23) gives

uh;max j C j ¼ ¼ q V0 r c; 2pRV 0 r c;

Z þ2 0

where uh;max is the maximum swirl velocity which occurs at n ¼ 1. Vatistas [17] suggested that the best value of the parameter a is 0.75. However, that value is not consistent with a constant circulation outside the vortex which requires a ¼ 1. When n ¼ 1:

r c uh;max  C=ð2pÞ

(

C T ¼ 2Xq 1 þ 0:796178r 2c; exp 

Vatistas [17] developed an alternative formulation for the swirl in a generic vortex, which is given by

ð32Þ

Combining Eqs. (9), (12) and (30), the integral forms are found as

 ðaþ1Þ aþ1 4 ¼n a þ n4

ð29Þ

2.3. Joukowsky model with a Delery vortex

When a is constant, the torque coefficient is given by

"Z

ð26Þ

The integral forms are

CQ ¼

and

(Z

aþ1 a þ n4

dC Q ¼ 2ð1  aÞjqn2

CT ¼

dC Q ¼

r c;

qn

and, for constant a,

The elemental thrust and torque coefficients, considering the Rankine vortex, are given by inserting Eq. (17) in Eqs. (6) and (11), respectively:

(



j

ð25Þ

The elemental thrust and torque coefficients are given by inserting Eq. (25) in Eqs. (6) and (11), respectively. Thus:

Wood [15] analyzed a HAWT with B ¼ 1 as X decreases to zero, using the BG analysis, modified to give positive torque on a stationary rotor. He demonstrated that the modified BG torque coefficient never exceeds 0.55 and the power is less than the Glauert optimum. The axial and swirl velocity distributions were taken from Okulov and Sorensen [9]. The BG model approaches the Joukowsky model at large X, and has two huge advantages that are closely connected: (a) it is not necessary to postulate a separate vortex core, and (b) when r # 0; uh # 0. Thus the BG model should avoid the problems with the Joukowsky model at low X due to its assumption of constant circulation. The swirl velocity distribution behind a BG rotor from Wood [15] is

666

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

uh; ¼ 2

wr  p p2 þ r 2

ð35Þ

where w ¼ w=V 0 is the normalized velocity of the vortex sheet formed behind the HAWT, and p ¼ ð1  wÞ=X is the dimensionless linear pitch of the helical vortex. p is constant with radius, which is the distinguishing feature of the BG rotor. The axial velocity at the rotor is

u ¼ 1 

wr 2 þ r2

ð36Þ

p2

and the dimensionless circulation q, is found by substituting Eq. (35) in Eq. (7):

q¼2

wr2 p p2 þ r 2

ð37Þ

Eqs. (35)–(37), that the velocity make the field and circulation now dependent on r for the whole rotor. Using Eqs. (6) and (11) in terms of the new formulations for u and uh; , yield the elemental thrust as

  wr 2 p wp dC T ¼ 4 2  2 X þ 2 2 p þ r p þ r

ð38Þ

and the elemental torque is

  wr2 wr2 p dC Q ¼ 4 1  2  2 2  2 p þ r p þ r

ð39Þ

ð40Þ


  2  pw p 2 2 1  w þ p ð1  2wÞ 1 þ ð1 þ p Þ ln 1 þ p2 1 þ p2

The velocity of the vortex sheet, w, in the BG optimization as described by Wood [15], is found as the root of

h

3

4

i

X X ð1  3wÞ þ ð1  wÞ ð8w  6w þ 1Þ h i2  2X 2 fw½ð6w  11Þw þ 6  1g  X 2 þ ð1  wÞ2 " # X2 ¼0  ð1  wÞð1  2wÞð1  4wÞ ln 1 þ ð1  wÞ2 2

Z

2

ðu1  V 0 Þ dA1 ¼ A1

A1

The terms

uh ur

and

uh1 u1 r 1

r 1 uh1 u1

2X  urh 2X   u u1

uh1 r1

ð42Þ

dA1

ð43Þ

can easily be replaced through mathemat-

u1;

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi Xq Xq  2ð1  u Þ u u

uiter¼1 1;

and riter¼1 . In this work 1; iter¼2

¼ u1; ðr ¼ 0Þ from Eq. (44), dr 1;

¼ dr  (at the rotor

¼ 0; plane) and riter¼1 1; for i ¼ 2 to Ns (Number of sections) do while error > TOL do iter ¼ iter þ 1 iter1

iter1 þ dr 1; Compute riter 1; ¼ r 1;

uiter 1; ,

;

using Eq. (45); uiter r iter

Compute the new dr 1; ¼ uiter riter dr  using the continuity 1; 1;

equation;

end while Compute ui1; ¼ uiter 1; . end for

Thus, the new axial induction factors at the rotor plane and in the wake are:

wr 2 þ r 2

p2

  q q b ¼ 1  u1; ¼ X þ 2 2r  u vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u   2 u q q q 1 1 t 2  Xþ 2  2X ð1  u Þ þ q  2r u u r 21; r2

ð46Þ

ð47Þ

For future use, it is noted that r c for the Rankine, Vatistas, and Delery vortices is the radius that maximizes both uh and the angular momentum flux. In the vortex breakdown analysis, rc will be used in the same manner for the BG model.

!

ical arrangements as shown by Sørensen and Kuik [11] using the continuity equation and the conservation of circulation, Eqs. (3) and (4) respectively, into Eq. (43), resulting in

Xq ¼1 þ u

iter¼2

Set initial values for uiter¼1 ; dr 1; 1;

a ¼ 1  u ¼

For all models of the rotor loading and the velocities immediately behind the rotor, special attention needs to be given to the wake expansion, since models based on vortex theory, represent the ideal far wake as a helical surface. It is necessary to find a consistent formulation for the axial velocity in the far wake, u1 . It is now shown that the axial momentum balance on the whole control volume (see Fig. 1) leads to a suitable relationship between axial and swirl velocities, which is given by

Z

Algorithm 1. Iteractive procedure to calculate u1; using Eq. (45).



iter iter1 Compute error ¼ dr 1;  dr 1; .

ð41Þ

2

ð45Þ

In this case, u1; is solved using a simple iteractive procedure, as described in Algorithm 1, where Eq. (45) is combined with the continuity equation ur  dr  ¼ u1 r 1; dr 1; , to calculate dr 1; .

iter

and the torque is

CQ ¼ 4

  q q u1; ¼ 1  X þ 2 2r  u vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u  2 u q q q 1 1 þt Xþ 2  2X ð1  u Þ þ q2 2  2 2r  u u r1; r 

Compute

By integration, Wood [15] found the turbine thrust as


  2  w p C T ¼ 4pw X  p  ðXp  wÞ ln 1 þ p2 1 þ p2

It is noteworthy that Eq. (44) is valid only if u and u1 are independent of r, leading to the removal of uh =ðurÞ and uh1 =ðu1 r 1 Þ from Eq. (43). In this study, the dependence on r is important, so a new formulation is used for the axial velocity in the far wake from Eq. (43) as

ð44Þ

3. Results and discussion The models described above are used in this section to determine the local and integral forms for the performance coefficients. Section 3.1 shows the blade element behavior of the formulations using all the core models and the BG analysis. In Section 3.2, the integral forms are compared to the results of Mikkelsen et al. [12], who studied the flow through rotors using a CFD actuator line model. In Section 3.3, the results of Okulov and Sørensen [9] are also utilized to assess the present approach. The BG analysis for infinite B is corrected by employing Prandtl’s approximation to

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

the tip loss function and the newer calculation of Wood et al. [20] for finite B. 3.1. Blade element results Here, the local formulations are evaluated. Thus, Fig. 2 depicts the axial induction factors for the BG model at the rotor plane and in the wake, respectively. It is evident that a and b are heavily influenced by the radial variation at low X (Fig. 2a and b). The same observation holds for the area expansion, which was calculated using the continuity equation (dA1 =dA ¼ u =u1; ), Fig. 3a. The physical consistency of the BG model is seen by all results converging to the values appropriate at the Betz-Joukowsky limit as X ! 1. Consequently, although the circulation changes strongly with r and X

667

at low X (Fig. 3b), the area expansion remains below the BetzJoukowsky limit. As previously described, the free vortex gives infinite swirl velocities on the axis, leading to infinitely high local thrust coefficient, dC T ðr # 0Þ ! 1. Fig. 4 confirms this behavior. These results were obtained using Eq. (6) with the Joukowsky model combined with Rankine, Vatistas [17] and Delery [18] vortices to regularize the swirl velocity using r c; ¼ 0:15. The calculations used b ¼ 0:7 at X ¼ 0:5, since the axial induction factor in the wake must be in the range 0 6 b 6 1. For X ¼ 7 the ideal optimum rotor value b ¼ 2=3 was used. Fig. 4a shows the strong deviation of the Joukowsky blade element results from the Betz-Joukowsky limit at X ¼ 0:5. This deviation is not necessarily a fault of the model. It is possible for wind turbines to exceed the Betz-Joukowsky limit

Fig. 2. Radial profiles of the axial induction factor (a) at the rotor plane and (b) in the wake for the BG model.

Fig. 3. (a) Radial profiles of the wake expansion and (b) the dimensionless circulation for the BG model.

Fig. 4. Elemental thrust coefficient distribution, (a) X ¼ 0:5 and (b) X ¼ 7:0.

668

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

on thrust, e.g. [24]. In the BG formulation, however, maximizing C P requires that C T remain below the limit of 8/9. The Joukowsky model gives dC T ! 1 as r # 0. On the other hand, for X ¼ 7 the Joukowsky model combined with the Vatistas [17], Delery [18] and Rankine vortexes avoids this problem, Fig. 4b. For the Vatistas vortex the parameters a ¼ 1 and j ¼ 0:75 were used. The third method to avoid the infinite swirl velocity near the axis, the Rankine vortex, has the elemental thrust given by Eq. (18). As for the other ‘‘core” models, dC T is zero on the axis. Even though the Rankine core is used, the solution does not converge at X ¼ 0:5 but does so at X ¼ 7, except around r  ¼ 0:2, Fig. 4b. The Vatistas [17] and Delery [18] swirl distributions, Eqs. (26) and (31), are similar. Differently, the BG model does not need a vortex core, and ensures a natural evolution for the swirl velocity with radius, causing dC T ðr # 0Þ ! 0 differently to the Joukowsky [16] and Delery [18] models. Fig. 5 shows the results for the local power coefficient. As expected, for low X the difference is considerable, while for high X the differences between the models does not significantly influence the power. The explanation for this is that the tangential velocity ðXr þ uh; =2) has a strong dependence on r at low X, but when X is large, the term Xr  overwhelms the effect of the swirl velocity, leading all models to converge. 3.2. Turbine thrust and power An assessment of the integral equations is presented in this subsection. Fig. 6a shows the non-linear relationship between the axial induction factors b=a, as a consequence of wake rotation in the Joukowsky model [16] without a vortex core, which clearly is not accurate at low X. The difference with the BG model (Fig. 6b) is notable. The results in Fig. 6b were obtained using the

area-averaged velocity calculated by the continuity equation. Fig. 7a shows the results for the circulation, q, obtained using Eq. (14). Note that the Joukowsky model gives the induced velocity field as independent of r, and thus uh becomes infinite as r # 0. On the other hand, uh in the BG model follows the circulation (Fig. 7b) remains finite and avoids the singularity. Fig. 7b was genR1 erated using the average qav ¼ 2 0 qr dr  , since q given by Eq. (37) is r-dependent. The infinite behavior for the Joukowsky model is also seen in Fig. 8a for C P . In contrast, the BG model gives C P (Fig. 8b) below Betz-Joukowsky limit, at all X. These results show that the BG model is well-behaved for all X. The results in Fig. 9a were calculated by assuming that the lifting part of the rotor disc extends from r  ¼ r c; to r ¼ 1 to avoid the singularity in r ¼ 0 in the Joukowsky model. In this case, a vortex core rc; ¼ 0:2 was used. Fig. 10 depicts the sensitive of rc; , where C T must decrease as shown through Eq. 9, in which C T / max ðuh Þ2 , and uh / 1=r c; . Fig. 9a shows consistent behavior for all models using the integrated form of the thrust coefficient, with reasonable agreement when compared with Mikkelsen et al. [12], who used both CFD actuator line model and an axi-symmetric CFD model. In Fig. 9b, all vortex models tend to the Betz-Joukowsky limit. The torque coefficient has a maximum of approximately C Q ¼ 0:5 at X  0:4 for the Joukowsky, Joukowsky with Rankine vortex, Vatistas and Delery formulations. For the BG model the maximum value obtained is C Q ¼ 0:5 at X  0:5. This result agrees with Wood [15]. Mikkelsen et al. [12] raised the possibility of vortex breakdown occurring immediately behind the rotor. Breakdown can cause a substantial increase in the vortex radius, rc , which would invalidate the velocity distributions that have been assumed in this and previous analyses and may set a limit on wind turbine

Fig. 5. Elemental power coefficient distribution, (a) X ¼ 0:5 and (b) X ¼ 7:0.

Fig. 6. (a) b=a as a function of X from the Joukowsky model without a vortex core. (b) bav =aav as a function of X from the BG model.

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

669

Fig. 7. (a) Circulation dependence on b=a for the Joukowsky model without a vortex core. (b) Circulation dependence on bav =aav from the BG model.

Fig. 8. (a) Power coefficient as a function of axial induction factors a and b. (b) Power coefficient as a function of the average axial induction factors aav and bav from the BG model.

Fig. 9. Turbine (a) thrust and (b) power coefficients for different X.

performance. According to Delery [18] and Mikkelsen et al. [12], breakdown appears when the circulation, q, is about 1.4 greater than u rc; . Alternatively, in terms of the swirl number, S, which is the normalized ratio of the angular momentum to axial momentum flux:

R1 S¼

qu uh 2pr dr > 1:4 R1 r c; 0 qu2 2pdr  0

ð48Þ

Thus, using the Joukowsky model, Eq. (48) assumes the form described by Delery [18] and Mikkelsen et al. [12],

SJ ¼

q > 1:4: rc; u

ð49Þ

For the Joukowsky model with a Rankine vortex, the swirl number is

SR ¼

  q 2 1  r c; : r c; u 3

ð50Þ

For the Joukowsky model with a Vatistas vortex,

2  !3 1þa q 4k 1 þ 1a 4 3 1þa 7 1 5 : ; SV ¼ ; ; 4 2F1 r c; u 4 4 4 r c; a 3r 2c;

ð51Þ

2 F 1 is the hypergeometric function documented in Chapter 15 of the NBS Handbook [25] where it is called Fða; b; c; zÞ, but F is reserved in this paper for the tip-loss factor discussed below. For the Joukowsky model with a Delery vortex,

SD ¼

   q 1:12071 1  0:79077r c; erf r c; u r c;

ð52Þ

670

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

commonly-used form of F [26]. Blade element calculations made with Prandtl’s model show good agreement with test data [27] at high X, however, Wood et al. [20] showed that it is inaccurate at low X and does not reproduce the increase in F above unity near the axis for any X. Shen et al. [28], for example, put the tip loss in the form F ¼ a=ab , where a and ab are the streamtube average axial interference factor, and the factor at the blade, respectively. Glauert [7] introduced the tip loss as a correction to a and a0 , through the axial and angular momentum equations. Matching the blade element thrust and torque with the streamtube momentum and angular momentum gives

  dT ¼ 4pqX2 r 3 a0b F 1 þ a0b dr

ð54Þ

and

dQ ¼ 4pqV 0 Xr 3 a0 ab Fð1  ab Þ dr

ð55Þ

Eqs. (54) and (55) become

Fig. 10. Sensitivity of the thrust coefficient C T to r c; .

where erf is the standard error function. For the BG model, " #  

2pw 2ðw  1Þ þ pð3w  2Þ p  ð1 þ p2 Þcot1 p 1 SBG ¼

r c; 2ðw  1Þ2 þ p p½ð4  3wÞw  2 þ wð1 þ p2 Þð3w  4Þcot1 p ð53Þ

Fig. 11 shows the results for all swirl numbers, Eqs. (49)–(53). The axial induction factor at the rotor plane was calculated using Eq. (2). For the Joukowsky, Rankine, Vatistas and Delery swirl numbers the values used by Mikkelsen et al. [12] were utilized, b ¼ f0:19; 0:35; 0:60g and X ¼ f0:5; 1:0; 2:0g. Two more points at X ¼ f0:25; 0:65g; b ¼ f0:088; 0:25g were also considered. Note that the BG flow never breaks down even at low X. For Joukowsky with Rankine, Vatistas and Delery swirl distributions, S is reduced so that breakdown only around X ¼ 0:65. Such results demonstrate that the swirl number is strongly influenced by rc; . 3.3. Finite number of blades So far the analysis has been of rotors with an infinite number of blades. The easiest extension to finite B is through the use of a tip-loss factor, F. Okulov and Sørensen [9] extended the BetzJoukowsky limit to finite B using Prandtl’s approximation [19], F P . F is formally defined as the ratio between the actual bound circulation and that of a rotor with B ¼ 1. F P is the most

dT ¼

1 qV 20 uh; ðXr þ uh; ÞFdA 2

ð56Þ

and

dQ ¼ quuh rFdA

ð57Þ

where Prandtl’s approximation F P is

FP ¼

2

p

cos1 ðef Þ

ð58Þ

in which

f ¼

B ð1  r  Þ 2

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ p2 p

ð59Þ

C T ; C Q and C P ¼ XC Q were solved using the trapezoidal method due to the complexity of Eqs. (38) and (39) when F is included. Results for B ¼ f1; 3; 1g were obtained using 100 blade elements in the interval of 0 6 X 6 15. Fig. 12 shows the results are in good agreement of those of Okulov and Sørensen [9]. The power coefficient, converges to the classical Glauert’s theory [7] for B ¼ 1 (Fig. 12a). The thrust coefficient, likewise increases to the BetzJoukowsky value of 8/9 as X increases, Fig. 12b. Fig. 12c shows that the torque coefficient has a maximum of 0.5 at X ¼ 0:5 for B ¼ 1; this value decreases for B ¼ 3 and 1. Okulov and Sørensen [9] did not present results for the torque coefficient. Glauert’s [7] results provide a comparison with B ¼ 1 at low X. In order to evaluate

Fig. 11. Swirl number, S, as a function of X.

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

671

Fig. 12. (a) Power coefficient, C P , (b) thrust coefficient, C T , and torque coefficient, C Q as function of X for different numbers of blades.

the accuracy of F P at low X, the tip-loss calculation method recently described by Wood et al. [20] is shown in Fig. 12. These calculations find F directly from the equations for the velocities induced by the helical trailing vortices. Three separate calculations are described in Wood et al. [20], of which Okulov’s equations were found to be the most accurate. F was generated only for B ¼ 3 due to the complexity in calculating the Goldstein’s function for finite B. The results show that F P leads to an error in C P of 0.44% at X ¼ 1:5, in relation to the result of Okulov and Sørensen [9], and 10% in comparison to the result of Wood et al. [20]. It is worth noting one further aspect of finite B: the effect of the increase in F near the axis of rotation documented by Wood et al. [20]. Adding F to the definition of S gives

Fig. 13. Tip-loss factor, F, and the ratio, SBG;3 =SBG;1 , as functions of radius, r  , for X ¼ f0:65; 7:0g.

SBG;B F  wFq > 1 if F > 1: ¼ SBG;1 1  wFq

ð60Þ

Thus it is possible that the BG model will lead to vortex breakdown for finite F as SBG;B is significantly larger than the infinite B near the axis, Fig. 13. However, a value of F significantly different from unity implies significant azimuthal variations in the velocities and the effects of these on vortex breakdown has, to the authors’ knowledge, not been determined.

4. Conclusions Analytical formulations have been developed for the performance coefficients of HAWTs considering the swirl in the flow downwind of the rotor. It is well known that the general magnitude of swirl increases with decreasing X, making that region a difficult one to model accurately and simply. The initial analysis considered the Joukowsky model of constant bound circulation modified by the vortex models of Rankine, Delery and Vatistas placed on the axis of rotation to avoid the singularity in the swirl velocity that would otherwise occur. Equations for the radiallyvarying thrust, torque and power were derived form general momentum theory using the induced velocity field from the vortex models. The Joukowsky model with a Rankine vortex produces unrealistic results for the elemental thrust near the axis. An alternative formulation – the BG model – was also shown to avoid the singularity on the axis and to provide consistent distributions of the thrust. This model produces no torque at zero X, but this deficiency can be simply removed. The torque coefficient has a maximum value, using Eq. (12), for all swirl velocity distributions described here. This value corresponds approximately to C Q ¼ 0:5, in agreement with the results of Wood [15]. The radially-varying axial and swirl velocities of the BG model was the only model found to be safe from vortex breakdown. The

672

J.R.P. Vaz, D.H. Wood / Energy Conversion and Management 126 (2016) 662–672

general, consistent behavior of the BG model is the most important result of this investigation. Finally, corrections for finite blade number using Prandtl’s tiploss factor [19] was included. The results were compared with classical Glauert theory [7] and those of Okulov and Sørensen [9]. The present calculations of the thrust, torque, and power coefficients are physically consistent. Using the more accurate tip loss calculation of Wood et al. [20] for a three-bladed rotor produced a significant change in all parameters at low X. Acknowledgments The authors would like to thank the NSERC Industrial Research Chairs program, the ENMAX Corporation, CNPq, CAPES, Eletronorte, PROPESP/UFPA for financial support. References [1] Ashuri T, Zaaijer MB, Martins JRRA, Zhang J. Multidisciplinary design optimization of large wind turbines – technical, economic, and design challenges. Energy Convers Manage 2016;123:56–70. [2] Göçmen T, van der Laan P, Réthoré PE, Diaz P, Chr.Larsen G, Ott S. Wind turbine wake models developed at the technical university of Denmark: a review. Renew Sustain Energy Rev 2016;60:752–69. [3] Sørensen JN. General momentum theory for horizontal axis wind turbines. Springer-Verlag London; 2016. [ISBN: 978-3-319-22113-7]. [4] van Kuik GAM, Sørensen JN, Okulov VL. Rotor theories by Professor Joukowsky: momentum theories. Prog Aerosp Sci 2015;73:1–18. [5] Wood DH. Including swirl in the actuator disk analysis of wind turbines. Wind Eng 2007;31(5):317–23. [6] Rio Vaz DATD, Vaz JRP, Mesquita ALA, Pinho JT, Brasil Junior AP. Optimum aerodynamic design for wind turbine blade with a Rankine vortex wake. Renew Energy 2013;55:296–304. [7] Glauert H. Aerodynamic theory. In: Durand WF, editor. Division L. Airplanes propellers, vol.4; 1935. p. 191–5 [reprinted, Dover, NewYork, 1963, chapter XI].

[8] Wilson RE, Lissaman PBS. Applied aerodynamics of wind power machines, tech rep. Oregon State University; 1974. [9] Okulov V, Sørensen JN. Refined Betz limit for rotors with a finite number of blades. Wind Energy 2008;11:415–26. [10] Goldstein S. On the vortex theory of screw propellers. Proc Roy Soc Lond, Ser A 1929;123:440–65. [11] Sørensen JN, van Kuik GAM. General momentum theory for wind turbines at low tip speed ratios. Wind Energy 2011;14:821–39. [12] Mikkelsen RF, Sarmast S, Henningson D, Sørensen JN. Rotor aerodynamic power limits at low tip speed ratio using CFD. J Phys: Conf Ser 2014:524. http://dx.doi.org/10.1088/1742-6596/524/1/012099. [13] Wood DH. Ideal wind turbine performance at very low tip speed ratio. J Renew Sustain Energy 2014;6:033115-1–033115-12. http://dx.doi.org/10.1063/ 1.4879276. [14] Goldstein S. On the limiting values for infinite pitch of a parameter occurring in airscrew theory. Proc Cambridge Philos Soc 1944;40:146–50. [15] Wood DH. Maximum wind turbine performance at low tip speed ratio. J Renew Sustain Energy 2015;7:053126. http://dx.doi.org/10.1063/1.4934805. [16] Joukowsky NE. Vortex theory of a rowing screw. Trudy Otdeleniya Fizicheskikh Nauk Obshchestva Lubitelei Estestvoznaniya; 1912. p. 16 [in Russian] [17] Vatistas GH. Simple model for turbulent tip vortices. J Aircraft 2006;43 (5):1577–9. [18] Delery J. Aspects of vortex breakdown. Prog Aerosp Sci 1994;30:1–59. [19] Prandtl L, Betz A. Four essays on the hydrodynamics and aerodynamics (original German title: Vier Abhandlungen zur Hydrodynamik und Aerodynamik). Göttinger Nachr.: Göttingen; 1927. [20] Wood DH, Okulov VL, Bhattacharjee D. Direct calculation of wind turbine tip loss. Renew Energy 2016;95:269–76. [21] Hansen MOL. Aerodynamics of wind turbines. 2nd ed. London, UK: Earthscan; 2008. [22] Okulov V, Sørensen JN, Wood DH. The rotor theories by Professor Joukowsky: vortex theories. Prog Aerosp Sci 2015;73:19–46. [23] Dobrev I, Maalouf B, Troldborg N, Massouh F. Investigation of the wind turbine vortex structure. In: 14th Int symp on applications of laser techniques to fluid mechanics Lisbon, Portugal, 07–10 July; 2008. [24] van Kuik GAM, Lignarolo LEM. Potential flow solutions for energy extracting actuator disc flows. Wind Energy 2016;19:1391–406. [25] Abramowtiz M, Stegun I. Handbook of mathematical functions. Dover; 1964. [26] Wald QR. The aerodynamics of propellers. Prog Aerosp Sci 2006;42:85–128. [27] Wilson RE. Aerodynamic behavior of wind turbines. In: Spera DA, editor. Wind turbine technology. New York, USA: ASME Press; 2009. [28] Shen WZ, Mikkelsen R, Sørensen JN, Bak C. Tip loss corrections for wind turbine computations. Wind Energy 2005;8:457–75.