The blade element theory (BET) divides the blades of a horizontal-axis wind turbine into a contiguous stack of radial elements, typically 30–50 in number. They are assumed to behave as airfoils whose lift and drag give the element’s thrust and torque. These are balanced against the axial and angular momentum changes in the annular streamtube flowing over the element derived from momentum theory (MT). The combined blade element/momentum theory (BEMT) is the workhorse of wind turbine aerodynamics. BEMT is described in all aerodynamics text books and is widely used in the initial, multidimensional design of blades, for which the potentially more accurate computational fluid dynamics modeling is prohibitively costly in computer time, e.g., Sessarego et al. (2015).

Wind turbine performance is usually considered in terms of the power and thrust coefficients, C _P and C _T , respectively, as functions of the tip speed ratio, λ = ΩR/U ₀, where Ω is the blade angular velocity, R is the tip radius, and U ₀ is the wind speed. BEMT is generally held to provide accurate estimates of C _P and C _T in comparison to wind tunnel and field tests, as shown in Spera (1994), Wood (2011), and Schmitz (2020). It is well known, however, that the relation between the induced axial velocity at the rotor, u, and in the far-wake, u _∞ , u = (1 + u _∞ )/2 in the one-dimensional version of the momentum theory starts to break down somewhere near the tip speed ratio for maximum power, λ _opt . This is the start of the high-thrust region: as λ increases further, C _T approaches and sometimes exceeds unity, the maximum value allowed by conventional momentum theory. The high-thrust region for wind turbines ends at the runaway tip speed ratio, λ _R , where C _P = 0 and C _T = C _TR is usually maximized if λ _R is large (Note that the subscript “R” denotes runaway values, whereas the script “R” is the rotor radius). Since modern turbines have λ _opt ∼ 7 − 9, the high-thrust region normally spans high values of λ: as a rule of thumb λ _R ≈ 2λ _opt . This approximate relation holds for multiblade windmills where λ _opt ∼ 1, as shown in Figure 5 of John et al. (2023), and the experiments of Krogstad and Adaramola (2012) and Krogstad and Eriksen (2013) that we consider below, for which λ _opt ∼ 6. Furthermore, it is normal for C _T to increase monotonically in this region when λ _opt is large.

Pratumnopharat and Leung (2011) and Schmitz (2020) documented a large number of modifications that have been made to the relationship between C _T and u for the high-thrust region. It is our contention that these modifications are ad hoc, have little or no physical basis, and rely on questionable experimental results obtained a century ago, as explained in Section 2.7 of Wood (2011). Unfortunately, there have been very few subsequent experiments on the high-thrust region that provide sufficient detail for BEMT revision. One of the few that does was by Limacher et al. (2022), who measured the thrust and all three velocities in the wake of a small model turbine by particle image velocimetry, for values of λ above and below λ _R = 7.9. They proposed a revision to the conventional momentum equation that will be described later. The other measurements in this category were by Krogstad and Adaramola (2012) and Krogstad and Eriksen (2013), who provided power, thrust, and limited wake measurements for a model wind turbine with λ _opt ≈ 6 and λ _R ≈ 12. Their results will be compared to the revised BEMT developed here.

Our further contention is that BET remains valid in the high-thrust region, but the conventional angular and axial momentum equations are incomplete. The most important omission is in the latter equation, which should contain a term related to the “wake vorticity.” This term balances the BE drag on stationary turbines but increasingly balances the BE lift (and thrust) as λ increases. It also decouples the approach to the runaway from the relationship between u and u _∞ . Our revised BEMT equations require no specific relation between the two velocities.

