Blade element momentum theory

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 BEM 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,1889). This article explains the application of the "Betz limit" to the efficiency of a ground-based wind turbine.

A development came in the form of Froude's blade element momentum theory (1878), later refined by Glauert (1926). Betz (1921) provided an approximate correction to momentum "Rankine–Froude actuator-disk" theory 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.^[1] 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:

$\frac{1}{2}\rho v^2 + P + \rho g h = \mathrm{const.}$

where $\rho$ is the density of the fluid, $v$ is the velocity of the fluid along a streamline, $P$ is the static pressure energy, $g$ is the acceleration due to gravity, and $h$ 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:

$\frac{1}{2}\rho v^2 + P = \mathrm{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_1$ and the pressure at 1 is $P_1$, and at point 2 the velocity of the fluid along the streamline is $v_2$ and the pressure at 2 is $P_2$, and no energy has been extracted from the fluid between points 1 and 2, then we have the following expression:

$\frac{1}{2}\rho v_1^2 + P_1 = \frac{1}{2}\rho v_2^2 + P_2$

Now let us return to our initial diagram. Consider pre-actuator flow. Far upstream, the fluid velocity is $v_{\infty}$; the fluid then expands as it approaches the rotor.^[2] In accordance with mass conservation, the mass flow rate must be constant. The mass flow rate, $\dot{m}$, through a surface of area $A$ is given by the following expression:

$\dot{m} = \rho Av$

where $\rho$ is the density and $v$ 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. Why a streamtube expands as it approaches an object is not explained in this document.

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:

$\frac{1}{2}\rho v_{\infty}^2 + P_{\infty} = \frac{1}{2}\rho \left(v_{\infty}(1 - a)\right)^2 + P_{D+}$

where we have written the fluid velocity at the rotor as $v_{\infty}(1 - a)$, where $a$ 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 let us 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 i.e. $P \rightarrow P_{\infty}$ far downstream. Assuming no further energy transfer, we can apply Bernoulli for downstream:

$\frac{1}{2}\rho \left(v_{\infty}(1 - a)\right)^2 + P_{D-} = \frac{1}{2}\rho v_w^2 + P_{\infty}$

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

$P_{D+} - P_{D-} = \frac{1}{2}\rho(v_{\infty}^2 - v_w^2)$

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 = \Delta PA$:

$\frac{1}{2}\rho(v_{\infty}^2 - v_w^2)A_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:

$F = \frac{\mathrm{d}p}{\mathrm{d}t} = \dot{m}(v_{\infty} - v_{w}) = \rho A_Dv_D(v_{\infty} - v_{w}) = \rho A_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:

$\mathrm{Power}_{ext} = Fv_D = 2a(1 - a)^2v_{\infty}^3\rho 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:

$\mathrm{Power} = \frac{1}{2}\rho A v^3$

where $\rho$ is the fluid density as before, $v$ is the fluid velocity, and $A$ 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, $\mathrm{Power}_{ext}$ is given by the following expression:

$\mathrm{Power}_{ext} = \mathrm{Power}\times 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 ($\mathrm{Power}_{ext}$). We can see that the power extracted is dependent on the axial induction factor. If we differentiate $\mathrm{Power}_{ext}$ with respect to $a$, we get the following result:

$\frac{\mathrm{d}\mathrm{Power}_{ext}}{\mathrm{d}a} = 2v_{\infty}^3\rho A_D \times \left((1 - a)^2 - 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 $a$ which yields maximum power extraction. This value is a third. 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 $a$, we will relate the change in angular momentum of the fluid with the tangential induction factor, $a'$.

Let us consider the following setup.^[1]

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:

$\frac{1}{2}\rho v_1^2 + P_1 = \frac{1}{2}\rho v_2^2 + P_2$

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^2$.

NOTE: We will in fact, be working in cylindrical co-ordinates for all aspects e.g. $\mathbf{F} = F_r\hat{\mathbf{r}} + F_{\theta}\hat{\theta} + F_z\hat{\mathbf{z}}$

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

Pre-rotor

$P_{\infty} + \frac{1}{2}\rho v_u^2 = P_{D+} + \frac{1}{2}\rho v_D^2$

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} + \frac{1}{2}\rho v_{\infty}^2 = P_{D+} + \frac{1}{2}\rho (v_{\infty}(1 - a))^2$

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.

It should be noted that, 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 \neq 0$ upstream. However, for the purpose of this analysis, that effect will be neglected.

Post-rotor

$P_{D-} + \frac{1}{2}\rho v_D^2 = P_{\infty} + \frac{1}{2}\rho v_w^2$

where $v_D$ is the velocity of the fluid just after interacting with the rotor. This can be written as $v_D^2 = v_{D,~r}^2 + v_{D,~\theta}^2 + v_{D,~z}^2$. 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-} + \frac{1}{2}\rho v_{D,~z}^2 = P_{\infty} + \frac{1}{2}\rho v_{w,~z}^2 = P_{D-} + \frac{1}{2}\rho (v_{\infty}(1 - a))^2$

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, z} = (1 - 2a)v_{\infty}$

