**UFRJ**

# Aerofólio S809 - 1 - s2 0 - s0960148112004156 - main

(Parte **1** de 2)

Advanced brake state model and aerodynamic post-stall model for horizontal axis wind turbines

R. Lanzafame, M. Messina* Department of Mechanical and Industrial Engineering, University of Catania, Viale A. Doria, 6, 95125 Catania, Italy art i cle i nfo

Keywords: BEM theory Centrifugal pumping Brake state model NREL wind rotor abstract

In scientiﬁc literature, when the aerodynamic design of a horizontal axis wind turbine is discussed, different brake state models are presented. The brake state models are implemented within a BEM code which is a 1-D numerical code, based on Glauert propeller theory, and able to predict HAWT performance. This code provides reliable results only if a proper brake state model and aerodynamic post-stall model are implemented.

In this research, the authors have produced a numerical code based on BEM theory in conjunction with an aerodynamic post-stall model, indispensable for taking into account radial ﬂow along the wind turbine blades (Himmelskamp Effect), and the brake state models by Buhl, combined with the calculation of Jonkman’s tangential induction factor.

In scientiﬁc literature, Shen’s brake state model is commonly implemented within 1-D numerical codes, based on BEM theory. Subsequently, a comparison with Shen’s brake state models was carried out. With the numerical code presented in this work, the power for an NREL wind rotor was predicted. With the numerical simulation, it was possible to notice when these different brake state model furnish results close to experimental data. 2012 Elsevier Ltd. All rights reserved.

1. Introduction

Thenumerical codes based on BEM theoryare powerful tools for the design and performance evaluation of HAWT. BEM theory is based on Glauert propeller theory [1] modiﬁed for wind turbines. Today, many researchers are developing numerical codes based on BEM theory [2e11]. Industry also utilizes these numerical codes to design HAWT. These numerical codes are 1-D codes and produce reliable results provided certain criticalities are resolved. These criticalities regard the correct representation of lift and drag coefﬁcients at high values of angle of attack, the implementation of a post-stall model (to take into account radial ﬂow along the blades) and the implementation of a brake state model (to correctly determine axial and tangential induction factors) [12e17].

This paper compares two different (the most accredited) brake statemodelstoevaluatetheperformanceofaHAWT.Thetwobrake state models are Shen’s [18,19] and Buhl’s [20,21] (here Buhl’s model is combined with Jonkman’s equation to determine the tangential induction factor).

First, a numerical code based on BEM theory [13] was developed, and a post-stall model [13,17] was implemented within the numerical code. Next, the two brake state models were compared, predicting the power curves for the NREL wind rotor [2].I n scientiﬁc literature [29], experimental measurements are reported for this wind turbine rotor. Finally, a comparison between the simulated and experimental power curves is performed.

2. BEM theory, post-stall model and brake state model

The numerical code, developed in [13] is a 1-D code for the design of a horizontal axis wind turbine. It has very fast computational times and provides great accuracy compared with experimental data. This code is based on Blade Element Momentum (BEM) theory, and can be implemented to design a wind rotor, and/ or evaluate its performance.

BEM theory based numerical codes subdivide the wind turbine rotor into annuli of dr thickness, the ﬂow of each sector being independent of adjacent circular sector ﬂows [17,23]. Applying the equations of momentum and angular momentum conservation, for

Abbreviations: BEM, Blade Element Momentum; HAWT, horizontal axis wind turbine; NREL, National Renewable Energy Laboratory; BSM, brake state model; 1- D, One-dimensional; 2-D, Two-dimensional. * Corresponding author. E-mail address: mmessina@diim.unict.it (M. Messina).

Contents lists available at SciVerse ScienceDirect

Renewable Energy journal homepage: w.els evier.com/locate/renene each inﬁnitesimal dr sector of the blade, the axial force and torque can be deﬁned (Eqs. (1) and (2)). The axial force on the blade element of width dr is:

sin2f NbðCLcos f þ CD sin fÞ c dr (1)

The torque on the blade element of width dr is:

sin f , u r ð1 þ a0Þ

cos f Nb ðCLsin f CD cos fÞ cr dr

Fig. 1 shows the axial and tangential forces (dN and dT) for an annulus of width dr.

From Eq. (2) wind rotor power can be evaluated, as reported in Eq. (3).

P ¼ Z u dT (3) while the power coefﬁcient is given by Eq. (4)

2.1. Post-stall model

Knowing the lift and drag coefﬁcients (CL and CD) is crucially important for assessing the forces and torques according to Eqs. (1)