The wake vorticity term was identified in the computational study of the flow through two-dimensional, equi-spaced cascades of airfoils by Golmirzaee and Wood (2023). Cascade elements are subjected to the same forces as wind turbine blade elements and have the same MT terms except for those due to expansion, rotation, and the vorticity shed from the blades as a consequence of the radial gradient of the bound vorticity. The wake vorticity term is introduced in the following section after the following preliminary observations on the induced velocities u and w in the circumferential direction. One of the key simplifications of MT is that the wake is characterized by circumferential averages of u and w. When the number of blades, N, is finite and λ is small, the velocities at the elements, which determine the lift and drag, may differ from their streamtube averages. This difference is usually accommodated by the use of “finite blade functions:” F _u = u/u _b and F _w = w/w _b , where the subscript “b” denotes a value at the blades. Nearly all BEMT codes for wind turbines use Prandtl’s tip loss factor, F _P to approximate F _u and F _w . Some of the inaccuracies of using F _P are documented in Wood et al. (2016) and Wood (2021), while Wood et al. (2021) reviewed the numerical methods that produce accurate estimates for F _u and F _w . Wood et al. (2016) showed that F _u , F _w → F _p and F _p → 1 as λ increases, suggesting that azimuthal variations become less important in the high-thrust region.

Our revisions to BEMT apply at all values of λ including zero, but we concentrate on the high-thrust region, where BEMT is most challenged. It is further noted that John et al. (2023) found conventional BEMT accurately predicted λ _R for a model of a low-speed, multi-bladed water pumping windmill where λ _R ≈ 2 and the runaway C _TR increased from 0.1 to 0.5 as N increased from 3 to 24. The implication is that the present revisions are not critical at low λ.

Our major aim is to establish the necessity to include the wake vorticity term in the axial momentum equation. This may not be the only change to MT necessary for accurate calculation of the high-thrust region, but we show that it is likely to be a very important one. We argue that more experimental information near the runaway is required to complete the revision of BEMT.

The rest of this paper is laid out as follows. Section 2 describes BEMT and our revision. We use a simple version of the revision in Section 3 to compute the power and thrust as the runaway is approached and compare these to wind tunnel measurements by Krogstad and Adaramola (2012) and Krogstad and Eriksen (2013) on a model three-bladed rotor. Section 4 focuses on further requirements for completing the revision of BEMT and the conclusions.

Using the definitions of the velocity at the blade element, U _rel , and the angle ϕ from Figure 1, the balance between the gradient of thrust, dT/dr, and the axial momentum flux for an N − bladed rotor at radius r is given as follows: d T d r = 1 2 N ρ U r e l 2 c C l ⁡ cos ⁡ ϕ + C d ⁡ sin ⁡ ϕ = π ρ 4 w 2 ̄ + 4 λ r w + v 2 ̄ − u 2 ̄ r , (1) where ρ is the density of air and c is the blade element chord. C _l and C _d are the lift and drag coefficients, respectively. They are functions of the angle of attack, α, which—according to Figure 1—is dependent on u _b and w _b . The term on the right involving u, w, and the radial velocity v is the axial momentum flux in the annular streamtube intersecting the N elements. It is important to note that w is the value at the blades, whereas the value immediately behind them is 2w. An overbar denotes a circumferential average over the streamtube for a quadratic quantity: overbars are not used on u or w on their own as their circumferential variation is important only in distinguishing values at the blades, which are subscripted as defined previously. u and w denote streamtube averages throughout this study. The last term, involving v 2 ̄ − u 2 ̄ , is called the “expansion” term because ∫ 0 ∞ v 2 ̄ − u 2 ̄ r d r = 0 , (2) which shows the expansion term re-distributes, rather than generates, T. Wood and Limacher (2021) showed that the expansion term gives the thrust due to the pressure acting on the expanding streamtubes upwind of the rotor. Eq. 2 ensures that v and u have the same magnitude at the rotor plane. Thus, significant flow expansion must occur when u is large.

The usual thrust equation from MT ignores the expansion term and contains terms in u instead of w on the right side of Eq. 1, which is the Kutta–Joukowsky form of the momentum flux. It, along with Eq. 2, was derived by Limacher and Wood (2021) using a control volume (CV) of radius much larger than R, with the inlet far upwind of the rotor where u = v = w = 0, and the outlet immediately behind the blades. Eqs 1, 2 hold for any amount of wake expansion at any λ. Their relation to the usual thrust equation is easily derived by ignoring circumferential variations and the expansion term. Wood and Hammam (2022) showed that the high−λ circumferentially uniform wakes of optimal rotors satisfy “helical symmetry:” pu = wr where p = (1 − u)/(w/r + λ) is the pitch of the helical vorticity. Applying these equations turns the right side of Eq. 1 (when v 2 ̄ − u 2 ̄ is ignored) into the usual form involving u: 4 π w w + λ r r → 4 π u 1 − u r , (3) which leads to a less general axial momentum equation than Eq. 1.

