# Notes on the blade element momentum theory (BEMT) for horizontal axis turbines

M.-A. Dufour

► To cite this version:

M.-A. Dufour. Notes on the blade element momentum theory (BEMT) for horizontal axis turbines. RDT/LHYMAR-BL015, Ifremer, Centre Manche Mer du Nord, 150 quai Gambetta, 62200 Boulogne-sur-Mer. 2025. hal-05385383

HAL Id: hal-05385383

https://hal.science/hal-05385383v1

Submitted on 28 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

RÉPUBLIQUE
FRANÇAISE

Liberté
Egalité
Fraternité

Ifremer

# Notes on the blade element momentum theory (BEMT) for horizontal axis turbines

Marc-Amaury DUFOUR¹

¹Laboratoire d'Hydrodynamique Marine, Boulogne-sur-mer

Révision :
Rev. A 24/11/2025

Diffusion : libre (internet)

Validé par : Grégory GERMAIN

RDT/LHYMAR-BL015

24/11/2025
# Contents

<table><tr><td>1</td><td>Introduction</td><td>3</td></tr><tr><td>2</td><td>Momentum theory (MT): actuator disc model</td><td>3</td></tr><tr><td>3</td><td>Blade element theory (BET)</td><td>14</td></tr><tr><td>4</td><td>Combining MT and BET equations</td><td>19</td></tr><tr><td>5</td><td>Example of application</td><td>21</td></tr><tr><td>6</td><td>Turbines performance assessment</td><td>25</td></tr><tr><td>7</td><td>Conclusion</td><td>26</td></tr><tr><td colspan="2">Bibliography</td><td>26</td></tr></table>

# List of Figures

<table><tr><td>2</td><td>Top view of the actuator disc model.</td><td>4</td></tr><tr><td>3</td><td>Diagram of fluid pressure against x-position.</td><td>4</td></tr><tr><td>4</td><td>Top view of the actuator disc model with the control surface Σ.</td><td>6</td></tr><tr><td>5</td><td>Top view of the actuator disc model with the control surfaces σ and Σ.</td><td>8</td></tr><tr><td>6</td><td>Diagram of fluid velocity against x-position.</td><td>9</td></tr><tr><td>7</td><td>Description of the surface of an annular element.</td><td>10</td></tr><tr><td>8</td><td>Diagram of fluid angular velocity against x-position.</td><td>11</td></tr><tr><td>9</td><td>Top view of the actuator disc model with fluid rotation in the wake.</td><td>11</td></tr><tr><td>10</td><td>Two visualizations of a single rotor blade with geometric notations definition.</td><td>15</td></tr><tr><td>11</td><td>Direct orthonormal bases for correct rotation matrices.</td><td>16</td></tr><tr><td>12</td><td>Blade profile with real-life values for velocity and profile axis angles.</td><td>17</td></tr><tr><td>13</td><td>NACA 0012 profile lift and drag coefficients polar curves.</td><td>22</td></tr></table>

# Acknowledgements

The author acknowledges the support of both IFREMER and Normandy region for co-funding the PhD project (Dufour 2025) during which this BEMT lecture has been given. The author expresses his gratitude towards Dr. Grégory Germain and Prof. Grégory Pinon for the trust they showed throughout the building of this work. The author also acknowledges the support of LOMC laboratory at Université Le Havre Normandie for hosting most part of the work. Lastly, IFREMER is also acknowledged for financing the present fixed term contract enabling to pursue BEMT developments.
# 1 Introduction

The blade element momentum theory, generally abbreviated as BEMT, combines two "competing" theories from the late 19th and beginning of 20th centuries. Those theories were first developed to better understand ship and aircraft propellers behaviour.

On the one hand, the momentum theory (MT) was first introduced by Rankine (1865). This theory was further developed at the beginning of the 20th century by Lanchester (1915), Joukowsky (1920) and Betz (1920). It is a model that relies on a macroscopic viewpoint. The whole streamtube of the flow passing through the turbine (or propeller) is studied. In this model, the rotor is seen as a disc affecting the flow both axially and tangentially. This is the reason why this theory is also known as the actuator disc model.

On the other hand, the blade element theory (BET) was first introduced by Froude (1920) and Drzewiecki (1920). It consists in slicing each blade into elements that can be analysed with a two-dimensional viewpoint. Hence, the analysis reduces to the lift and drag forces evaluation on each blade element. Analysis of the forces is based on tabulated polar curves of lift and drag coefficients and on the relative, also known as true, velocity assessment. This velocity corresponds to local velocity expressed in the blade profile reference frame. Eventually, the two theories (MT and BET) were combined in the late 1930's by Glauert (1935).

State-of-the-art BEMT formalism is presented in the following books by Burton et al. (2011), Sørensen (2016) and Branlard (2017). Recent research articles on the topic encompass the enhancement of the solving process (Ning 2013; Ledoux et al. 2021) and the generalisation of the theory for yawed or highly-loaded rotors (Liew et al. 2024). The reader is invited to have a look at those reference books and papers. The present report explains step by step the demonstrations required to obtain the BEMT equations. This is the transcription of a lecture given at Université Le Havre Normandie, inspired by the lessons given at Durham University (Ingram 2011) and at PennState University (Schmitz 2013). The final set of equations is written to match Ning (2013) formalism. As mentioned in the introductory paragraph, the BEMT can describe the effect of both turbine or propeller on the fluid. However, the present manuscript focuses on turbines with only key aspects specified for propellers modelling.

The momentum theory is studied in Sec. 2. It is sliced into annular elements to match the blade element theory viewpoint detailed afterwards in Sec. 3. A special attention is given to the axis systems definition. The choice of using direct orthonormal bases results in projections that are different from the ones generally described in the literature. This choice simplifies BEMT coding and ensures the compatibility with other frameworks if needed. MT and BET are combined to obtain the numerical formulation of BEMT equations in Sec. 4. The obtained BEMT is applied in Sec. 5 on the example proposed by Ingram (2011) for the reader to get familiar with the present notations. Lastly, some additional definitions on turbine performance assessment are provided in Sec. 6.

# 2 Momentum theory (MT): actuator disc model

In the momentum theory, the turbine action on the flow is modelled by an imaginary porous disc that extracts kinetic energy from the flow. This viewpoint justifies the "actuator disc model" name of this theory. Hence, an axial force ($F_x$), also known as thrust ($T = F_x$), is generated on the disc. Additionally, it is assumed that the porous disc triggers flow
rotation downstream of its position which is related to the turbine torque (Q) extraction.

## 2.1 Axial force or thrust

Fig. 2 represents the streamtube around the porous disc seen from the top. Four x-positions of interest are specified: 1, 2, $x_D$ (disc position), 3 and 4. Positions 2 and 3 are limits that tend to the disc position. This means that $2 \xrightarrow{2<x_D} x_D$ and $3 \xrightarrow{3>x_D} x_D$.

<figure><img src="image_3.png" /></figure>

Figure 2: Top view of the actuator disc model. Positions 2 and 3 are respectively limits to the disc x-position ($x_D$) at the left and at the right.

Upstream conditions are located at x-position 1 and downstream conditions at 4. MT hypothesis are reminded here to indicate some of its inner limitations. First, boundary conditions and porous disc effect on fluid pressure are specified. Fluid pressure variations are then represented as function of x-position in Fig. 3.

* Velocity continuity at the disc position: $V_2 = V_3$;

* Constant pressure at infinity: $p_1 = p_4$;

* Between positions 2 and 3 the turbine (porous disc) is extracting energy from the fluid. It is modelled by a drop in fluid pressure: $p_2 > p_3 = p_2 - \Delta p$ with $\Delta p > 0$.

<figure><img src="image_4.png" /></figure>

Figure 3: Diagram of fluid pressure against x-position.

The second set of hypothesis is needed to ensure that the Bernoulli relation (Eq. (1)) can be applied in the present case.

* Fluid is inviscid;

* Flow is incompressible. It is modelled with a fluid of constant density ($\rho$);
* Flow is permanent, meaning independent of time, and irrotational;

* Gravity effects are negligible compared to other terms;

* The fluid is only subjected to conservative forces. So no viscous forces can be involved. This is valid because the fluid is inviscid. Moreover, the discontinuity in the local pressure force between x-positions 2 and 3 is not conservative as it cannot be derived from a continuous scalar potential. As a consequence, Bernoulli's principle cannot be directly applied between positions 1 and 4.

Under the previous assumptions, the Bernoulli's principle on a streamline reads:

$$ \frac{\partial}{\partial l} \left( \frac{V^2}{2} + \frac{p}{\rho} + gz \right) = 0, \qquad (1) $$

where $l$ represents the 1-dimensional space variable on the streamline, $V$ the local velocity, $p$ the local fluid pressure, $\rho$ the fluid density, and $g$ is the acceleration of gravity. Bernoulli's principle (Eq. (1)) is applied on the streamline $1 \to 2$ shown in Fig. 2.

$$ \int_{1}^{2} \frac{\partial}{\partial l} \left( \frac{V^2}{2} + \frac{p}{\rho} + gz \right) dl = \int_{1}^{2} 0 dl \qquad (2) $$