A problem which might cause numerical instability is linked to the mathematical description of the airfoil lift coefﬁcient. The airfoil’s experimental data is 2-D and taken from wind tunnel measurements. Furthermore, due to rotation, the boundary layer is subjected to Coriolis and centrifugal forces which alter the 2-D airfoil characteristics. This is especially pronounced in stall. It is thus often necessary to extrapolate existing airfoil data into deep stall and include the effect of rotation [24e28].

Owing to radial ﬂow along the turbine blades, mathematical equations describing lift coefﬁcient have to overestimate experi- mental CL values within a precise range of the angle of attack. Centrifugal pumping is a phenomenon which describes radial air ﬂow along blades [29,30]. This ﬂow slows the ﬂow detaching from the airfoil, causing an increase in airfoil lift. To take into account radial ﬂow along a rotating blade in scientiﬁc literature, many authors modify the CL distribution [29,31,32].

The curve ﬁt can be applied to any airfoil in the same aerodynamic region (the fully stalled one), because in this region the beneﬁts due to radial ﬂow are greater [3].

To increase CL distribution in the fully stalled region, a new approach was implemented. As shown in Fig. 2,a ﬁfth-order log- arithmic polynomial ðCL ¼ Pifcosti*½lnðaÞ igÞ was adopted for the ‘attached ﬂow region’ ( 2 a 18 ), while for the ‘dynamic stall region’ (18 a 90 ) the function CL ¼ 2CLmax*sin(a)*cos(a)w as adopted. This last function intersects the logarithmic polynomial curve at about 1/2e1/3 of its descendent part. This method furnishes the correct amount of increase in CL in the fully stalled region, and permits the 1-D numerical code, to take into account radial ﬂow along the blades.

Analogously to lift coefﬁcient, two different mathematical functions have been implemented to describe the drag coefﬁcient.

A ﬁfth-order logarithmic polynomial ðCD ¼ Pifcosti*½lnðaÞ igÞ was adopted for the ‘attached ﬂow region’ ( 2 a 18 ), while for

*sin2(a) was adopted. The “cos ti”, in the CD logarithmic polynomial, have been evaluated through the least squares method, starting from the CD experimental data [21]. Also the CDmax has been obtained from CD experimental data.

Nomenclature

A rotor area [m2] Re Reynolds number [e] a angle of attack [ ] f incoming ﬂow direction angle [ ] u angular velocity [s 1] a axial induction factor [-] a0 tangential induction factor [e] r local blade radius [m]

V0 wind velocity far upstream [m/s]

V1 relative airfoil velocity [m/s] dL lift [N] dD drag [N] dR resultant force from lift and drag [N] dN normal rotor force [N] dT tangential rotor force [N] c airfoil chord [m] r air density [kg/m3]

CL airfoil lift coefﬁcient [e] CD airfoil drag coefﬁcient [e] CLmax CL at a ¼ 45 .

Nb number of blades [e] F tip loss factor [e]

CN normal force coefﬁcient [e] lr local tip speed ratio [e] C torque [Nm]

P mechanical power [W]

Fig. 1. Forces acting on the airfoil.

2.2. Brake state model

A brake state model is a set of mathematical equations implemented within a 1-D numerical code, based on the Blade Element Momentum theory, to design and evaluate the performance of horizontal axis wind turbines.

The brake state model implements different mathematical expressions to evaluate the tangential (a0) and axial (a) induction factors. In the numerical code presented in this paper, Eqs. (5)e(7) were implemented [13].

The numerical stability of the mathematical code depends on tangential (a0) and axial (a) induction factors. Before selecting these mathematical expressions, many simulations have been carried out, implementing different mathematical expressions for the tangential and axial induction factors. In all the simulations the results were not good as those presented in this paper (see Fig. 3), and in some cases, the numerical code does not converge to the solution, but diverges to an inﬁnite loop of calculations.

In Eqs. (5) and (6) the two mathematical expressions implemented in this code for the evaluation of the axial induction factor are reported: for a < 0.4:

cNb

2pr ðCLcosf þ CDsinfÞ

In Eq. (7) the mathematical expression implemented in this code for the evaluation of the tangential induction factor is reported [21]:

where F is the Prandtl tip loss factor, as reported in [29,31].

The proposed post-stall model, in conjunction with the brake state model, has been validate through the comparison between the simulated and experimental data on the mechanical power for the NREL wind turbine (see Fig. 3). The results of the numerical code proposed in this work (including the post-stall model described in Subsection 2.1, and the brake state model described in Subsection 2.2), are very close to experimental data.

3. Comparison between experimental and simulated data

The numerical code produced by the authors, is implemented to predict the power curve of the NREL wind rotor [2]. The radius of

Fig. 2. Graphic visualization of the post-stall model.

Wind Speed [m/s]

NREL wind rotor

Experimental dataSimulated data Fig. 3. Experimental mechanical power for the NREL wind turbine.