The equation for the torque, Q, derived from the same CV has the similar and familiar form: d Q d r = 1 2 N ρ U r e l 2 c C l ⁡ sin ⁡ ϕ − C d ⁡ cos ⁡ ϕ r = 4 π ρ w − u w ̄ r 2 . (4) For steady operation, the power extracted from the wind is the product of the rotor torque, obtained by summing Eq. 4 over all elements, and Ω: in the non-dimensional form, C _P = λC _Q , where C _Q is the standard torque coefficient.

We will assume that the BE parts of Eqs 1, 4 involving the lift and drag do not need revision. The validity of this “airfoil assumption” is difficult to assess. Simulations of the 2D cascade flow suggest it is conservative, because the lift:drag ratio for a cascade element is slightly greater than that on the corresponding airfoil, Golmirzaee and Wood (2023), but we know of no other direct assessment of the assumption for wind turbines.

We now revise the right side of Eqs 1, 4 based on the cascade simulations of Golmirzaee and Wood (2023). The revision is different from that of Limacher et al. (2022), who simplified the equations by assuming the shed vorticity was strain-free. They did not make a direct connection to blade element drag, and no consideration was given to the angular momentum equation. Their model and its relation to the current revision are considered in the following sections. For an infinite cascade of identical and equi-spaced airfoils, the BE or left sides of Eqs 1, 4 are unchanged except for the removal of N. For a cascade element, Golmirzaee and Wood (2023) show that Eq. 1 becomes 1 2 U r e l 2 c C l ⁡ cos ⁡ ϕ + C d ⁡ sin ⁡ ϕ = 2 S w 2 ̄ + ∫ − S / 2 S / 2 1 − u Ω z y d y , (5) when the density is removed. The cascade element spacing, S, is related to the circumferential distance between blade elements by S = 2πr/N, and there is no term corresponding to λrw in Eq. 1 as the cascade is stationary. Ω _z is the (transverse) vorticity, distinguished from the angular velocity by its subscript, and y is the normal co-ordinate, whose origin is the blade quarter-chord. The integral, which we call the “wake vorticity” integral or term, has no counterpart in the current form of the BEMT equation.

Eqs 1.5 and 1.6 of Liu et al. (2015) and Eq. (9.1.20) of Wu et al. (2015) give important constraints on wake vorticity: ∫ − S / 2 S / 2 Ω z d y = 0 and ∫ − S / 2 S / 2 1 − u Ω z d y = 0 . (6) The former, which relies on the applicability of the slender flow approximation Ω _z ≈ − ∂u/∂y, prevents Ω _z from contributing to the bound vorticity of the blade because it forces the area integral of Ω _z over the wake to be zero. The latter is exact for any wake independent of the validity of the slender flow approximation. It has two important consequences. First, the choice of origin for y has no effect on the wake vorticity integral. Second, it removes the integral containing Ω _z from the y−direction force balance and, by implication, from the angular momentum equation for BEMT because the integrand is (1 − u)Ω _z x. From here on, the phrase “wake vorticity term” will be used only for the integral in Eq. 5 or the corresponding integral in the axial momentum equation for BEMT. The exact y−direction cascade equivalent of Eq. 4 without the density is as follows: 1 2 U r e l 2 c C l ⁡ sin ⁡ ϕ − C d ⁡ cos ⁡ ϕ = 2 S w − u w ̄ , (7) showing that the cascade analysis retains the non-linear terms from Eqs 1, 4. Golmirzaee and Wood (2023) found that the spatial variations in the velocities caused only small differences between u w ̄ and uw or between w 2 ̄ and w ². It is reasonable, then, to assume that significant circumferential variations in u and w for BEMT arise only through the shed vorticity, as explained previously for F _u and F _w . Furthermore, if F _u and F _w are close to unity, then it is reasonable to also ignore any circumferential variation in u ², w ², and uw in the high-thrust region.