$$ \text{Eq. (2)} \Rightarrow \frac{V_2^2}{2} + \frac{p_2}{\rho} - \left( \frac{V_1^2}{2} + \frac{p_1}{\rho} \right) = 0 \Leftrightarrow \frac{V_2^2 - V_1^2}{2} + \frac{p_2 - p_1}{\rho} = 0 \qquad (3) $$

Injecting the hypothesis $V_2 = V_3$ and $p_1 = p_4$ into Eq. (3) leads to Eq. (4).

$$ \frac{V_3^2 - V_1^2}{2} + \frac{p_2 - p_4}{\rho} = 0 \qquad (4) $$

In the same way, Bernoulli's principle (Eq. (1)) is applied on the streamline $3 \to 4$.

$$ \frac{V_4^2 - V_3^2}{2} + \frac{p_4 - p_3}{\rho} = 0 \qquad (5) $$

Eventually, Eq. (4) and Eq. (5) are added to obtain a pressure-velocity relation (Eq. (6)).

$$ \text{Eq. (4) + Eq. (5)} \rightarrow \frac{V_4^2 - V_3^2 + V_3^2 - V_1^2}{2} + \frac{p_4 - p_3 + p_2 - p_1}{\rho} = 0 \qquad (6) \\ \Rightarrow \Delta p = p_2 - p_3 = \frac{1}{2} \rho (V_1^2 - V_4^2) $$

By definition, the pressure force on a given surface equates to the differential of the pressure times the surface itself. Hence, the pressure force applied on the porous disc by the fluid is $d\vec{F}_x = \Delta p dS \vec{e}_x$ with $dS$ a small surface element of the porous disc. This pressure force is carried by $\vec{e}_x$ unit vector as the pressure varies only along $x$ axis. So, by injecting Eq. (6) into the pressure force definition, the axial force becomes:

$$ dF_x = (p_2 - p_3) dS = \frac{1}{2} \rho (V_1^2 - V_4^2) dS. \qquad (7) $$
This relation for the thrust is the basis of the awaited result. However, the reader must note that the velocity $V_4$, present in Eq. (7), is not known. Then, finding an expression of $V_4$ that depends on $V_1$ and possibly on some other parameter is the objective of the following demonstrations.

## Definition of axial induction coefficient (a)

First, the axial induction coefficient (a) is defined. It is a ratio that characterizes the velocity variation between $V_2$ and $V_1$. It is assumed that $V_1 \neq 0$.

$$a = \frac{V_1 - V_2}{V_1} \qquad (8)$$

Then, using Eq. (8), $V_2$ is expressed depending on the axial induction coefficient (a).

$$aV_1 = V_1 - V_2 \Leftrightarrow V_2 = V_1 - aV_1 \Leftrightarrow V_2 = V_1(1 - a) \qquad (9)$$

## Relation between upstream ($V_1$) and downstream velocity ($V_4$)

The following developments establish the relation: $V_4 = V_1 (1 - 2a)$, providing a way to evaluate Eq. (7). It may be the “hardest” part of demonstrating the BEMT equations.

### Mass flow conservation

Fig. 4 represents the decomposition of $\Sigma$ surface. This closed surface is composed of three specific surfaces: $\Sigma = S_1 \cup S_{ext} \cup S_4$. The surface $S_{ext}$ is defined as a streamtube. A streamtube is a surface made of lines that all are tangential to the flow at any point. In other words, the vectors normal to this surface ($\vec{n}$) are orthogonal to the velocity at any considered location on $S_{ext}$. Consequently, on $S_{ext}$ surface, $\vec{V} \cdot \vec{n} = 0$. The upstream section is $S_1$ and the downstream one is $S_4$. For these two surfaces, their normal vectors ($\vec{n}$) are respectively opposed or aligned to the velocity vector ($\vec{V}$). Note that, by definition, the normal vector ($\vec{n}$) points towards the outer of the control surface. The integral mass flow rate ($D_m$) is defined by:

$$D_m = \iint_{\Sigma} \rho \vec{V} \cdot \vec{n} dS. \qquad (10)$$

<figure><img src="image_5.png" /><figcaption>Figure 4: Top view of the actuator disc model with the control surface $\Sigma$. This control surface is defined by: $\Sigma = S_1 \cup S_{ext} \cup S_4$.</figcaption></figure>
The integral mass flow rate ($D_m$) is decomposed following the definition of the special surface $\Sigma$. Note also that, at a given position $x_0$, velocity is constant in the plane $(x_0, y, z)$ because of the axial symmetry of the problem and of the uniform velocity distribution assumed at this point. Hence, the mass flow rate can be evaluated as:

$$D_m = \iint_{S_1 \cup S_{ext} \cup S_4} \rho \vec{V} \cdot \vec{n} dS = \iint_{S_1} \rho \vec{V}_1 \cdot \vec{n}_1 dS + \iint_{S_{ext}} \rho \underbrace{\vec{V} \cdot \vec{n}}_{=0} dS + \iint_{S_4} \rho \vec{V}_4 \cdot \vec{n}_4 dS \\ = - \iint_{S_1} \rho V_1 dS + \iint_{S_4} \rho V_4 dS \\ = \rho(-S_1 V_1 + S_4 V_4). \qquad (11)$$

Eq. (10) represents the **algebraic** integral mass flow rate. This means that the mass that goes **out** is counted positively while the mass that goes **in** is counted negatively. Hence, for a **closed surface of integration**, like the $\Sigma$ surface (Fig. 4), the value of the integral is always zero. To understand that the algebraic mass flow rate evaluated along a closed surface equates to zero, a valuable analogy is the flow in a pipe. In this case, the closed control surface is defined by the pipe boundaries ($S_{ext}$) and the input and output sections. No choice: everything that goes in, goes out. So:

$$D_m = 0. \qquad (12)$$

Equating Eq. (11) and Eq. (12) provides immediately: $\rho(S_1 V_1 - S_4 V_4) = 0$. Then, as the density $\rho$ is a constant different from zero, it reads:

$$S_1 V_1 = S_4 V_4. \qquad (13)$$

It is the conservation of the cross-section mass, or volume, flow rate. As the density ($\rho$) is constant, the flow conservation between the input ($S_1$) and output ($S_4$) surfaces can be expressed both in terms of mass or volume. For the reader more familiar with Navier-Stokes equations, Eq. (13) is a consequence of the continuity equation under the assumptions provided at the beginning. Let $\dot{m}$ be the mass flow rate through any cross-section of the streamtube. It corresponds to the mass that goes through the cross-section during a given amount of time. It can be defined through $S_1$ surface at $x$-position 1:

$$\dot{m} = \rho S_1 V_1. \qquad (14)$$

The mass flow conservation is valid for any cross-section of the streamtube. Remembering that the subscript $D$ denotes the quantities of interest at the porous disc $x$-position, the disc surface is denoted $S_D = S_2 = S_3$. This surface is normal to the flow and corresponds to a cross-section of the streamtube. The velocity magnitude at the disc $x$-position is $V_D = V_2 = V_3$. Then, according to Eq. (13) and Eq. (14), it reads:

$$\dot{m} = \rho S_1 V_1 = \rho S_D V_D = \rho S_4 V_4. \qquad (15)$$
## Momentum conservation of a closed system

Let $\vec{P}$ be the momentum quantity of a closed system. It equates to the contained mass times its velocity. To derive a momentum balance, the point of view is changed from the one previously developed. It is now assumed that the surface $\Sigma$ is not fixed. It travels with the flow. It contains at each time $(t)$ the exact same set of fluid particles. In other words, the fluid particles are tracked in time by the moving surface. This moving surface can be seen as a nest. Hence, the particles contained in the surface define a closed system. As a consequence, all the particles that are contained within $\Sigma$ at time $t$, were contained within $\sigma$ at time $t - dt$. Fig. 5 details this new viewpoint. The view in Fig. 5 shows $\sigma$ surface associated to the closed system of fluid particles at time $t - dt$. After a small time increment $(dt)$, all those particles are advected and they are now contained in the $\Sigma$ surface at time $t$. Hence, the momentum quantity associated to each time becomes:

$$ \text{Over } \sigma: \vec{P}(t - dt) = \vec{P}_{\sigma} = \vec{P}_{1} + \vec{P}_{m}, \qquad (16) $$

$$ \text{Over } \Sigma: \quad \vec{P}(t) = \vec{P}_{\Sigma} = \vec{P}_{m} + \vec{P}_{4}. \qquad (17) $$

<figure><img src="image_6.png" /><figcaption>Figure 5: Definition of a closed system as a given number of fluid particles that are tracked. They are contained within $\Sigma$ at time $t$, and were contained within $\sigma$ at time $t - dt$.</figcaption></figure>

To obtain the momentum balance, the principle is to differentiate the momentum of the closed system between the times $t - dt$ and $t$. By convention, when differentiating $\vec{P}(t) - \vec{P}(t - dt)$, a positive result is obtained if the fluid gained momentum, negative otherwise. This is linked to the definition of $\vec{P}$ which quantifies the momentum contained within the fluid closed system. Then, by differentiating with the opposite sign, $\vec{P}(t - dt) - \vec{P}(t)$, the result is positive if the porous disc gained/absorbed momentum. The time-step $dt$ is small enough so that the all fluid particles contained in $\sigma \setminus \Sigma$, corresponding to $\vec{P}_1$, travel at velocity $\vec{V}_1$ and that all the particles contained in $\Sigma \setminus \sigma$, corresponding to $\vec{P}_4$, travel at velocity $\vec{V}_4$. Then, using Eq. (16) and Eq. (17), it reads:

