Rotor Aerodynamic Modeling using Blade Element Momentum (BEM) Methods
Andrew Swift
Director, Texas Wind Energy Institute
Texas Tech University
Lubbock, Texas 79409
February 2011
Computer programs like PROP and the NREL equivalent, WT-PERF are rotor design and simulation tools based on the Glauert Blade Element-Momentum (BEM) Method. This method requires the division of rotor into sections, or elements, and knowledge of detailed rotor geometry and airfoil characteristics for each section. A computer algorithm is then used for solution.
The Glauert BEM method uses a set of equations that lend themselves to spread sheet calculation to provide first order rotor analysis and design. The method is named after H. Glauert, a British Aerodynamicist working in the 1920-1930’s who suggested a rotor blade could be divided into radial sections and each section analyzed independently. Momentum Theory, needed to calculate the induced velocity at the rotor, is derived using actuator disc theory from Rankine and Froude and from Newton’s Second Law, impulse-momentum formulation for steady flow in a control volume, as described in most undergraduate texts on Newtonian Dynamics. Equation development and spread sheet algorithm, as well as an example simulation, is provided in this paper.
Author’s Note: Thanks to Anant Jain for his valuable assistance in editing and preparing this document.
Table of Contents
Section 1: Nomenclature
Section 2: Definitions of non-dimensional terms
Section 3: Momentum Theory
3.1 Actuator Disc Theory
3.2 Summary of important results from Momentum Theory
Section 4: Glauert Annulus – Blade Element Theory
4.1 Assumptions and Equation Development
4.2 Summary of important results
Section 5: Instructions to Create BEM code using Microsoft Excel
5.1 Algorithm
5.2 Example 1(based on data for the GE 1.5 MW Wind Turbine)
Section 6: Airfoil Data for Modeling
Section 6.1 Lift and Drag Data
Section 6.2 Post-Stall Modeling
Section 7: BEM Corrections and Refinements
Section 7.1: Introduction
Section 7.2: Description of Modification Theories
Section 8: Comparisons to “QuickProp”
To be added
References:
Section 1: Nomenclature:
BEM Rotor Analysis Method ©Swift, A., Texas Tech University, 11/14/18Page 1
P – Power
CP – Power Coefficient
ρ – Air Density
A – Swept Area of Rotor
Aa – Area of Elemental Annulus
A* - Area of the airfoil
V – Free Stream Wind Velocity
VRel – Relative Wind Velocity
Vrotor – Velocity of the wind at the rotor
ΔV - Change in Velocity
T – Thrust
Q – Torque
CT – Coefficient of Thrust
Cn – Normal Force Coefficient
Cin-plane– In-Plane Force Coefficient
a –Induced Velocity Factorat the rotor disc
b – Induced Velocity Factor down-stream of the rotor
σ – Solidity of the Rotor
B – Number of Blades
n – Number of blade elements
r – Selected Blade Radius
Δr – Elemental Blade Section
R – Total Radius of Rotor
c – Local Chord Length at radial station r
c.67 – Chord Length at 67% or 2/3 R
ctip – Chord Length at the rotor tip,r=R
λ – Tip Speed Ratio
ω – Angular Velocity of Rotor
α –Angle of Attack
φ – Inflow Angle
θ – Pitch Angle
L – Lift Force
D – Drag Force
Cl – Lift Coefficient
Cd – Drag Coefficient
δ – Lift Curve Slope
BEM Rotor Analysis Method ©Swift, A., Texas Tech University, 11/14/18Page 1
Section 2: Definitions of non-dimensional terms:
Power Coefficient:
……….. (2.1)
Thrust Coefficient:
……….. (2.2)
Lift Coefficient:
..……… (2.3)
Drag Coefficient:
..……… (2.4)
Normal Force Coefficient:
..……… (2.5)
In-Plane Force Coefficient:
..……… (2.6)
Tip Speed Ratio:
……….. (2.7)
Rotor Solidity:
…………(2.8)
Section 3: Momentum Theory
3.1 Actuator Disc Theory
This is the simplest model for rotor analysis and forms the basis for using momentum theory to describe rotor performance based on the induced velocity. The method and resulting relationships are used in propeller and helicopter analysis as well as wind turbine design and analysis.The figure below describes the flow:
Figure 3.1: Actuator Disc model of a wind turbine of Swept Area (A)
Assumptions:
- Steady, homogeneous wind inflow
- No flow obstructions either upstream or downstream
- Uniform flow velocity at the disc
- Flow is incompressible
- No flow rotation produced by the disc
- Infinite number of rotor blades at the disc
Power is extracted at the disc and the reduced velocity at the disc causes the stream tube to expand downwind of the rotor due to the conservation of mass and the assumption ofincompressibility. The power output is the product of thrust times the velocity of the air stream at the rotor disc; Power = Thrust Force × Velocity;
....…… (3.1)
Where, using Newton’s Second Law for momentum change in steady flow,
……… (3.2)
And,
……… (3.3)
By combining eqns. 3.1 and 3.2, the expression for power output of the actuator disc in Figure 3.1 can be stated as:
…….. (3.4)
Where, a is the axial induction factor at the rotor and b is the axial induction factor downstream of the rotor.
Also, power can be definedas the rate of change of kinetic energy in the wind due to the rotor slowing the wind:
…..… (3.5)
..…… (3.6)
Equating eqns. 3.4 and 3.6, the relation between the induced velocity factors at the rotor disc and far downstream can be obtained:
..…… (3.7)
The above eqns. can be used to relate induced velocity factor with the coefficient of thrust and co-efficient of power. Using definitions for CP and CT (section 2: eqns. 2.1 and 2.2), one can derive using eqns. 3.2, 3.6 and 3.7, the following relationships from Momentum Theory:
……… (3.8)
………. (3.9)
Where:
The induction factor at the rotor must be less than 0.5 to avoid a flow reversal and a turbulent wake, where momentum theory is not applicable. The Betz limit comes from eqn. 3.8 for a = 1/3; Cp = 16/27 or 59%. The rotor induction factor for a given rotor operating state is generally unknown and must be determined, as shown in Section 4: Glauert Annulus – Blade Element Theory.
3.2 Summary of important results from Momentum Theory:
1)Relation between a and b (induction factors):
2)Power Coefficient:
3)Thrust Coefficient:
Section 4: Glauert Annulus – Blade Element Theory
4.1Assumptions and Equation Development
Using the actuator disc model, Glauert assumed that the rotor could be analyzed by looking at each radial section of the blade independently, and that the sections did not interact with one another in axial flow
Figure 4.1: Radial sections of a rotor
A typical blade element is shown in Figures 4.2 a and b operating at radial station r with the related velocity vectors shown. V(1-a) is the wind velocity, slowed by the axial induction factor a at the rotor, which is unknown and must be calculated, and ωris the in-plane wind velocity due to the rotation of the rotor.
Figure 4.2(a) 3-Dimensional View of Rotor Blade Element Δrat radial station r.
Figure 4.2 (b) Section view of blade element Δr operating at radial station r
From the blade element diagram, Fig 4.2 b one can write:
……… (4.1.a)
...……. (4.1.b)
Applying the sine trigonometric principle to angle φ(Fig. 4.2 b), one has
……. (4.2.a)
The definitions of Lift and Drag forces are given by (refer section 2: eqns. 2.3 and 2.4):
……. (4.2.b)
……. (4.2.c)
Where, A* is the area of the rotor blade element Δr and;
……. (4.2.d)
Referring to Fig 4.2 b again, one can write the thrust force for a blade element as the sum of the normal components of lift and drag forces (along the direction of the thrust force), which is given by:
……. (4.2.e)
Inserting eqns. 4.2.b and 4.2.c into eqn. 4.2.e, one has:
……. (4.2.f)
The total rotor thrust will be the incremental thrust force for ALL the blades. Therefore one has to multiply eqn. 4.2.fbyB (number of blades on the rotor). Also, applying eqn. 4.2.d to eqn. 4.2.f, the equation for incremental thrust becomes:
...……. (4.3)
Where Cn is the axial force coefficient, and is given by:
……… (4.4)
Next, the incremental thrust forces for each blade element(eqn. 4.3) are equated using momentum theoryto each annular element of the rotor swept area:
……… (4.5.a)
From Momentum Theory Eqn. 3.7, one can set b = 2a and solve.
The area Aa (rotor annular ring area) on the right hand side of the equation is different from A* (airfoil area) and A (rotor swept area). The rotor annular ring area Aa is given by:
Aa=2πr•Δr……… (4.5.b)
Inserting eqns. (3.7), (4.2.a) and (4.5.b) into eqn. (4.5.a), we get:
……… (4.6.a)
Which can be solved for a and rewritten as:
……..(4.6.b)
Next, we derive the eqn. for calculating the torque applied on the rotor. Refer again to Fig. 4.2 b. The torque exerted on the rotor by blade element Δr at radial station r is:
……… (4.7.a)
..……. (4.7.b)
Then, following the incremental thrust equation development procedure, one has:
.…….. (4.7.c)
Where, CIN-PLANE is the tangential force coefficient acting on the rotor, and is given by the relation:
………. (4.8)
We know that:
ΔP = ωΔQ ...……. (4.9)
Therefore using eqn. 4.7.c and eqn. 4.9 the incremental power at each element of the blade, is given by:
..…….. (4.10.a)
…...…. (4.10.b)
Where, n = number of blade elements
Using, eqn. 4.10.b and eqn. 2.1, the power coefficient is given by:
…..……. (4.11)
4.2 Summary of important results:
1)Induction Factor:
; For each blade element …(4.6.a)
2)Thrust acting on the rotor:
....…. (4.12.a)
……. (4.12.b)
3)Torque exerted on the rotor:
……. (4.13.a)
……. (4.13.b)
4)Power produced by the rotor:
……. (4.14.a)
……. (4.14.b)
Section 5: Instructions to Create BEMcode using Microsoft Excel
Note 1: This methodology was adapted for wind turbines by Wilson and Lissaman, and called the PROP code, Ref. [4].
5.1 Algorithm
The steps to create a BEM code ( or PROP code) are outlined below:
1)Divide the rotor blade into n elements, usually 10 and calculate Δr.
2)Calculate the radial station, r, for the center of each element (Usually the blades do not start at the center of rotation and an adjustment must be made).
3)Start with the first blade element from center of the rotor.
4)Assume a value for the induction factor a,usually 0 or 1/3
5)Calculate the inflow angle φ.
6)Compute angle of attack α.
7)Import the numerical values of coefficients of lift and drag (CL and CD) from a reference table or appropriate equations.
8)Calculate the value of new induction factor anew.
9)Check with the assumed value and iterate to convergence (MS Excel solver add-in).
10)Calculate the values of incremental thrust (ΔT) and incremental torque (ΔQ).
11)Repeat steps (3) to (10) for all the elements of the blade.
12)Sum up the values of incremental thrust (ΔT) and incremental torque (ΔQ) to calculate T and Q.
13)Calculate the thrust coefficient
14)Calculate total power P and compute power coefficient Cp.
Details of how to apply this algorithm are presented in the example that follows.
5.2 Example 1(based on data for the GE 1.5 MW Wind Turbine)
Specifications:
- Turbine rotor diameter: 77 m
- 3 Blades
- Chord = 1.5 meters (assumed constant for this analysis)
- Pitch angle: 0 degrees (assumed constant for this analysis)
- Rated power: 1.5 MW
- Rotor speed: 21 RPM (assumed constant for this analysis)
Solution (Usually use SI units so that power is calculated directly in Watts):
1)Open a new workbook in Microsoft ExcelSheet and save the file as BEM.Basic.
2)Create Table 1 consisting following parameters with their numerical values.
- Wind Velocity: 10 m/s
- Rotor Angular Velocity: 2.2 rad/s (21 rpm)
- Tip Speed Ratio: = 8.47
- Number of Blades: 3
- Rotor Radius: 38.5 m
- Number of Sections: 10
- Chord: 1.5 m (constant)
- Pitch: 0(rad or deg.)(constant)
- Sea Level Air Density: 1.23 kg/m3
Tabular Form:
Wind Velocity (m/s) / Rotor Angular Velocity (rad/s) / Tip Speed Ratio / No. of Blades / Rotor Radius (m) / No of Sections / Chord (m) / Pitch (rad or deg.) / Air Density (kg/ m^3)10 / 2.2 / Calculate using Eqn. 2.7 = 8.47 / 3 / 38.5 / 10 / 1.5 / 0 / 1.23
Table 1: Independent Parameters
3)Create Table 2 (on the same sheet as above) for data of the following aerodynamic parameters
Note 2: The BEM equations developed previously are renumbered for convenience and the original equation number is written in “[ ]” brackets next to the new equation number.
- Initial Induction Factor (a): a = 0.333
- Section Number (SN): usually 10 or more
- Section radial station (r):
..…….. (5.2)
- Inflow Angle (φ):
……… (5.3) / [4.1.b]
- Relative Velocity (Vrel):
……. (5.4) / [4.2.a]
- Angle of Attack (α):
……… (5.5) / [4.1.a]
Note 3: Before moving onto the next step, please import the numerical values of coefficients of lift and drag (CL and CD) from a reference table (provided for in-class example and relevant assignments).
- Normal Force Coefficient (CN):
..……… (5.6) / [4.4]
- In-plane Force Coefficient (CIN-PLANE):
………. (5.7) / [4.8]
- New Induction Factor (anew):
…….. (5.8) / [4.6.b]
- Error between a and anew:
………. (5.9)
Note 4: Use the add-in option in MS EXCEL to add SOLVER feature under DATA tab. After installation, use SOLVER feature to drive the error between a and anewto a value near zero.
- Incremental Thrust (ΔT):
...……. (5.10) / [4.3]
- Incremental Torque (ΔQ):
..…….. (5.11) / [4.7.c]
Tabular Form:
Ind.Vel.
Factor (a) / Sect. No. (SN) / Sect. Radius r
(m) / Inflow Angle
φ
(rad or deg.) / Rel. Vel. Vrel (m/s) / Angle Attackα
(rad or deg.) / CN / CIN-PL / New Ind. Factor (anew) / Error bet.
a and anew / Thrust
T
(N) / TorqueQ (N-m)
.33 / 1 / Eqn.5.2 / Eqn.5.3 / Eqn.5.4 / Eqn.5.5 / Eqn.5.6 / Eqn.5.7 / Eqn.5.8 / Eqn.5.9 / Eqn.5.11 / Eqn.5.12
.33 / 2…to / Eqn.5.2 / Eqn.5.3 / Eqn.5.4 / Eqn.5.5 / Eqn.5.6 / Eqn.5.7 / Eqn.5.8 / Eqn.5.9 / Eqn.5.11 / Eqn.5.12
.33 / 10 / Eqn.5.2 / Eqn.5.3 / Eqn.5.4 / Eqn.5.5 / Eqn.5.6 / Eqn.5.7 / Eqn.5.8 / Eqn.5.9 / Eqn.5.11 / Eqn.5.12
Table 2: Dependent Parameters
4)Create Table 3 (usually on the same sheet as above) to find numerical values of Power, Power Coefficient and Thrust Coefficient.
- Total Thrust: Sum up all the values of the eleventh column of Table 2
……… (5.12) / [4.12.b]
- Total Torque: Sum up all the values of the twelfth column of Table 2
………. (5.13) / [4.13.b]
- Power:
………. (5.14) / [4.9]
- Power Coefficient:
……….. (5.15) / [2.1]
- Thrust Coefficient:
……….. (5.16) / [2.2]
Tabular Form using SI units:
∑ Thrust (N) / ∑ Torque(Nm) / Power (W) / CP / CTEqn. (5.12) / [4.12.b] / Eqn. (5.13) / [4.13.b] / Eqn. (5.14) / [4.9] / Eqn. (5.15) / [2.1] / Eqn. (5.16) / [2.2]
Table 3: Performance Parameters
Note 5: All the steps shown above can be repeated choosing a new wind speedto obtain the Power-Curve or Cp vs Tip Speed Ratio curve for incremental wind velocities and related Tip Speed Ratios.
Section 6: Airfoil Data for Modeling
Section 6.1 Lift and Drag Data
In order to model rotor airfoils, lift and drag data are required. Usually the data is presented as a table look-up as a function of angle of attack and Reynolds number. The data usually comes from computer codes or airfoil tests in a wind tunnel.
In this work a nominal airfoil is presented to substitute for actual airfoil data. Figure 7.1 below shows the lift and drag chart for the nominal airfoil.
Pre-stall behavior is modeled analytically with the following equations.
Nominal Airfoil Pre-Stall
CL = 0.9*2π*α
Cd = 0.01 + 0.5 α2
Post-stall behavior requires a model, as described next.
Section 6.2Post-Stall Modeling
What follows is a description of the Viterna-Corrigan post stall model. It is based on an empirical analysis of airfoil performance beyond stall specifically for wind turbine applications.
Most airfoil data stops at values beyond stall, since aircraft typically do not operate in that region. Wind turbines, especially fixed pitch, stall controlled turbines depend on the post stall behavior of the turbine, making this data necessary for analysis.
[The following text is from:
Eggleson and Stoddard, “Wind Turbine Engineering Design”, (1987); Pgs. 60-61]
For a Nominal Airfoil the following approximate data apply:
Pre-Stall:CL = 0.9*2π*α ; Cd = 0.01 + 0.5 α2
Stall: αstall = αs = 15 degrees = 0.26 radians
Post-Stall: Use equations 2.106, 107 and 108 above where:
Cd (stall)= Cds = .01+(0.26)2 = 0.044
CL (stall) = Cls = .9 * 2π*(0.26) = 1.48
SECTION 7: BEM CORRECTIONS AND REFINEMENTS
7.1 Introduction
The BEM theory discussed above has certain limitations towards calculation of accurate results for wind turbine performance simulation.For example, the assumption of no rotation in the air flow is not valid in a physical model. To lift this assumption we need to introduce a tangential induction factor, which is described later in this section. The assumption of actuator disc theory assumes a rotor disc with infinite number of blades, which is contradictory to an actual turbine rotor. In order to reduce the inaccuracy, the Prandtl tip loss and hub loss models can be applied to the formulations described in the earlier sections in this report. Furthermore, the BEM theory does not remain valid when the axial induction factor goes beyond 0.5 and the turbine enters the “turbulent wake state”. Figure7.1(below) shows empirical values for CTand that values of thrust coefficient begin to diverge from momentum theory at values above a = 0.4.
Figure 7.1: Glauert Correction for BEM theory for different tip loss factor values.
7.2 Descriptions of Modifications for BEM models
The brief descriptions of the theories that refine and improve BEM theory and modeling are given below:
1)Tangential Induction Factor: The tangential induction factor is introduced to add the wake rotation characteristics into the modeling algorithm, which can be defined as the normalized measure of the induced tangential velocity that acts on the rotor in the direction opposite to its rotation. The figure below shows the modified tangential velocity component.
Figure 7.2: Modified Tangential Velocity Component
The tangential induction factor is given by the following relation
………. (7.1)
Due to the addition of tangential induction factor to our algorithm, the formulation for some dependent parameters (Section 4: Table 2 in MS Excel Sheet) have to be changed. Apart from these changes one also is required to add a separate column for tangential induction factor next to the axial induction factor and perform the same iterative operations on the tangential induction factor as conducted to calculate the axial induction factor using incremental torque values in a similar fashion; i.e. assuming a’ = 0.33 and then recalculating every parameter which is influenced by the introduction of tangential induction factor. The following parameters are modified and are described below with their corresponding old and new formulation:
- Inflow Angle: The inflow angle can now be expressed in terms of tangential induction factor. Previously, it was expressed in terms of axial induction factor only. The following expressions for the inflow angle can be written in reference to Figure 7.1:
………. (4.1.b)
………. (7.1.a)
- Relative Velocity: Again referring to Figure 7.1, one can write the trigonometric sine and cosine relationships for relative velocity as:
……. (4.2.a)