Blade element momentum theory explained

Blade element momentum theory is a theory that combines both blade element theory and momentum theory. It is used to calculate the local forces on a propeller or wind-turbine blade. Blade element theory is combined with momentum theory to alleviate some of the difficulties in calculating the induced velocities at the rotor.

This article emphasizes application of blade element theory to ground-based wind turbines, but the principles apply as well to propellers. Whereas the streamtube area is reduced by a propeller, it is expanded by a wind turbine. For either application, a highly simplified but useful approximation is the Rankine–Froude "momentum" or "actuator disk" model (1865,^[1] 1889^[2]). This article explains the application of the "Betz limit" to the efficiency of a ground-based wind turbine.

Froude's blade element theory (1878)^[3] is a mathematical process to determine the behavior of propellers, later refined by Glauert (1926). Betz (1921) provided an approximate correction to momentum "Rankine–Froude actuator-disk" theory ^[4] to account for the sudden rotation imparted to the flow by the actuator disk (NACA TN 83, "The Theory of the Screw Propeller" and NACA TM 491, "Propeller Problems"). In blade element momentum theory, angular momentum is included in the model, meaning that the wake (the air after interaction with the rotor) has angular momentum. That is, the air begins to rotate about the z-axis immediately upon interaction with the rotor (see diagram below). Angular momentum must be taken into account since the rotor, which is the device that extracts the energy from the wind, is rotating as a result of the interaction with the wind.

Rankine–Froude model

The "Betz limit," not yet taking advantage of Betz' contribution to account for rotational flow with emphasis on propellers, applies the Rankine–Froude "actuator disk" theory to obtain the maximum efficiency of a stationary wind turbine. The following analysis is restricted to axial motion of the air:

In our streamtube we have fluid flowing from left to right, and an actuator disk that represents the rotor. We will assume that the rotor is infinitesimally thin.^[5] From above, we can see that at the start of the streamtube, fluid flow is normal to the actuator disk. The fluid interacts with the rotor, thus transferring energy from the fluid to the rotor. The fluid then continues to flow downstream. Thus we can break our system/streamtube into two sections: pre-acuator disk, and post-actuator disk. Before interaction with the rotor, the total energy in the fluid is constant. Furthermore, after interacting with the rotor, the total energy in the fluid is constant.

Bernoulli's equation describes the different forms of energy that are present in fluid flow where the net energy is constant, i.e. when a fluid is not transferring any energy to some other entity such as a rotor. The energy consists of static pressure, gravitational potential energy, and kinetic energy. Mathematically, we have the following expression:

	1
	2

\rhov²+P+\rhogh=const.

where

\rho

is the density of the fluid,

is the velocity of the fluid along a streamline,

is the static pressure energy,

is the acceleration due to gravity, and

is the height above the ground. For the purposes of this analysis, we will assume that gravitational potential energy is unchanging during fluid flow from left to right such that we have the following:

	1
	2

\rhov²+P=const.

Thus, if we have two points on a streamline, point 1 and point 2, and at point 1 the velocity of the fluid along the streamline is

v₁

and the pressure at 1 is

P₁

, and at point 2 the velocity of the fluid along the streamline is

v₂

and the pressure at 2 is

P₂

, and no energy has been extracted from the fluid between points 1 and 2, then we have the following expression:

	1
	2

\rho

	2
v
	1

+P₁=

	1
	2

\rho

	2
v
	2

+P₂

Now let us return to our initial diagram. Consider pre-actuator flow. Far upstream, the fluid velocity is

v_infty

; the fluid velocity then decreases and pressure increases as it approaches the rotor. In accordance with mass conservation, the mass flow rate through the rotor must be constant. The mass flow rate,

•

, through a surface of area

is given by the following expression:

•

=\rhoAv

where

\rho

is the density and

is the velocity of the fluid along a streamline. Thus, if mass flow rate is constant, increases in area must result in decreases in fluid velocity along a streamline. This means the kinetic energy of the fluid is decreasing. If the flow is expanding but not transferring energy, then Bernoulli applies. Thus the reduction in kinetic energy is countered by an increase in static pressure energy.

So we have the following situation pre-rotor: far upstream, fluid pressure is the same as atmospheric,

P_infty

; just before interaction with the rotor, fluid pressure has increased and so kinetic energy has decreased. This can be described mathematically using Bernoulli's equation:

	1
	2

\rho

	2
v
	infty

+P_infty=

	1
	2

\rho\left(v_infty(1-a)\right)²+P_D+