$$ \begin{aligned} d\vec{P} = \vec{P}(t - dt) - \vec{P}(t) &= \vec{P}_m + \vec{P}_1 - (\vec{P}_m + \vec{P}_4) = \vec{P}_1 - \vec{P}_4 \\ &= m_1 \vec{V}_1 - m_4 \vec{V}_4 \qquad (18) \\ &= \rho S_1 V_1 dt \vec{V}_1 - \rho S_4 V_4 dt \vec{V}_4. \end{aligned} $$

Injecting the result of the mass flow conservation (Eq. (15)) into Eq. (18), yields:
$$d\vec{P} = \dot{m}dt\vec{V}_1 - \dot{m}dt\vec{V}_4 = \dot{m}(\vec{V}_1 - \vec{V}_4)dt. \qquad (19)$$

According to Newton's laws of motion, the derivative of the momentum with respect to time equates to the force exerted on the closed system. The only force to which the system is submitted is the porous disc action, already defined as $\vec{F}_x$. Hence:

$$\vec{F}_x = \frac{d\vec{P}}{dt}. \qquad (20)$$

So, as all vectors are carried by $\vec{e}_x$ unit vector, Eq. (21) is obtained by combining Newton's law (Eq. (20)) with the result from the momentum differential (Eq. (19)).

$$F_x = \dot{m}(V_1 - V_4) \qquad (21)$$

The result from Bernoulli's principle (Eq. (7)), integrated over the disc surface ($S_D$), and the result from the momentum conservation (Eq. (21)), with the mass flow rate ($\dot{m}$) being expressed at the disc $x$-position, are equated. It leads to Eq. (22).

$$\begin{align}
\frac{1}{2}\rho(V_1^2 - V_4^2)S_D = \dot{m}(V_1 - V_4) &\Leftrightarrow \frac{1}{2}\rho S_D \vec{V}_4 (V_1 + V_4) = \rho S_D \vec{V}_2 (V_1 - V_4) \tag{22} \\
&\Leftrightarrow \frac{1}{2}(V_1 + V_4) = V_2 \nonumber
\end{align}$$

It is interesting to note that the velocity at the disc $x$-position ($V_D = V_2 = V_3$) is equal to the average between the upstream velocity ($V_1$) and the downstream one ($V_4$). Following this finding, a diagram of the velocity as function of the $x$-position is presented in Fig. 6.

<figure><img src="image_7.png" /><figcaption>Figure 6: Diagram of fluid velocity against $x$-position.</figcaption></figure>

The axial induction coefficient definition (Eq. (9)) is used to transform $V_2$ as a function of $V_1$. Consequently, the desired expression for $V_4$ as function of $V_1$ and the axial induction coefficient ($a$) is obtained in Eq. (23).

$$\frac{1}{2}(V_1 + V_4) = V_1(1 - a) \Leftrightarrow \frac{1}{2}V_4 = V_1\left(\frac{1}{2} - a\right) \Leftrightarrow V_4 = V_1(1 - 2a) \qquad (23)$$
## Final thrust force expression

The result on $V_4$ velocity (Eq. (23)) is used within Eq. (7) to eventually obtain an expression for $dF_x$ that only depends upon the surface considered (dS), the upstream velocity ($V_1$) and the axial induction coefficient (a):

$$ 
\begin{aligned}
dF_x &= \frac{1}{2} \rho (V_1^2 - V_1^2 (1 - 2a)^2) dS \\
&= \frac{1}{2} \rho (V_1^2 - V_1^2 (1 - 4a + 4a^2)) dS \\
&= \frac{1}{2} \rho V_1^2 4a (1 - a) dS.
\end{aligned} 
\quad (24) $$

The final step is to express this force on a surface (dS) that is adapted to match the blade element theory (BET) formalism that will be seen and used afterwards. A relevant choice is to look at an annular element, or ring. By definition, a segment of a ring of width (dr) and angular spacing (dθ), as depicted in Fig. 7a, has a small surface of: $rdrd\theta$. Once integrated on the $2\pi$ range of the θ angles, the annular element provides the surface of the cross-section of an annular streamtube: $dS_{el} = r2\pi dr$ as shown in Fig. 7b.

<figure><img src="image_8.png" /><figcaption>Figure 7: Description of the surface of an annular element. The small surface of a segment of the annular element (a) and the complete surface of the ring (b).</figcaption></figure>

As a consequence, the final equation for the thrust force on the cross-section of an annular streamtube ($dF_x$) is given by:

$$ dF_x = \frac{1}{2} \rho V_1^2 4a (1 - a) dS_{el} = \rho \pi 4a (1 - a) V_1^2 r dr. \quad (25) $$

## 2.2 Torque evaluation using the rotating annular streamtube

The produced torque (Q) is evaluated via the actuator disc action in the tangential direction. The studied system is an annular element located of radius $r$ at x-position 3. It is assumed that the passing blades generate an opposite rotation of the fluid downstream of the disc x-position ($x_D$). It is represented by the actuator disc swirling action. Specific regions are identified regarding the fluid angular velocity. Those fluid angular velocity variations against x-positions are evidenced in Fig. 8.

* From 1 → 2 there is no fluid rotation;

* From 3 → 4 the fluid rotates at an angular velocity $\omega > 0$, of rotation vector $-\omega e_x$;
* At the specific disc x-position ($x_D$), the fluid rotates at $\omega/2$, an “average” angular velocity between upstream and downstream values, see Fig. 8.

<figure><img src="image_9.png" /><figcaption>Figure 8: Diagram of fluid angular velocity against x-position.</figcaption></figure>

Following Fig. 8, the studied annular element at the x-position 3 rotates with a tangential velocity of $-r\omega\vec{e}_\theta$. This phenomenon is illustrated with the help of Fig. 9. The streamline representation highlights the fluid angular velocity in the wake. It may disturb the reader compared to what is previously presented in Fig. 2. As a matter of fact, Fig. 2 focuses on the axial force using Bernoulli's principle. The formulation of this principle on a streamline is independent from the path followed. It only relies on the starting and end points. Hence, Fig. 2 represents a first approximation of the fluid motion that enables to correctly evaluate Bernoulli's principle.

<figure><img src="image_10.png" /><figcaption>Figure 9: A streamline (red) is represented in this top view of the actuator disc model. It highlights the fluid angular velocity due to the disc influence.</figcaption></figure>

## Tangential induction coefficient and hypothesis