For wind turbines, the wake vorticity is largely radial and is associated with the element drag. The shed vorticity arises from radial gradients in the bound circulation and lies predominantly in the radial direction. The wake expansion will cause the shed vorticity to have a radial component, so the main difference between wake and shed vorticity is that the latter can follow the local streamlines and so be force-free, whereas the former cannot. A further distinction is that an ideal rotor has shed but no wake vorticity, whereas a real rotor can have a wake without shed vorticity in the unlikely event of constant bound circulation along the blades.

We now test the plausibility of the revised BEMT. The experiments for comparison are those previously mentioned by Krogstad and Adaramola (2012) and Krogstad and Eriksen (2013) on a three-bladed rotor with R = 1.5 m. The airfoil (S826) lift and drag data, extracted from Bartl et al. (2019), cover the range of Reynolds number, Re, of the experiments. C _l (α, Re) and C _d (α, Re) were found by linear interpolation in α and log(Re).

The results for power and thrust using N _BE = 40 are shown in Figure 4. Table 1 lists the variation in C _P and C _T for N _BE = {20, 40, 80} for λ = 6 which is close to λ _opt , and λ = 12, near λ _R . All subsequent results were obtained with N = 40. The experimental results have not been corrected for blockage, which was about 10% (Sarlak et al., 2016). Using computational modeling of the experiment, Sarlak et al. (2016) found that C _P was under-estimated by around 0.05, and the corrected λ _R value was closer to 12 (their Fig. 13). Figure 14 of Sarlak et al. (2016) suggests that a blockage correction would reduce C _T by approximately the same amount. It appears that applying blockage corrections to the experimental data would generally improve the accuracy of the simple model but no correction has been applied. As expected for a BEMT analysis without finite blade corrections, C _P is over-estimated near λ _opt . The largest measured C _P at λ = 6.164 was 0.4481, which is about 10% lower than the value at λ = 6 shown in Table 1. C _T is generally under-estimated, but the trends of power and thrust and the estimate for λ _R ≈ 12 are encouraging. Furthermore, the simple model appears to have no problems in computing through λ _R and with C _T values well in excess of unity.

The two methods for calculating u are shown in Figure 5 at λ = 10, which is the closest to λ _R that measurements were made. Eq. 4 gave higher values of u across the rotor, so these were used to compute C _P and C _T . These should be lower than the values measured one radius downwind by an amount that is difficult to estimate. It is, however, noteworthy that the high values near the edge of the rotor have been reproduced, at least qualitatively. There is considerable discrepancy between the two calculations of u toward the hub, and it is likely that the large negative values are unphysical. They result from helical symmetry, which forces u to have the same sign as w and Γ. Furthermore, it is not clear whether zero or slightly negative u measured downwind of the rotor and nacelle implies similar values at the rotor computed without considering the nacelle.

Negative w indicates that a part of the rotor is operating as a propeller, whose flow normally contracts by a small amount rather than expanding significantly. In the present case, the “propeller region” is near the hub where we would expect v ² to be small. Assuming v ² = 0 and including u ² in the MT calculation forces u = 0. This can be seen by starting with d Q d T * = r ⁡ tan ⁡ ϕ = r 1 − u w + λ r = p , (12) where the first equality comes from dividing the BE part of Eq. 4 by that of Eq. 10 and the second from the blade element velocity triangle in Figure 1. Alternatively, using the right sides of Eqs 4, 10 with v ² = 0 gives d Q d T * = r 1 − u w + λ r − u 2 / w , (13) and so u = 0 for consistency. If u is set to zero in the BEMT calculations whenever w < 0, C _P , C _T , and λ _R increase; at λ = 12, for example, C _P increases to 0.248, C _T to 1.405, and λ _R to be closer to 14. The change is largely due to the reduction in magnitude of the negative w value, which is not shown. Detailed measurements of all velocity components emerging from the rotor would be a valuable aid for improving the modeling. The development of the calculated induced velocities (without setting u = 0 for negative w) with λ is shown in Figure 6. At λ _opt , w is everywhere positive but has a propeller region by λ = 10. The figure suggests the discrepancy between the two methods of determining u becomes more important as λ increases, and there is an increasing part of the rotor where w and Γ are negative. It appears that the runaway is approached more by the overall rotor torque going to zero rather than the torque on all elements being zero. This is consistent with the development of the vortex pitch shown in Figure 7). p decreases with λ, but the very small value at the runaway from Eq. 11 is confined to the tip region.