the rotor is 5.03 m and rotates at 72 r/min. The blade section is the S809 airfoil. The experiments were carried out in the world’s biggest wind tunnel at NASA Ames.

The rotor blades are twisted and tapered. Power control is passive and occurs by deep stalling a section of the wind turbine blades. One method to maintain almost constant the wind turbine power, while the wind speed varies, is that to design the blades so that they work in the deep stall region, and power production is limited by these aerodynamic conditions (see Fig. 3, for the wind speed varying from 10 m/s to 20 m/s).

Experimental measurements of power as wind speed varies were taken from scientiﬁc literature [29].

In Fig. 3 the comparison between the experimental and simulated data is shown. It is possible to notice how the 1-D numerical code proposed in this work (with the post-stall model, and the brake state model described in Subsections 2.1 and 2.2) furnishes reliable results.

4. Numerical simulation and comparison of brake state models

The combined JonkmaneBuhl [20,21] brake state model is implemented within the numerical code developed by the authors (see Eqs. (5)e(7)), and is compared with Shen’s brake state model [18,19].

with

Y2 ¼ 4Fsinfcosf sF1CT

Figs. (4) and (5) show the brake state models. In [34], Glauert reported the experimental results showing that the trust coefﬁcient equation CN ¼ 4a(1 a) is not valid if the axial induction factor exceeds 0.4. Glauert [34] gave a correction for determining the axial induction factor, when a > 0.4, valid only for F ¼ 1. If the losses at the tip of the blade are taken into account (F < 1), the correction proposed by Buhl [20] must be considered and implemented. Fig. 4 shows Buhl’s brake state model in conjunction with experimental data taken from [23].

This correction is needed to eliminate the numerical instability which occurs when the Glauert correction is implemented in conjunction with tip losses.

In 2005, Shen et al. [18] proposed a new tip loss correction model to predict physical behavior close to the tip. The local thrust coefﬁcient is replaced by a linear relation when the axial induction value is greater than a critical value (a 1/3). Fig. 5 shows Shen’s brake state model in conjunction with experimental data taken from [23].

Implementing the numerical code, the BEM computation is carried out using 20 blade elements distributed uniformlyalong the blade. Comparative axial induction factors are ﬁrst computed. In Fig. 6, the axial induction factor is plotted as a function of radius.

The axial induction values are almost identical for the inner part of the blade but diverge when approaching the tip. This value is lower than the value obtained with the Glauert model (reported also in [18]), but greater than the value obtained with Shen’s brake state model. A greater a value implies less predicted power.

Fig. 7showsthepredictedpowercurves fortheNRELwindrotor.

It shows the predicted power curve of the simpliﬁed Glauert model [18], the predicted power curve of Shen’s model [18], and the

Fig. 4. Buhl’s brake state model.

Fig. 5. Shen’s brake state model.

predicted power curve obtained with the code developed in this work. All the predicted power curves are compared with experimental data.

Notice how in Fig. 7, the power predicted in this work is very close to the experimental data for wind speeds varying from 5 to 20 m/s. For wind speeds greater than 20 m/s, Shen’s BSM predicts a better curve.

In future research, a new strategy to implement both the BSMs presented here, will be taken into account. Each BSM will be implemented for different wind speed ranges to maximize the correlation between experimental and simulated data.

5. Conclusions

The authors produced a numerical code based on BEM theory in conjunction with an aerodynamic post-stall model, indispensable for taking into account radial ﬂow along the wind turbine blades, and the brake state models by Buhl combined with Jonkman’s tangential induction factor.

This brake state model was compared with that of Shen et al. to predict the power curve for an NREL wind rotor for which experimental mechanical power measurements are reported in scientiﬁc literature.

The comparison highlighted two different behaviors for the two brake state models which in this work better predict the power curves at low and middle wind speeds, whereas Shen’s works better at high wind velocities.

The advantages of the developed method are those of a 1-D numerical code: very little computational weight, the possibility to effect many simulations in a very little time, the possibility to evaluate different geometrical conﬁgurations of the wind turbine in order to obtain high power coefﬁcient, maximize the Annual Energy Production.

The disadvantage of this numerical code is its precision. These disadvantages can be overcome with a ﬁnal 3-D CFD simulation.

References

[1] Glauert H. The elements of airfoil and airscrew theory. Cambridge University

Press; 1926. [2] Pinheiro Vaz JR, Pinho JT, Amarante Mesquita AL. An extension of BEM method applied to horizontal-axis wind turbine design. Renewable Energy 2011;36:1734e40. [3] McWilliam M, Crawford C. The behavior of ﬁxed point iteration and Newton-

(Parte **1** de 2)