As previously mentioned, at $x_D$, the fluid tangential velocity is: $-r\omega/2\vec{e}_\theta$. In order to quantify the amount of angular velocity gained by the fluid due to the disc swirling effect on the flow, the tangential induction coefficient ($a'$) is defined as follow:

$$a' = \frac{\omega}{2\Omega}. \qquad (26)$$

The hypothesis related to the present analysis are summarized:
* Fluid is inviscid;

* Flow is incompressible. It is modelled with a fluid of constant density ($\rho$);

* Flow is permanent, meaning independent of time, and irrotational;

* Angular velocity of the blades is denoted $\Omega > 0$. Rotation vector is: $\Omega \vec{e}_x$;

* Fluid angular velocity in the wake is denoted $\omega$ with the rotation vector: $-\omega \vec{e}_x$;

* The studied system is a fluid annular element of width $dr$ and zero thickness, located in the wake at $x$-position 3. The inner radius of this annular element is $R_1$ and the outer radius is $R_2$. Its average radius is $r$ as shown in Fig. 7;

* This zero thickness annular element has a homogenous (constant between $R_1$ and $R_2$) surface mass ($\zeta(t)$) expressed in $kg \cdot m^{-2}$, which depends upon the time ($t$). Hence, the mass of the element is: $m_{el} = S_{el}\zeta(t)$. Note that the differential notation is not applied as, at first, the annular element is defined by $R_1$ and $R_2$ boundaries;

* Use of polar coordinates: $y = r \cos(\theta)$ and $z = r \sin(\theta)$.

## Torque definition

Newton's law of motion applied to angular momentum provides the torque ($Q$) definition:

$$ \vec{Q} = \frac{d\vec{L}_{el}}{dt}, \qquad (27) $$

with $\vec{L}_{el}$ being the net angular momentum of the wake annular element. It is reminded that this wake annular element is located at $x$-position 3, meaning its rotation vector is: $-\omega \vec{e}_x$. The net angular momentum of such an annular element becomes:

$$ \vec{L}_{el} = -\omega I_{el} \vec{e}_x. \qquad (28) $$

So, the next step is to evaluate the moment of inertia of the annular element ($I_{el}$).

## Moment of inertia of an annular element

The objective is to evaluate and simplify the moment of inertia ($I_{el}$) of a thin annular element of width ($dr$). The question that arises is: can this moment of inertia be simplified as the moment of inertia of a line circle ($m_{el}r^2$)?

The small mass ($dm$) linked to the small surface of an annular element segment, as represented in Fig. 7a, is:

$$ dm = \zeta(t)rdrd\theta. \qquad (29) $$

Eq. (30) is the definition of the moment of inertia around the rotation axis $\vec{e}_x$. The hypothesis on geometry and mass are immediately applied.

$$ I_{el} = \iint (y^2 + z^2) dm = \zeta \int_{R_1}^{R_2} \int_{0}^{2\pi} (r^2 \cos^2(\theta) + r^2 \sin^2(\theta)) r dr d\theta \qquad (30) $$

Then, the variables are grouped under each related integral symbol, leading to Eq. (31).
$$ \frac{I_{el}}{\zeta} = \int_{R_1}^{R_2} r^3 dr \int_{0}^{2\pi} \underbrace{(\cos^2(\theta) + \sin^2(\theta))}_{=1} d\theta = 2\pi \int_{R_1}^{R_2} r^3 dr \quad (31) $$

To find a primitive, a trick is to rely on derivation knowledge. In particular, $(\frac{r^4}{4})' = r^3$, so :

$$ \int_{R_1}^{R_2} r^3 dr = \left[ \frac{r^4}{4} \right]_{R_1}^{R_2} = \frac{R_2^4}{4} - \frac{R_1^4}{4} = \frac{R_2^4 - R_1^4}{4} = \frac{(R_2^2 - R_1^2)(R_2^2 + R_1^2)}{4}. \quad (32) $$

The results of Eq. (31) and Eq. (32) are gathered to obtain the moment of inertia of the annular element:

$$ I_{el} = 2\pi\zeta \frac{(R_2^2 - R_1^2)(R_2^2 + R_1^2)}{4} = \frac{1}{2}\zeta\pi (R_2^2 - R_1^2)(R_2^2 + R_1^2). \quad (33) $$

The objective is to evidence the mass of the annular element in Eq. (33). Recalling the hypothesis, the mass of the annular element is: $m_{el} = S_{el}\zeta(t)$, with $\zeta$ mass density per unit surface. So, the only thing to do is evaluating the annular element surface ($S_{el}$). This surface can be interpreted as the difference between the surfaces of the largest and of the smallest discs.

$$ S_{el} = \pi R_2^2 - \pi R_1^2 = \pi (R_2^2 - R_1^2) \quad (34) $$

The surface of the annular element (Eq. (34)) is injected into Eq. (33). The moment of inertia of the annular element is eventually obtained in Eq. (35).

$$ I_{el} = \frac{1}{2} m_{el} (R_2^2 + R_1^2) \quad (35) $$

The moment of inertia of the annular element is known. However, the considered annular element is of small width ($dr$). Hence, the inner and outer radii can be written as follow:

$$ R_1 = r - \frac{dr}{2}, \quad (36) $$

$$ R_2 = r + \frac{dr}{2}. \quad (37) $$

So the annular element moment of inertia becomes:

$$ I_{el} = \frac{1}{2} m_{el} \left( r^2 + r dr + \frac{dr^2}{4} + r^2 - r dr + \frac{dr^2}{4} \right) = \frac{1}{2} m_{el} \left( 2r^2 + \frac{dr^2}{2} \right). \quad (38) $$

Only first order contributions are considered, neglecting $dr^2$ term. Thus, Eq. (39) is obtained, proving that, at the first order of approximation, the moment of inertia of an annular element of small width ($dr$) is equal to the moment of inertia of a circle.

$$ I_{el} = m_{el} r^2 \quad (39) $$
Now that the moment of inertia of the annular element is established, the next step is to evaluate the torque contained in such a rotating annular fluid element.

## Torque of an annular fluid element

By injecting the moment of inertia of the annular element (Eq. (39)) and the net angular momentum definition (Eq. (28)) into Newton's law of motion (Eq. (27)), the torque of the annular element becomes:

$$ \vec{Q} = -\omega \frac{dI_{el}}{dt} \vec{e}_x = -\omega r^2 \frac{dm_{el}}{dt} \vec{e}_x. \qquad (40) $$

Is is reminded that the mass is the only time-dependent quantity as the definition of the mass density per unit surface highlights it. Then, $dm_{el}/dt$ corresponds to the mass flow passing through the annular element. This term is easily evaluated using the mass flow definition of Sec. 2.1. The velocity continuity, stating that: $V_2 = V_3$, is reminded.

$$ \frac{dm_{el}}{dt} = S_{el}V_3\rho = \rho V_2 S_{el} = \rho V_1 (1 - a) S_{el} \qquad (41) $$

Thus, the mass flow through this annular element (Eq. (41)) is injected into torque expression (Eq. (40)) to obtain Eq. (42). It is important to recall that the annular element is of small width ($dr$). Hence, the annular element surface is a small quantity depending on the width: $dS_{el} = 2\pi rdr$ (Fig. 7b). The differential notation is used.

$$ d\vec{Q} = -\omega r^2 \rho V_1 (1 - a) dS_{el} \vec{e}_x = -2\pi\omega\rho V_1 (1 - a) r^3 dr \vec{e}_x \qquad (42) $$

In the end, to obtain the torque of a disc annular element (Eq. (43)), the remaining unknowns should be the induction coefficients, be they axial ($a$) or tangential ($a'$). So, the definition of the tangential induction coefficient (Eq. (26)) is injected into Eq. (42). Moreover, as the studied system has been the fluid annular element, the torque associated to its solid counterpart (the disc annular element) is of opposite sign. It reads:

$$ d\vec{Q} = 4\pi\Omega\rho V_1 a' (1 - a) r^3 dr \vec{e}_x. \qquad (43) $$

# 3 Blade element theory (BET)

Contrary to the momentum theory, the blade element theory (BET) must appear clearer to people who are used to lift and drag forces on lifting profiles. The blade element theory consists in establishing expressions for both the thrust ($dF_x = dT$) and the torque ($dQ$) for all the blade elements contributing at a given local radius $r$ as shown for one blade in Fig. 10a. This viewpoint consists in “slicing” each blade of the rotor in a given number of sections ($N_s$). As blade geometry, including the hub radius ($R_h$) and blade length ($L_b$), are properties of the studied turbine, the elements width ($dr$) is defined by: $dr = L_b/N_s$. The angle related to turbine rotation is $\theta$. The chord ($c$) of a profile is defined as the distance between its leading and trailing edges, see Fig. 10b. The twist plus pitch angle ($\gamma$) is represented for the tip blade profile in Fig. 10b. Fig. 10 gathers all those properties represented for a single blade of a rotor. Axial velocity at the rotor $x$-position ($V_2 = V_1(1 - a)$) is represented as well as tangential velocity at rotor tip due to rotation. Both velocities include the induction factors that are function of the local blade section.
<figure><img src="image_11.png" /><figcaption>(a) Back view of a blade with hub radius ($R_h$), blade length ($L_b$) and $r_b \in [0, L_b]$ and local radius $r = R_h + r_b$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Side view of a blade with the twist plus pitch angle ($\gamma$) and chord ($c$) defined at the tip blade profile.</figcaption></figure>

Figure 10: Back (a) and side (b) views of a single rotor blade. Note that the axial (a) and tangential (a') inductions are dependent upon the radius (r) at which they are evaluated. It means that, for instance, axial velocity ($V_1(1-a)$) is different at each blade section.

## 3.1 Axis systems definition

Fig. 11 depicts the axis systems with the viewpoint of reference: projections in the first quarter of direct orthonormal bases. It is of the utmost importance because it enables to define all the projections in a mathematically consistent manner. As shown in Fig. 10a, $\theta$ angle represents the turbine rotation along $\vec{e}_x$ axis. The associated representation in Fig. 11a leads to the definition of the rotation matrix in Eq. (44). In addition, some rotations related to a given blade section are defined in Fig. 11b. Those rotations along the radial axis ($\vec{e}_r$) are described by the rotation matrix of second kind defined in Eq. (45). The angle associated to the profile twist plus pitch angle is $\gamma$. The angle describing the relative/true velocity vector is $\phi$. While Fig. 11 schematics are very helpful to accurately establish the rotation matrices (Eq. (44) and Eq. (45)), they do not represent the physical case.

$$ \begin{pmatrix} \vec{e}_x \\ \vec{e}_r \\ \vec{e}_\theta \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos(\theta) & \sin(\theta) \\ 0 & -\sin(\theta) & \cos(\theta) \end{pmatrix} \begin{pmatrix} \vec{e}_x \\ \vec{e}_y \\ \vec{e}_z \end{pmatrix} = M_1(\theta) \begin{pmatrix} \vec{e}_x \\ \vec{e}_y \\ \vec{e}_z \end{pmatrix} \quad (44) $$

$$ \begin{pmatrix} \vec{e}_{x\phi} \\ \vec{e}_r \\ \vec{e}_\phi \end{pmatrix} = \begin{pmatrix} \cos(\phi) & 0 & -\sin(\phi) \\ 0 & 1 & 0 \\ \sin(\phi) & 0 & \cos(\phi) \end{pmatrix} \begin{pmatrix} \vec{e}_x \\ \vec{e}_r \\ \vec{e}_\theta \end{pmatrix} = M_2(\phi) \begin{pmatrix} \vec{e}_x \\ \vec{e}_r \\ \vec{e}_\theta \end{pmatrix} \quad (45) $$
<figure><img src="image_13.png" /></figure>

(a) Bases definition for rotations along $\vec{e}_x$ axis.

<figure><img src="image_14.png" /></figure>

(b) Bases definition for rotations along $\vec{e}_r$ axis.

<caption>Figure 11: Direct orthonormal bases associated to the fixed Cartesian reference frame ($O, \vec{e}_x, \vec{e}_y, \vec{e}_z$), the blades rotating frame ($O, \vec{e}_x, \vec{e}_r, \vec{e}_\theta$), the blade profile axis ($O, \vec{e}_{x\gamma}, \vec{e}_r, \vec{e}_\gamma$) and the relative/true velocity ($O, \vec{e}_{x\phi}, \vec{e}_r, \vec{e}_\phi$).</caption>

As previously mentioned, Fig. 11b does not represent the real-life vectors. Hence, the priority is to represent the physical case and express those vectors using the presented coordinate systems. Due to the turbine rotation, assumed in the mathematically positive (counterclockwise) direction, the velocity seen by the blade profile, called relative or true velocity, contains a component along $-\vec{e}_\theta$. The $\vec{e}_x$ velocity component is known as $V_2$ following Sec. 2. From these observations, the approximate value of the true velocity angle ($\phi$) is easily guessed: $\phi \in [+\pi/2, +\pi]$. It is evidenced by Fig. 12 with the velocity vectors highlighted in red. Moreover, as shown in Fig. 12, each velocity component is denoted with a subscript related to its direction: respectively $V_x$ and $V_\theta$. As both components are known as functions of the induction coefficients following Sec. 2 developments, the relative/true velocity vector ($\vec{V}_t$) can be defined as:

$$ \vec{V}_t = V_t \vec{e}_\phi = V_x \vec{e}_x + V_\theta \vec{e}_\theta = V_1(1-a)\vec{e}_x - r\Omega(1+a')\vec{e}_\theta, \quad (46) $$

with r the radius of the studied blade section. Furthermore, as the computation will require the $\phi$ angle evaluation, it is obtained by means of the velocity components. Using Fig. 11b and Eq. (46), the following trigonometric relations are easily obtained:

$$ \begin{aligned} \sin(\phi) &= \frac{V_x}{V_t} & \cos(\phi) &= \frac{V_\theta}{V_t} & \tan(\phi) &= \frac{V_x}{V_\theta} \tag{47a-c} \\ &= \frac{V_1(1-a)}{V_t}, & &= -\frac{r\Omega(1+a')}{V_t}, & &= -\frac{V_1(1-a)}{r\Omega(1+a')}. \end{aligned} $$

When using one of the above relations (Eq. (47a) to Eq. (47c)) to find $\phi$ angle, a special care is taken on the definition domains of the inverse trigonometric functions. Those domains are reminded in Eq. (48a) to Eq. (48c). The choice of either $\arcsin$ (Eq. (48a)), $\arccos$ (Eq. (48b)) or $\arctan$ (Eq. (48c)) functions to obtain $\phi$ angle depends on the respective signs of $V_x$ and $V_\theta$ components. In the case of a turbine, $\arccos$ function (Eq. (48b))
should be preferred as being defined whatever $V_\theta$ sign is, noting that $V_t > 0$ by construction. It then means that $\phi$ angle value is independent of the rotation direction (clockwise or counterclockwise). In the case of a propeller, assuming $V_x$ is positive, arccos function (Eq. (48b)) is also the best option. However, if $V_x$ is negative, arctan function (Eq. (48a)) should be used with a special care on the sign of $V_\theta$ component.

$$ \arcsin : [-1, 1] \rightarrow [-\frac{\pi}{2}, \frac{\pi}{2}] \quad \arccos : [-1, 1] \rightarrow [0, \pi] \quad \arctan : \mathbb{R} \rightarrow [-\frac{\pi}{2}, \frac{\pi}{2}] \quad (48a-c) $$

Based on Fig. 11b or Fig. 12, the angle of attack ($\alpha$) is defined as the angle between profile axis and true velocity, leading to Eq. (49). It is different between turbine and propeller cases due to the orientation of the profile suction, or upper, surface. The suction surface is oriented towards $\vec{e}_{x\gamma}$ in the turbine case, as shown in Fig. 11b and Fig. 12. However, the profile suction surface points towards $-\vec{e}_{x\gamma}$ for a propeller with $V_x > 0$.

$$ \alpha(\phi) = \begin{cases} \gamma - (-\pi + \phi) = \pi + \gamma - \phi & \text{if turbine,} \\ 2\pi - (\pi + \gamma - \phi) = \pi + \phi - \gamma & \text{if propeller.} \end{cases} \quad (49) $$

Finally, to complete the blade profile description in real-life conditions, the local twist plus pitch angle is adapted. It is positioned to obtain small but positive angle of attack values, ensuring that profiles work in their region of highest performance. For both turbine or propeller rotating in the mathematically positive direction (counterclockwise, Fig. 11a), it is inferred that the twist plus pitch angle is: $\gamma \in [-\pi/2, 0]$. Hence, Fig. 12 represents the real-life angles for both velocity and profile of a turbine case, see the blade profile upper surface orientation. Compared to Fig. 11b, Fig. 12 thus provides the most common angle of attack representation. It directly justifies Eq. (49) and simplifies its understanding.

<figure><img src="image_15.png" /><figcaption>Figure 12: Blade profile representation with real-life values for velocity and profile twist plus pitch. The true velocity ($\vec{V}_t$) components, $V_x > 0$ and $V_\theta < 0$, are highlighted in red. The local twist plus pitch angle ($\gamma$) is negative.</figcaption></figure>

## 3.2 Aerodynamic loads analysis to obtain thrust and torque

With the axis systems being defined, the loads on the blade profile can now be assessed. By construction, the drag vector is carried by the unit vector of the true velocity ($\vec{e}_\phi$). In
the case of turbines, the unit vector that carries the lift is $\vec{e}_L = -\vec{e}_{x\phi}$. It is chosen to point towards the profile suction surface for small angle of attack ($\alpha$) values, as represented by Fig. 12. However, if a propeller is considered with its upper surface positioned towards the other side, then $\vec{e}_L = +\vec{e}_{x\phi}$. Following the classic lift and drag definition (Katz et al. 2001; Giovannini et al. 2025) and focusing on the turbine case, the forces on the studied blade element are defined by:

$$d\vec{L} = dL\vec{e}_L = -dL\vec{e}_{x\phi} = -\frac{1}{2}\rho cdrV_t^2 C_L(\alpha(\phi))\vec{e}_{x\phi}, \quad (50)$$

$$d\vec{D} = dD\vec{e}_{\phi} = \frac{1}{2}\rho cdrV_t^2 C_D(\alpha(\phi))\vec{e}_{\phi}, \quad (51)$$

with $c$ being the chord of the local blade element and $c \times dr$ the characteristic surface of the studied blade element. Then, the goal is to evaluate the associated axial and tangential forces. Those forces are defined by Eq. (52) resorting to the inverse of $M_2$ rotation matrix (Eq. (45)), the lift (Eq. (50)) and drag (Eq. (51)) definitions. As those forces are element forces of **one** blade element among the other blades contributions, a subscript 1 is added to the notation. Again, the lift direction is specific to a turbine case.

$$ \begin{pmatrix} dF_{x,1} \\ 0 \\ dF_{\theta,1} \end{pmatrix} = M_2^{-1}(\phi) \begin{pmatrix} -dL \\ 0 \\ dD \end{pmatrix} = M_2(-\phi) \begin{pmatrix} -dL \\ 0 \\ dD \end{pmatrix} \\ = \frac{1}{2} \rho cdr V_t^2 \begin{pmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \end{pmatrix} \begin{pmatrix} -C_L(\alpha(\phi)) \\ 0 \\ C_D(\alpha(\phi)) \end{pmatrix} \quad (52) \\ = \frac{1}{2} \rho cdr V_t^2 \begin{pmatrix} -\cos(\phi)C_L(\alpha(\phi)) + \sin(\phi)C_D(\alpha(\phi)) \\ 0 \\ \sin(\phi)C_L(\alpha(\phi)) + \cos(\phi)C_D(\alpha(\phi)) \end{pmatrix} $$

The definition that links a force with the associated torque yields:

$$d\vec{Q}_1 = r\vec{e}_r \wedge d\vec{F}_{\theta,1} = r dF_{\theta,1}(\underbrace{\vec{e}_r \wedge \vec{e}_{\theta}}_{=\vec{e}_x}), \quad (53)$$

with the cross/vector product denoted $\wedge$. Let $B$ be the number of blades. The axial force (or thrust, $d\vec{F}_x$), and torque ($d\vec{Q}$) summed for the contributions of each blade at a given radius ($r$) are eventually obtained. It reads:

$$d\vec{F}_x = \sum_{k=1}^{B} d\vec{F}_{x,1} = B d\vec{F}_{x,1} = \frac{B}{2} \rho cdr V_t^2 [-\cos(\phi)C_L(\alpha(\phi)) + \sin(\phi)C_D(\alpha(\phi))] \vec{e}_x, \quad (54)$$

$$d\vec{Q} = \sum_{k=1}^{B} d\vec{Q}_1 = B d\vec{Q}_1 = r \frac{B}{2} \rho cdr V_t^2 [\sin(\phi)C_L(\alpha(\phi)) + \cos(\phi)C_D(\alpha(\phi))] \vec{e}_x. \quad (55)$$
# 4 Combining MT and BET equations

This part is based upon Ning (2013) article to match his formalism. It adds the possibility to rely on a root-finding algorithm for numerical resolution.

## 4.1 Equating axial force expressions

As stated in the subsection title, the objective is to equate Eq. (25) (MT) and Eq. (54) (BET) that both describe the axial force, or thrust, for a disc annular element at a given radius ($r$). The idea is to find a relation for the axial induction coefficient ($a$) that only depends on the true velocity angle ($\phi$). To achieve this goal, Eq. (47a) is used to transform the true velocity into a term depending only on $a$ and $\phi$: $V_t = V_1(1-a)/\sin(\phi)$.

$$ \frac{B}{2} \rho c \dot{r} \frac{V_1^2 (1-a)^2}{(\sin(\phi))^2} [-\cos(\phi)C_L(\alpha(\phi)) + \sin(\phi)C_D(\alpha(\phi))] = \rho \pi 4a (1-a) V_1^2 r \dot{\alpha} \quad (56) $$

So:

$$ \frac{Bc}{8\pi r(\sin(\phi))^2} [-\cos(\phi)C_L(\alpha(\phi)) + \sin(\phi)C_D(\alpha(\phi))] = \frac{a}{1-a}. \quad (57) \\ = \kappa(\phi) $$

Eventually, the axial induction ($a$) relation with the true velocity angle ($\phi$) reads:

$$ a(\phi) = \frac{\kappa(\phi)}{1 + \kappa(\phi)}. \quad (58) $$

## 4.2 Equating torque expressions

The same principle is applied to obtain the tangential induction coefficient ($a'$) as a function of the true velocity angle ($\phi$). To this end, equating torque expressions from Eq. (43) (MT) and Eq. (55) (BET) provides the desired relation.

$$ r \frac{B}{2} \rho c \dot{r} V_t^2 [\sin(\phi)C_L(\alpha(\phi)) + \cos(\phi)C_D(\alpha(\phi))] = 4\pi\Omega\rho V_1 a' (1-a) r^3 \dot{\alpha} \quad (59) $$

Contrary to the previous subsection, the $V_t^2$ term is differently rewritten. On the one hand, the result of Eq. (47a), $V_t = V_1(1-a)/\sin(\phi)$, is injected on the left hand side of Eq. (59). It enables to simplify both $(1-a)$ terms and isolate the tangential induction coefficient ($a'$) on the right hand side. On the other hand, the remaining true velocity term on the left hand side is replaced by the result of Eq. (47b), $V_t = -r\Omega(1+a')/\cos(\phi)$, to eventually obtain a function only depending on $a'$ and $\phi$.

$$ \frac{Bc}{8\pi r^2 \Omega} \frac{-r\Omega(1+a')}{\cos(\phi)} \frac{V_1(1-a)}{\sin(\phi)} [\sin(\phi)C_L(\alpha(\phi)) + \cos(\phi)C_D(\alpha(\phi))] = V_1(1-a)a' \quad (60) $$

Eq. (60) leads to:
$$ -\frac{Bc}{8\pi r \cos(\phi) \sin(\phi)} \underbrace{[\sin(\phi)C_L(\alpha(\phi)) + \cos(\phi)C_D(\alpha(\phi))]}_{=\kappa'(\phi)} = \frac{a'}{1+a'}. \quad (61) $$

The final relation giving the tangential induction ($a'$) as a function of $\phi$ angle is then:

$$ a'(\phi) = -\frac{\kappa'(\phi)}{1 + \kappa'(\phi)}. \quad (62) $$

## 4.3 Solving BEMT equations numerically

### 4.3.1 Immediate iteration process: fixed-point algorithm

A direct iteration process is to guess a starting value for $\phi$. Let this value be $\phi_0$. Then, the axial ($a$, Eq. (58)) and tangential ($a'$, Eq. (62)) inductions are evaluated using this $\phi_0$ value. The obtained inductions provide a true velocity vector through Eq. (46). As a consequence, the associated angle ($\phi$) is assessed using Eq. (47b). The induction coefficients are modified in turn using Eq. (58) and Eq. (62), possibly with a relaxation factor to enhance the convergence process. Those iterations go on and on till reaching a threshold on the difference between $\phi$ angle values of two consecutive iterations. The threshold is user-defined and typically is of a tenth of a percent of the absolute value.

### 4.3.2 Ning's root-finding algorithm

Both $a$ and $a'$ are expressed as a function of $\phi$. For a root-finding process to work, an equation linking $a$, $a'$ and $\phi$ is required. It can be seen as a closing relation. The equation performing this task is Eq. (47c) with, on the left hand side, $\phi$-dependent terms and, on the right hand side, induction-dependent ($a$, $a'$) terms. It reads:

$$ \frac{\sin(\phi)}{\cos(\phi)} = \frac{V_x}{V_\theta} \Leftrightarrow \frac{\sin(\phi)}{\cos(\phi)} = -\frac{V_1(1 - a(\phi))}{r\Omega(1 + a'(\phi))}. \quad (63) $$

According to Ning (2013), Eq. (63) is rewritten to enhance its mathematical convergence.

$$ \frac{\sin(\phi)}{V_1(1 - a(\phi))} + \frac{\cos(\phi)}{r\Omega(1 + a'(\phi))} = \underbrace{\frac{\sin(\phi)}{V_x} - \frac{\cos(\phi)}{V_\theta}}_{=f(\phi)} = 0 \quad (64) $$

So the problem on $a$ and $a'$ reduces to a one-dimensional problem which can be solved with any root-finding algorithm like the Newton-Raphson method. In other words, the question becomes: when does $\phi \mapsto f(\phi)$ equate to zero? When Eq. (64) is solved providing a $\phi$ angle value, $a$ and $a'$ are respectively evaluated by means of Eq. (58) and Eq. (62). However, a limitation of Eq. (64) formulation is that it prevents from accounting for multiple Reynolds number per profile in the polar curves. This is due to re-opening the closing relation by making the induction coefficients ($a$ and $a'$) dependent on themselves, and not only on $\phi$. It is the case because the inductions depend upon the lift and drag coefficients which, if they are dependent upon the local Reynolds number, are dependent upon the inductions to assess the local value of the Reynolds number. A way to circumvent this limitation is to use an estimate of the local velocity without induction to compute the local Reynolds number.
# 5 Example of application

The notations developed in this document may appear unusual to the reader. To help getting familiar with the present BEMT, an example is detailed in this section. It is based on the configuration provided in Sec. 9 of Ingram (2011).

## 5.1 Input data

The studied rotor is composed of three blades. Their geometry is reported in Tab. 1. The hub radius and blade length respectively equate to: $R_h = 0.2$ m and $L_b = 4.8$ m. The original blade geometry (Ingram 2011) can be reconstructed by applying the local radius definition: $r = R_h + r_b$. Furthermore, the twist plus pitch angle ($\gamma$) values given in Tab. 1 agree with the formalism detailed in Sec. 3, specially Fig. 11b and Fig. 12.

<table><caption>Table 1: Blade geometry with local radius ($r_b$), chord ($c$) normalized by blade length ($L_b$) and twist plus pitch angle ($\gamma$).</caption><thead><tr><th>r<sub>b</sub>/L<sub>b</sub> [-]</th><th>c/L<sub>b</sub> [-]</th><th>γ [°]</th></tr></thead><tbody><tr><td>0</td><td>0.146</td><td>-29</td></tr><tr><td>0.167</td><td>0.148</td><td>-15.7</td></tr><tr><td>0.375</td><td>0.092</td><td>-5.1</td></tr><tr><td>0.583</td><td>0.063</td><td>-0.9</td></tr><tr><td>0.792</td><td>0.048</td><td>1.3</td></tr><tr><td>1</td><td>0.04</td><td>2.6</td></tr></tbody></table>

The rotor blade profile is the NACA 0012. Exact two-dimensional geometry is retrieved from Abbott et al. (1959) and Carmichael (2001). The tabulated polar curves are obtained using XFoil software by Drela et al. (1987) and expanded at high angles of attack by means of Viterna et al. (1982) method. The resulting lift ($C_L$) and drag ($C_D$) polar curves are shown in Fig. 13. Some discrepancies are observed compared to the polar presented by Ingram (2011). They can be explained by the uncertainties on parameters crucial to compute the polar curves, namely: Reynolds number (approximately $10^6$) and “critical N-factor” ($N_{crit}$) piloting boundary layer transition, set close to wind tunnel values.

The rotor operates at a tip speed ratio of height: TSR = 8. The inflow velocity is set as: $V_1 = 7 \text{ m·s}^{-1}$. To simplify the solving process, it is now assumed that the drag contribution is zero, meaning $C_D = 0$. It is worth noting that it is only performed to ease the present example. However, drag plays an important role in the correct assessment of rotor performance. It should be accounted for in real-life analysis. Hence, this assumption simplifies $\kappa$ (Eq. (57)) and $\kappa'$ (Eq. (61)) terms respectively used to define axial ($a$, Eq. (58)) and tangential ($a'$, Eq. (62)) inductions:

$$ \kappa(\phi) = -\frac{Bc}{8\pi r} \times \frac{\cos(\phi)C_L(\alpha(\phi))}{(\sin(\phi))^2}, \quad (65) $$

$$ \kappa'(\phi) = \frac{Bc}{8\pi r} \times \frac{\sin(\phi)C_L(\alpha(\phi))}{\cos(\phi)\sin(\phi)}. \quad (66) $$
<figure><img src="image_16.png" /><figcaption>Figure 13: NACA 0012 profile lift (blue, circle) and drag (orange, diamond) coefficients polar curves as function of angle of attack (α). Computed using XFoil and Viterna's method.</figcaption></figure>

## 5.2 Iterations by hand on the last section

To illustrate the BEMT solving process (Sec. 4), iterations are conducted by hand on the tip blade section, corresponding to the last line of Tab. 1. The solving process is based on the fixed-point algorithm as it is easier to handle by hand. No relaxation, meaning a relaxation coefficient of 1, is performed because of the convergence stability in the present case. The process requires knowing the tangential velocity component. To this end, the rotor angular velocity is immediately evaluated using the given operating TSR value:

$$ \Omega = \frac{TSR \times V_1}{R_h + L_b} \approx \frac{8 \times 7}{5} = 11.2 \, \text{rad} \cdot \text{s}^{-1}. \qquad (67) $$

In addition, Eq. (65) and Eq. (66) evidence a common term ($B_c/8\pi r$) depending on the studied section. This term can be interpreted as a quarter of the local solidity which is the ratio of the blade elements surface divided by the surface of the complete annular element. It can be immediately evaluated using Tab. 1 values: $r = 5$ m and $c \approx 0.19$ m.

$$ \frac{B_c}{8\pi r} \approx \frac{3 \times 0.19}{8 \times 3.1415 \times 5} \approx 4.58 \times 10^{-3} \qquad (68) $$

### Initialization

As previously mentioned, the tip blade section is studied. It is where the strongest tangential velocity component occurs. Hence, by remembering Fig. 12 schematic, the true velocity angle can be inferred to obtain a starting point: $\phi_0 = 170^{\circ}$. The associated angle of attack is assessed using Eq. (49) in the turbine case.

$$ \alpha_0 \approx 180 + 2.6 - 170 = 12.6^{\circ} \qquad (69) $$

### Iteration 1

The initialization value of the angle of attack is used to estimate the lift coefficient value from the NACA 0012 polar curves (Fig. 13): $C_L \approx 1.2$. All the parameters needed to assess axial and tangential induction factors have been gathered. So the axial induction of the first iteration ($a_1$) is evaluated using Eq. (65) and Eq. (58):
$$ \kappa_1 \approx -4.58 \times 10^{-3} \times \frac{\cos(170^\circ) \times 1.2}{(\sin(170^\circ))^2} \approx 1.86 \times 10^{-1}, \qquad (70) $$

$$ a_1 \approx \frac{1.86 \times 10^{-1}}{1 + 1.86 \times 10^{-1}} \approx 1.57 \times 10^{-1}. \qquad (71) $$

The considered quantities have no dimension. So no unit is specified. Similarly, the tangential induction of the first iteration ($a'_1$) is assessed with Eq. (66) and Eq. (62).

$$ \kappa'_1 \approx 4.58 \times 10^{-3} \times \frac{\sin(170^\circ) \times 1.2}{\cos(170^\circ) \times \sin(170^\circ)} \approx -5.8 \times 10^{-3} \qquad (72) $$

$$ a'_1 \approx -\frac{-5.8 \times 10^{-3}}{1 + (-5.8 \times 10^{-3})} \approx 5.8 \times 10^{-3} \qquad (73) $$

Axial ($V_x$) and tangential ($V_\theta$) velocity components are evaluated using those induction coefficients and Eq. (46). The norm of the true velocity is: $V_t = \sqrt{V_x^2 + V_\theta^2}$. Hence, as advised in Sec. 3, the true velocity angle ($\phi_1$) is assessed using Eq. (47b):

$$ \phi_1 \approx \arccos \left( \frac{-5 \times 11.2 \times (1 + 5.8 \times 10^{-3})}{56.6} \right) \approx 174^\circ. \qquad (74) $$

The resulting angle of attack for this first iteration is then:

$$ \alpha_1 \approx 180 + 2.6 - 174 = 8.6^\circ. \qquad (75) $$

## Iteration 2

The exact same process is conducted for the second iteration. The lift coefficient is evaluated using the angle of attack of the previous iteration and the NACA 0012 polar curves (Fig. 13): $C_L \approx 0.9$. The true velocity angle from previous iteration ($\phi_1$) is used to compute the axial ($a_2$) and tangential ($a'_2$) inductions of this second iteration.

$$ \kappa_2 \approx -4.58 \times 10^{-3} \times \frac{\cos(174^\circ) \times 0.9}{(\sin(174^\circ))^2} \approx 3.98 \times 10^{-1} \qquad (76) $$

$$ a_2 \approx \frac{3.98 \times 10^{-1}}{1 + 3.98 \times 10^{-1}} \approx 2.85 \times 10^{-1} \qquad (77) $$

$$ \kappa'_2 \approx 4.58 \times 10^{-3} \times \frac{\sin(174^\circ) \times 0.9}{\cos(174^\circ) \times \sin(174^\circ)} \approx -4.4 \times 10^{-3} \qquad (78) $$

$$ a'_2 \approx -\frac{-4.4 \times 10^{-3}}{1 + (-4.4 \times 10^{-3})} \approx 4.4 \times 10^{-3} \qquad (79) $$

$$ \phi_2 \approx \arccos \left( \frac{-5 \times 11.2 \times (1 + 4.4 \times 10^{-3})}{56.5} \right) \approx 174.9^\circ \qquad (80) $$

$$ \alpha_2 \approx 180 + 2.6 - 174.9 = 7.7^\circ \qquad (81) $$
Iteration 3

A third iteration is conducted by hand. The lift coefficient becomes: $C_L \approx 0.85$. Axial ($a_3$) and tangential ($a'_3$) inductions are evaluated using $\phi_2$ angle from previous iteration.

$$ \kappa_3 \approx -4.58 \times 10^{-3} \times \frac{\cos(174.9^\circ) \times 0.85}{(\sin(174.9^\circ))^2} \approx 5.07 \times 10^{-1} \quad (82) $$

$$ a_3 \approx \frac{5.07 \times 10^{-1}}{1 + 5.07 \times 10^{-1}} \approx 3.37 \times 10^{-1} \quad (83) $$

$$ \kappa'_3 \approx 4.58 \times 10^{-3} \times \frac{\sin(174.9^\circ) \times 0.85}{\cos(174.9^\circ) \times \sin(174.9^\circ)} \approx -4.0 \times 10^{-3} \quad (84) $$

$$ a'_3 \approx \frac{-4.0 \times 10^{-3}}{1 + (-4.0 \times 10^{-3})} \approx 4.0 \times 10^{-3} \quad (85) $$

$$ \phi_3 \approx \arccos \left( \frac{-5 \times 11.2 \times (1 + 4.0 \times 10^{-3})}{56.4} \right) \approx 175.3^\circ \quad (86) $$

$$ \alpha_3 \approx 180 + 2.6 - 175.3 = 7.3^\circ \quad (87) $$

At this point, the threshold on the difference between two consecutive values of the true velocity angle ($\phi$) is not reached. However, the reader must now be familiar with the BEMT process. The converged values for all the blade sections are given in the next subsection.

## 5.3 Numerically converged induction values

As observed during the iterations manually performed on the tip blade section, minor changes in true velocity angle ($\phi$) still generate significant changes in associated inductions. This phenomenon is all the more visible on the axial induction ($a$). So, to obtain the numerically converged values presented in Tab. 2, the threshold on the difference between two consecutive values of the true velocity angle ($\phi_k$ and $\phi_{k-1}$) is set as:

$$ \left| \frac{\phi_k - \phi_{k-1}}{\phi_{k-1}} \right| \times 100 \le 0.05. \quad (88) $$

<table><caption>Table 2: Converged induction values ($a$) and ($a'$) for each blade section identified by its normalized radius ($r_b/L_b$). Are also reported the true velocity angle ($\phi$), the angle of attack ($\alpha$), the lift coefficient ($C_L$) and the iteration number reaching convergence ($n_{Ite}$).</caption><thead><tr><th>r<sub>b</sub>/L<sub>b</sub> [-]</th><th>φ [°]</th><th>α [°]</th><th>C<sub>L</sub> [-]</th><th>a [-]</th><th>a' [-]</th><th>n<sub>Ite</sub> [-]</th></tr></thead><tbody><tr><td>0</td><td>125.9</td><td>25.1</td><td>0.63</td><td>0.192</td><td>0.829</td><td>13</td></tr><tr><td>0.167</td><td>157.8</td><td>6.5</td><td>0.77</td><td>0.299</td><td>0.076</td><td>7</td></tr><tr><td>0.375</td><td>168.5</td><td>6.4</td><td>0.77</td><td>0.333</td><td>0.021</td><td>2</td></tr><tr><td>0.583</td><td>172.2</td><td>6.9</td><td>0.81</td><td>0.345</td><td>0.010</td><td>4</td></tr><tr><td>0.792</td><td>174.2</td><td>7.1</td><td>0.82</td><td>0.358</td><td>0.006</td><td>5</td></tr><tr><td>1</td><td>175.6</td><td>7.0</td><td>0.82</td><td>0.385</td><td>0.004</td><td>6</td></tr></tbody></table>
For each blade section, the converged state is reported in Tab. 2. The number of iterations needed to reach the threshold (Eq. (88)) is specified in the last column ($n_{Ite}$). The starting point angle ($\phi_0$) has been kept at 170° for each section. Finally, although the order of magnitude of the values presented in Tab. 2 are well inline with the ones reported by Ingram (2011), some discrepancies are noticed. The values obtained here have been checked multiple times, always converging towards the same results.

# 6 Turbines performance assessment

To get an overview of a turbine behaviour, additional quantities are needed, namely global loads integrated on the rotor. Those integrated loads can be evaluated using the BEMT results. Global loads are then used to define normalized coefficients that characterize the turbine performance for a broad range of operating conditions.

## 6.1 Rotor loads

The BEMT results provide induction coefficients ($a$ and $a'$), thus true velocities ($\vec{V}_t$), for each blade section, meaning for each radius ($r$). Those radius-dependent velocities enable to assess local loads at each blade section. Global loads on the turbine are then obtained by summing the local contributions for each section and for each blade. Consequently, rotor axial force ($F_x$), or thrust, and torque ($Q$) are respectively defined as:

$$F_x = \sum_{k=1}^{B} \sum_{j=1}^{N_s} dF_{x,1}(a_j, a'_j), \qquad (89)$$

$$Q = \sum_{k=1}^{B} \sum_{j=1}^{N_s} dQ_1(a_j, a'_j), \qquad (90)$$

with $dF_{x,1}$ evaluated using Eq. (52) and $dQ_1$ assessed using Eq. (53). Based upon the rotor loads defined by Eq. (89) and Eq. (90), normalized coefficients can be defined to generalize the performance of the turbine to other inflow velocity conditions ($V_1$).

## 6.2 Performance coefficients

The performance coefficients of the turbine pursue two main objectives. The first one is to know the mechanical power produced by the turbine rotor to be able to dimension the electrical and mechanical parts or, conversely, to adapt the wind turbine to the electro-mechanical components. To this end, the turbine power coefficient ($C_P$, Eq. (91)) is defined as the ratio between the power extracted by the turbine and the power carried by the upstream wind in a disc of equivalent surface:

$$C_P = \frac{Q\Omega}{\frac{1}{2}\rho\pi(R_h + L_b)^2(V_1)^3}. \qquad (91)$$

The second objective is to evaluate the axial force, or thrust, to be able to dimension the turbine foundations. Hence, the turbine thrust coefficient ($C_T$, Eq. (92)) is defined as the ratio between the thrust (axial) force exerted by the wind on the turbine and the axial force exerted by the wind on a disc of equivalent surface:
$$C_T = \frac{F_x}{\frac{1}{2}\rho\pi(R_h + L_b)^2(V_1)^2} \qquad (92)$$

## 6.3 Reduced velocity or tip speed ratio

Eventually, the objective is to characterize the turbine performance with respect to as much operating points as possible. Each operating point is characterized by a given tuple made of upstream flow condition ($V_1$) and turbine angular velocity ($\Omega$) It has been shown in the literature (Burton et al. 2011; Schmitz 2013; Sørensen 2016) that power and thrust coefficients of a given turbine remain similar when a parameter called the reduced velocity, or tip speed ratio (TSR), is kept constant. This parameter compares the blades tip tangential velocity to the mean upstream flow velocity, then representing a given operating point. The TSR is defined as follow:

$$TSR = \frac{(R_h + L_b)\Omega}{V_1} \qquad (93)$$

To summarize, a preliminary evaluation of a turbine model consists in producing its characteristic curves of $C_P$ and $C_T$ as function of TSR.

## 7 Conclusion

To summarize, the demonstration of the BEMT equations for horizontal axis rotors is performed. The study focuses on turbines, extracting energy from the fluid, but specific elements dedicated to propellers analysis are also specified. Strategies to solve BEMT equations are shown in this document. The presented formalism, based on direct orthonormal bases, enables to combine the present approach with other frameworks, although leading to slightly cumbersome notations. To get familiar with these notations, an example is detailed showing the applicability of the method and the convergence process. Finally, general aspects on rotor performance analysis are defined. The presented BEMT can be enhanced by using a correction to account for a finite number of blades in the momentum theory (Burton et al. 2011; Sørensen 2016; Branlard 2017) and account for yawed or highly-loaded rotors (Liew et al. 2024).

## References

* Abbott, I. H. and A. E. von Doenhoff (1959). *Theory of Wing Sections: Including a Summary of Airfoil Data*. Correc. rep. of : McGraw-Hill Book, 1949. New York: Dover.

* Betz, A. (1920). "Das Maximum Der Theoretisch Moglichen Auswendung Des Windes Durch Windmotoren". In: *Zeitschrift fur gesamte Turbinewesen* 26.

* Branlard, E. (2017). *Wind Turbine Aerodynamics and Vorticity-Based Methods: Fundamentals and Recent Applications*. Vol. 7. Research Topics in Wind Energy. Cham: Springer International Publishing. DOI: 10.1007/978-3-319-55164-7.

* Burton, T., N. Jenkins, D. Sharpe, and E. Bossanyi (May 2011). *Wind Energy Handbook*. 1st ed. Wiley. DOI: 10.1002/9781119992714.
Carmichael, R. (Oct. 2001). "Algorithm for Calculating Coordinates of Cambered NACA Airfoils at Specified Chord Locations". In: *1st AIAA, Aircraft, Technology Integration, and Operations Forum*. Los Angeles, CA, U.S.A.: American Institute of Aeronautics and Astronautics. DOI: 10.2514/6.2001-5235.

Drela, M. and M. B. Giles (Oct. 1987). "Viscous-Inviscid Analysis of Transonic and Low Reynolds Number Airfoils". In: *AIAA Journal* 25.10, pp. 1347-1355. DOI: 10.2514/3.9789.

Drzewiecki, S. (1920). *Théorie générale de l'hélice: hélices aériennes et hélices marines*. Gauthier-Villars et cie.

Dufour, M.-A. (Oct. 2025). "Surface Wave Effects on Horizontal Axis Tidal Turbines Behaviour: An Experimental and Numerical Study". PhD thesis. Le Havre: Université Le Havre Normandie.

Froude, W. (July 1920). *On the Elementary Relation Between Pitch, Slip, and Propulsive Efficiency*. NACA Technical Memorandum NACA-TM-1. National Advisory Committee for Aeronautics.

Giovannini, A. and C. Airiau (2025). *Aérodynamique fondamentale*. 2e éd. Toulouse: Cépaduès éditions.

Glauert, H. (1935). "Airplane Propellers". In: *Aerodynamic Theory*. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 169-360. DOI: 10.1007/978-3-642-91487-4_3.

Ingram, G. (Oct. 2011). *Wind Turbine Analysis Using the Blade Element Momentum Method*. Tech. rep. Durham, UK: School of Inngineering, Durham University.

Joukowsky, N. E. (1920). "Windmill of the NEJ Type". In: *Transactions of the Central Institute for Aero-hydrodynamics of Moscow* 1, p. 57.

Katz, J. and A. Plotkin (2001). *Low-Speed Aerodynamics*. 2nd Edition. Cambridge University Press.

Lanchester, F. W. (May 1915). "A CONTRIBUTION TO THE THEORY OF PROPULSION AND THE SCREW PROPELLER". In: *Journal of the American Society for Naval Engineers* 27.2, pp. 509-510. DOI: 10.1111/j.1559-3584.1915.tb00408.x.

Ledoux, J., S. Riffo, and J. Salomon (Jan. 2021). "Analysis of the Blade Element Momentum Theory". In: *SIAM Journal on Applied Mathematics* 81.6, pp. 2596-2621. DOI: 10.1137/20M133542X.

Liew, J., K. S. Heck, and M. F. Howland (Aug. 2024). "Unified Momentum Model for Rotor Aerodynamics across Operating Regimes". In: *Nature Communications* 15.1, p. 6658. DOI: 10.1038/s41467-024-50756-5.

Ning, A. S. (July 2013). "A Simple Solution Method for the Blade Element Momentum Equations with Guaranteed Convergence: Guaranteed Solution of the BEM Equations". In: *Wind Energy*. DOI: 10.1002/we.1636.

Rankine, W. J. M. (1865). "On the Mechanical Principles of the Action of Propellers". In: *Transactions of the Institution of Naval Architects* 6.

Schmitz, S. (2013). *AERSP 583 Wind Turbine Aerodynamics*. https://www.e-education.psu.edu/aersp583/. University Website. (Visited on 05/17/2023).

Sørensen, J. N. (2016). *General Momentum Theory for Horizontal Axis Wind Turbines*. Vol. 4. Research Topics in Wind Energy. Cham: Springer International Publishing. DOI: 10.1007/978-3-319-22114-4.

Viterna, L. A. and D. C. Janetzke (Sept. 1982). *Theoretical and Experimental Power from Large Horizontal-Axis Wind Turbines*. Tech. rep. DOE/NASA/20320-41; NASA-TM-82944. NASA Lewis Research Center. DOI: 10.2172/6763041.