The comparison with the experimental results in the previous section shows that the wake vorticity term in the revised BEMT provides a plausible explanation for the occurrence of a runaway as a consequence of the blade element drag, and a reasonable description of the λ − dependence of the thrust and power as the runaway is approached. In addition, the high values of the calculated induced axial velocity near the blade tip as the runaway is approached agree qualitatively with the experiment. These high values arise from the effective decoupling of the thrust from the axial velocity through the wake vorticity term which we added to BEMT.

What is not clear is whether the particular form of the wake vorticity term is sufficient, and ignoring the effects of the shed vorticity is justified. The thrust is consistently under-predicted (Figure 4), suggesting that additional shed vorticity terms are needed in the axial momentum equation, as argued by Limacher et al. (2022). We also note that Ebert and Wood (2002) found the tip vortex had sufficient negative angular momentum at runaway to balance the kinetic energy deficit in the remainder of the wake. If this is a general feature, then the angular momentum equation would require revision. Figure 6 shows that the wake considered here had a substantial deficit in kinetic energy at runaway, which must be balanced elsewhere in the wake to prevent it from being extracted by the rotor. It will be necessary to explore the link between this balance and the wake vorticity term before the revision of BEMT can be completed. This would require measurements of the three-dimensional distribution of all three velocity and vorticity components immediately behind the rotor—a task that is far from trivial.

Runaway at high λ may be associated with the blade element lift and drag changing from the two-dimensional values used here. There is a long history of three-dimensional corrections to airfoil data, e.g., Bangga et al. (2023), which should be investigated in future experiments. On the other hand, high λ is likely to reduce the effective blade sweep caused by the large v velocities as the rotor flow expands significantly.

A further issue is a consequence of the large u value; the expansion integral in Eq. 2 requires high levels of the radial velocity, v, near runaway. The calculations of Wood and Limacher (2021) suggest that v is smaller than u near the axis of rotation where v must be zero. v is then comparable to u near the blade tip and greater in the external flow. In the present context, this suggests a complex wake development. On the other hand, it is emphasized that the simple revision to BEMT is based on using a control volume for the axial and angular momentum equations that ends immediately behind the blades. Thus, no relation between the axial velocity leaving the rotor and in the far-wake is required. Because the angular and axial momentum fluxes must be constant downwind of the turbine, there are only two ways the wake can influence the present model: in determining the spanwise variations induced by its vorticity distribution and through the shed vorticity as noted previously. Gupta and Leishman (2005) imply that the wake at high λ (small p) becomes unstable, and this results in the so-called “vortex ring” and “turbulent wake” states and the breakdown of BEMT. The importance of the wake vorticity term suggests otherwise, unless large amounts of radial vorticity are associated with these states. This seems unlikely, however, as the “free wake” model of Gupta and Leishman (2005) does not include the radial vorticity associated with the blade drag. The complexity of the wake at high thrust is also suggested by the simulations of Martínez-Tossas et al. (2022) but again without considering the effect of the wake vorticity term.

Two general aspects of the wake vorticity term and the current analysis are worth mentioning. Eq. 12 gives the vortex pitch as the ratio of the blade element torque to reduced thrust, which also holds for the ideal actuator disc model of Wood and Hammam (2022), for which dT = dT*. Thus, it is not surprising that the runaway is associated with p ↓ 0. Second, the circulation is found from the axial momentum equation, which has significant contributions from the blade element drag. It follows that blade element circulation, in general, is determined partly by drag. This is also a consequence of the angular momentum (Eq. 4) having no wake vorticity term to cancel the drag. A similar result holds for elements of a cascade (Golmirzaee and Wood, 2023).

In conclusion, it is clear that there is still much to learn about wind turbine operation in the high-thrust region. Detailed and high-quality experiments are needed to guide the further development of models for numerical simulation. The present contribution suggests the wake vorticity term in the axial momentum equation will be an important part of improved models.