where we have written the fluid velocity at the rotor as

v_infty(1-a)

, where

is the axial induction factor. The pressure of the fluid on the upstream side of the actuator disk is

P_D+

. We are treating the rotor as an actuator disk that is infinitely thin. Thus we will assume no change in fluid velocity across the actuator disk. Since energy has been extracted from the fluid, the pressure must have decreased.

Now consider post-rotor: immediately after interacting with the rotor, the fluid velocity is still

v_infty(1-a)

, but the pressure has dropped to a value

P_D-

; far downstream, pressure of the fluid has reached equilibrium with the atmosphere; this has been accomplished in the natural and dynamically slow process of decreasing the velocity of flow in the stream tube in order to maintain dynamic equilibrium (i.e.

P → P_infty

far downstream. Assuming no further energy transfer, we can apply Bernoulli for downstream:

	1
	2

\rho\left(v_infty(1-a)\right)²+P_D-=

	1
	2

\rho

	2
v
	w

+P_infty

where

v_w=

The velocity far downstream in the Wake

Thus we can obtain an expression for the pressure difference between fore and aft of the rotor:

P_D+-P_D-=

	1
	2

	2
\rho(v
	infty

	2)
v
	w

If we have a pressure difference across the area of the actuator disc, there is a force acting on the actuator disk, which can be determined from

F=\DeltaPA

	1
	2

	2
\rho(v
	infty

	2)A
v
	D

where

A_D

is the area of the actuator disk. If the rotor is the only thing absorbing energy from the fluid, the rate of change in axial momentum of the fluid is the force that is acting on the rotor. The rate of change of axial momentum can be expressed as the difference between the initial and final axial velocities of the fluid, multiplied by the mass flow rate:

	dp
	dt

•

(v_infty-v_w)=\rhoA_Dv_D(v_infty-v_w)=\rhoA_D(1-a)v_infty(v_infty-v_w)

Thus we can arrive at an expression for the fluid velocity far downstream:

v_w=(1-2a)v_infty

This force is acting at the rotor. The power taken from the fluid is the force acting on the fluid multiplied by the velocity of the fluid at the point of power extraction:

Power_ext=Fv_D=2a(1-

	3\rho
a)
	infty

A_D

Maximum power

Suppose we are interested in finding the maximum power that can be extracted from the fluid. The power in the fluid is given by the following expression:

Power=

	1
	2

\rhoAv³

where

\rho

is the fluid density as before,

is the fluid velocity, and

is the area of an imaginary surface through which the fluid is flowing. The power extracted from the fluid by a rotor in the scenario described above is some fraction of this power expression. We will call the fraction the power co-efficient,

C_p

. Thus the power extracted,

Power_ext

is given by the following expression:

Power_ext=Power x C_p

Our question is this: what is the maximum value of

C_p

using the Betz model?

Let us return to our derived expression for the power transferred from the fluid to the rotor (

Power_ext

). We can see that the power extracted is dependent on the axial induction factor. If we differentiate

Power_ext

with respect to

, we get the following result:

	dPower_ext
	da

	3\rho
2v
	infty

A_D x \left((1-a)²-2a(1-a)\right)

If we have maximised our power extraction, we can set the above to zero. This allows us to determine the value of

which yields maximum power extraction. This value is a

1/3

. Thus we are able to find that

C_P~max=16/27

. In other words, the rotor cannot extract more than 59 per cent of the power in the fluid.

Blade element momentum theory

Compared to the Rankine–Froude model, Blade element momentum theory accounts for the angular momentum of the rotor. Consider the left hand side of the figure below. We have a streamtube, in which there is the fluid and the rotor. We will assume that there is no interaction between the contents of the streamtube and everything outside of it. That is, we are dealing with an isolated system. In physics, isolated systems must obey conservation laws. An example of such is the conservation of angular momentum. Thus, the angular momentum within the streamtube must be conserved. Consequently, if the rotor acquires angular momentum through its interaction with the fluid, something else must acquire equal and opposite angular momentum. As already mentioned, the system consists of just the fluid and the rotor, the fluid must acquire angular momentum in the wake. As we related the change in axial momentum with some induction factor

, we will relate the change in angular momentum of the fluid with the tangential induction factor,

Consider the following setup.

We will break the rotor area up into annular rings of infinitesimally small thickness. We are doing this so that we can assume that axial induction factors and tangential induction factors are constant throughout the annular ring. An assumption of this approach is that annular rings are independent of one another i.e. there is no interaction between the fluids of neighboring annular rings.

Bernoulli for rotating wake

Let us now go back to Bernoulli:

	1
	2

\rho

	2
v
	1

+P₁=

	1
	2

\rho

	2
v
	2

+P₂

The velocity is the velocity of the fluid along a streamline. The streamline may not necessarily run parallel to a particular co-ordinate axis, such as the z-axis. Thus the velocity may consist of components in the axes that make up the co-ordinate system. For this analysis, we will use cylindrical polar co-ordinates

(r,~\theta,~z)

. Thus

v²=

	2
v
	r

	2
v
	\theta

	2
v
	z

NOTE: We will in fact, be working in cylindrical co-ordinates for all aspects e.g.

F=F_r\hat{r

} + F_\hat + F_z\hat

Now consider the setup shown above. As before, we can break the setup into two components: upstream and downstream.

Pre-rotor

P_infty+

	1
	2

\rho

	2
v
	u

=P_D++

	1
	2

\rho

	2
v
	D

where

v_u

is the velocity of the fluid along a streamline far upstream, and

v_D

is the velocity of the fluid just prior to the rotor. Written in cylindrical polar co-ordinates, we have the following expression:

P_infty+

	1
	2

\rho

	2
v
	infty

=P_D++

	1
	2

\rho(v_infty(1-a))²

where

v_infty

and

v_infty(1-a)

are the z-components of the velocity far upstream and just prior to the rotor respectively. This is exactly the same as the upstream equation from the Betz model.

As can be seen from the figure above, the flow expands as it approaches the rotor, a consequence of the increase in static pressure and the conservation of mass. This would imply that

v_r ≠ 0

upstream. However, for the purpose of this analysis, that effect will be neglected.

Post-rotor

P_D-+

	1
	2

\rho

	2
v
	D

=P_infty+

	1
	2

\rho

	2
v
	w

where

v_D

is the velocity of the fluid just after interacting with the rotor. This can be written as

	2
v
	D

	2
v
	D,~r

	2
v
	D,~\theta

	2
v
	D,~z

. The radial component of the velocity will be zero; this must be true if we are to use the annular ring approach; to assume otherwise would suggest interference between annular rings at some point downstream. Since we assume that there is no change in axial velocity across the disc,

v_D,~z=(1-a)v_infty

. Angular momentum must be conserved in an isolated system. Thus the rotation of the wake must not die away. Thus

v_\theta

in the downstream section is constant. Thus Bernoulli simplifies in the downstream section:

P_D-+

	1
	2

\rho

	2
v
	D,~z

=P_infty+

	1
	2

\rho

	2
v
	w,~z

=P_D-+

	1
	2

\rho(v_infty(1-a))²

In other words, the Bernoulli equations up and downstream of the rotor are the same as the Bernoulli expressions in the Betz model. Therefore, we can use results such as power extraction and wake speed that were derived in the Betz model i.e.

v_w,=(1-2a)v_infty

Power=2a(1-

	3\rho
a)
	infty

A_D

This allows us to calculate maximum power extraction for a system that includes a rotating wake. This can be shown to give the same value as that of the Betz model i.e. 0.59. This method involves recognising that the torque generated in the rotor is given by the following expression:

\deltaQ=2\pir\deltar x \rhoU_infty(1-a) x 2a'r\omega

with the necessary terms defined immediately below.

Blade forces

Consider fluid flow around an airfoil. The flow of the fluid around the airfoil gives rise to lift and drag forces. By definition, lift is the force that acts on the airfoil normal to the apparent fluid flow speed seen by the airfoil. Drag is the forces that acts tangential to the apparent fluid flow speed seen by the airfoil. What do we mean by an apparent speed? Consider the diagram below:

The speed seen by the rotor blade is dependent on three things: the axial velocity of the fluid,

v_infty(1-a)

; the tangential velocity of the fluid due to the acceleration round an airfoil,

a'\omegar

; and the rotor motion itself,

\omegar

. That is, the apparent fluid velocity is given as below:

v=\omegar(1+a')\hat{\theta

} + v_(1 - a)\hat

Thus the apparent wind speed is just the magnitude of this vector i.e.:

|v|²=(\omegar(1+a'))²+(v_infty(1-a))²=W²

We can also work out the angle

\phi

from the above figure:

\sin\phi=

	v_infty(1-a)
	W

Supposing we know the angle

\beta

, we can then work out

\alpha

simply by using the relation

\alpha=\phi-\beta

; we can then work out the lift co-efficient,

c_L

, and the drag co-efficient

c_D

, from which we can work out the lift and drag forces acting on the blade.

Consider the annular ring, which is partially occupied by blade elements. The length of each blade section occupying the annular ring is

\deltar

(see figure below).

is given by the following expression:

\deltaL=

	1
	2

\rhoNW²c x c_{L(\alpha)\delta}r

where

c_L

is the lift co-efficient, which is a function of the angle of attack, and

is the number of blades. Additionally, the drag acting on that part of the blades/airfoils with chord

is given by the following expression:

\deltaD=

	1
	2

\rhoNW²c x c_{D(\alpha)\delta}r

Remember that these forces calculated are normal and tangential to the apparent speed. We are interested in forces in the

\hat{z

} and

\hat{\theta}

axes. Thus we need to consider the diagram below:

Thus we can see the following:

\deltaF_\theta=\deltaL\sin\phi-\deltaD\cos\phi

\deltaF_z=\deltaL\cos\phi+\deltaD\sin\phi

F_\theta

is the force that is responsible for the rotation of the rotor blades;

F_z

is the force that is responsible for the bending of the blades.

Recall that for an isolated system the net angular momentum of the system is conserved. If the rotor acquired angular momentum, so must the fluid in the wake. Let us suppose that the fluid in the wake acquires a tangential velocity

v_\theta=2a'\omegar

. Thus the torque in the air is given by

|\delta{Q

}| = \rho(2\pi r\delta r)U_(1 - a)\times(2\Omega a'r^2)By the conservation of angular momentum, this balances the torque in the blades of the rotor; thus,

	1
	2

\rho

	2Nc(c
W
	l\sin\phi

-c_{d\cos\phi)r\delta}r=\rho(2\pir\deltar)U_infty(1-a) x (2\Omegaa'r²⁾

	1
	2

	2Nc(c
W
	l\sin\phi

-c_d\cos\phi)=4\piU_infty(1-a) x \Omegaa'r²

Furthermore, the rate of change of linear momentum in the air is balanced by the out-of-plane bending force acting on the blades,

\deltaF_z

. From momentum theory, the rate of change of linear momentum in the air is as follows:

\deltaF_z=\rho(2\pir\deltar)U_infty(1-a) x (v_infty-v_w)

which may be expressed as

\deltaF_z=\rho(4\pir\delta

	2
r)U
	infty

a(1-a)

Balancing this with the out-of-plane bending force gives

	1
	2

	2Nc(c
W
	l\cos\phi

+c_d\sin\phi)=\rho(4\pir\delta

	2
r)U
	infty

a(1-a)

Let us now make the following definitions:

C_y=c_l\sin\phi-c_d\cos\phi

C_x=c_l\cos\phi+c_d\sin\phi

So we have the following equations:

Let us make reference to the following equation which can be seen from analysis of the above figure:

Thus, with these three equations, it is possible to get the following result through some algebraic manipulation:

	a
	1-a

	C_x\sigma_r
	4\sin^2\phi

We can derive an expression for

in a similar manner. This allows us to understand what is going on with the rotor and the fluid. Equations of this sort are then solved by iterative techniques.

Assumptions and possible drawbacks of BEM models

Assumes that each annular ring is independent of every other annular ring.^[6]
Does not account for wake expansion.
Does not account for tip losses, though correction factors can be included.^[7]
Does not account for yaw, though it can be made to do so.
Based on steady flow (non-turbulent).

References

Rankine . William . William Rankine . 6 April 1865 . On the Mechanical Principals of the Action of Propellers . . 6 . 13 . Hathi Trust.
Froude . Robert . 12 April 1889 . On the Part Played in Propulsion by Differences in Fluid Pressure . . 30 . 390 . Hathi Trust.
Froude . William . 11 April 1878 . The Elementary Relation between Pitch, Slip, and Propulsive Efficiency . Inst. Naval Architects . 19 . 47 . Hathi Trust.
Wilson. Robert E.. Lissaman. Peter B.S.. 1974. Applied Aerodynamics of Wind Power Machines. NASA Sti/Recon Technical Report N. 75 . 22669 . 1974STIN...7522669W .
Wind Energy Handbook: Burton, Jenkins
http://www.stanford.edu/~eena/windpower07.pdf{{dead link|date=November 2016 |bot=InternetArchiveBot |fix-attempted=yes }}
Buhl. M.L. Jr. . 2005-08-01. New Empirical Relationship between Thrust Coefficient and Induction Factor for the Turbulent Windmill State. en. NREL/TP–500–36834, 15016819. 10.2172/15016819.