$\mathrm{Power} = 2a(1 - a)^2v_{\infty}^3\rho 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:

$\delta\mathbf{Q} = 2\pi r\delta r \times \rho U_{\infty}(1 - a) \times 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'\omega r$; and the rotor motion itself, $\omega r$. That is, the apparent fluid velocity is given as below:

$\mathbf{v} = \omega r(1 + a') \hat{\mathbf{\theta}} + v_{\infty}(1 - a)\hat{\mathbf{z}}$

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

$|\mathbf{v}|^2 = (\omega r (1 + a'))^2 + (v_{\infty}(1 - a))^2 = W^2$

We can also work out the angle $\phi$ from the above figure:

$\sin\phi = \frac{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 $\delta r$ (see figure below).

The lift acting on those parts of the blades/airfoils each with chord $c$ is given by the following expression:

$\delta L = \frac{1}{2}\rho NW^2 c \times c_L(\alpha)\delta r$

where $c_L$ is the lift co-efficient, which is a function of the angle of attack, and $N$ is the number of blades. Additionally, the drag acting on that part of the blades/airfoils with chord $c$ is given by the following expression:

$\delta D = \frac{1}{2}\rho NW^2 c \times 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{\mathbf{z}}$ and $\hat{\theta}$ axes. Thus we need to consider the diagram below:

Thus we can see the following:

$\delta F_{\theta} = \delta L\sin\phi - \delta D\cos\phi$

$\delta F_z = \delta L\cos\phi + \delta D\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'\omega r$. Thus the torque in the air is given by

$|\mathbf{\delta{Q}}| = \rho(2\pi r\delta r)U_{\infty}(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,

$\frac{1}{2}\rho W^2Nc(c_l\sin\phi - c_d\cos\phi)r\delta r = \rho(2\pi r\delta r)U_{\infty}(1 - a)\times(2\Omega a'r^2)$

$\frac{1}{2} W^2Nc(c_l\sin\phi - c_d\cos\phi) = 4\pi U_{\infty}(1 - a)\times\Omega a'r^2$

Furthermore, the rate of change of linear momentum in the air is balanced by the out-of-plane bending force acting on the blades, $\delta F_z$. From momentum theory, the rate of change of linear momentum in the air is as follows:

$\delta F_z = \rho(2\pi r\delta r)U_{\infty}(1 - a)\times(v_\infty - v_w)$

which may be expressed as

$\delta F_z = \rho(4\pi r\delta r)U^2_{\infty}a(1 - a)$

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

$\frac{1}{2} W^2Nc(c_l\cos\phi + c_d\sin\phi) = \rho(4\pi r\delta r)U^2_{\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:

$\frac{1}{2} W^2NcC_y = 4\pi U_{\infty}(1 - a)\times\Omega a'r^2$

(1)

$\frac{1}{2}\rho W^2NcC_x = 4\pi\rho\left[(a'\Omega r)^2 + U^2_{\infty}a(1 - a)\right]r$

(2)

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

$\sin\phi = \frac{U_{\infty}}{W}(1 - a) \rightarrow \sin^2\phi = \left(\frac{U_{\infty}}{W}(1-a)\right)^2$

(3)

Thus, with these three equations, it is possible to get the following result through some algebraic manipulation:^[1][2]

$\frac{a}{1 - a} = \frac{C_x\sigma_r}{4\sin^2\phi}$

We can derive an expression for $a'$ 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^[3]
• Does not account for wake expansion
• Does not account for tip losses, though correction factors can be included.^[2]
• Does not account for yaw, though it can be made to do so
• Based on steady flow (non-turbulent)

References

1. Wind Energy Handbook: Burton, Jenkins
2. http://cdn.intechopen.com/pdfs/16241/InTech-Aerodynamics_of_wind_turbines.pdf
3. http://www.stanford.edu/~eena/windpower07.pdf