# Calcul des efforts aérodynamiques agissant sur les pales d'une petite éolienne

Zine Labidine Mahri *, Mohamed Salah Rouabah et Said Zid

Département de Génie climatique,
Université de Mentouri, Constantine, Algérie

(reçu le 09 Mai 2007 – accepté le 25 Juin 2007)

**Résumé** - La modélisation aérodynamique de l'hélice d'un aérogénérateur est une étape essentielle dans la conception de cette machine. Son objectif est de calculer les efforts aérodynamiques qui s'exercent sur les pales, de déterminer les paramètres optimaux de ces pales et d'estimer la puissance maximale extraite par cette machine. Cette modélisation est une tache complexe à cause de la rotation de l'hélice qui crée un écoulement tridimensionnel de l'air autour des pales. Dans ce travail, deux théories sont combinées, la première théorie est celle de l'écoulement axial (the axial momentum theory), la deuxième est la théorie de l'élément de la pale (blade element theory), qui tient compte de la composante de rotation de l'écoulement de l'air engendré par la rotation des pales. La combinaison de ces deux approches a permis de résoudre les équations des efforts et du couple aérodynamiques par méthode itérative. La convergence de l'algorithme de résolution a été accomplie pour plusieurs profils aérodynamiques.

**Abstract** - The aerodynamic modeling of the wind turbine blades constitutes one of the most important processes in the design of the turbine. Its objective is to compute the aerodynamic loads, to estimate the energetic performance and to determine the blade optimal parameters. This aerodynamic modeling is a complicated task due to airflow around the blades, produced by the rotation of the rotor. In this work, two aerodynamic theories are used, the first one is the axial momentum theory and the second is the blade element theory. In this later theory the effect of wake rotation is included. When both theories are combined the sets of equations obtained by the two approaches can be solved, using iterative method, to obtain the aerodynamic loads and torque. The algorithm used to solve this problem has converged for many aerodynamic profiles.

**Mots clés:** Energie éolienne - Aérodynamique - Analyse numérique - Mécanique des fluides.

## 1. INTRODUCTION

La conception d'un aérogénérateur repose essentiellement sur les théories aérodynamiques; car cette machine utilise des pales caractérisées par une géométrie spéciale ayant un profil aérodynamique. Le dimensionnement de cette forme géométrique a un grand impact sur son rendement énergétique et par conséquent sur sa rentabilité économique.

Pour les petites éoliennes, où le coût du kilowatt reste encore élevé et le site d'implantation est généralement imposé, la modélisation aérodynamique prend une importance particulière. Des calculs très poussés doivent être menés dans le but de déterminer les différents paramètres de la pale, comme les distributions de la longueur de sa corde, l'angle de vrillage ainsi que l'épaisseur le long de cette pale.

La modélisation aérodynamique de la pale est une tache complexe, cette complexité est due à son mouvement rotatif qui crée un écoulement tridimensionnel de l'air autour de cette pale. Cette modélisation a pour objectif de déterminer les paramètres optimaux de la pale (permettant de capter une puissance maximale) et ensuite de calculer les efforts aérodynamiques agissant sur les pales ainsi que la puissance captée.

Ce calcul, qui permet de déterminer les efforts aérodynamiques, est nécessaire aussi bien pour estimer les performances énergétiques de l'aérogénérateur, que pour prévoir les problèmes de structure.

<footnote>* zlmahri@hotmail.com</footnote>
Pour rendre cette modélisation possible, plusieurs simplifications sont à envisager.

En effet, deux théories sont combinées, la première théorie est celle de l'écoulement axial (the axial momentum theory). Dans cette approche, un modèle unidimensionnel est utilisé pour exprimer l'énergie extraite du vent, en supposant que l'air est incompressible. La deuxième est la théorie de l'élément de la pale (blade element theory), qui tient compte de la composante de rotation de l'air due à la rotation des pales.

Si la première théorie ne peut pas révéler seule les informations nécessaires pour dimensionner les pales [1], sa combinaison avec la deuxième théorie permet d'avoir tous les renseignements utiles.

La méthode de calcul développée dans ce travail est basée sur la solution simultanée des équations des efforts et des couples, obtenues par les deux théories.

Dans cette méthode, on commence par donner des valeurs initiales (arbitraire) aux facteurs d'interférences axial et tangentiel. Ces deux facteurs permettront de déterminer l'angle d'écoulement de l'air ainsi que les deux efforts aérodynamiques axial et tangentiel. Ces paramètres permettront à leur tour de calculer de nouveau les deux facteurs d'interférences.

Les étapes de calcul précédentes seront répétées suivant une approche itérative jusqu'à convergence entre les anciennes et les nouvelles valeurs des facteurs d'interférences.

Une formulation détaillée des deux théories est présentée dans les paragraphes suivants.

## 2. THEORIES AERODYNAMIQUES

Dans ce travail, deux théories sont utilisées d'une manière combinée, la première théorie est celle de l'écoulement axial (the axial momentum theory). Dans cette approche, un modèle unidimensionnel est utilisé pour exprimer l'énergie extraite du vent, en supposant que l'air est incompressible.

La deuxième est la théorie de l'élément de la pale (blade element theory), qui tient compte de la composante de rotation de l'air due à la rotation des pales.

### 2.1 Théorie de l'écoulement axial (the axial momentum theory)

Dans cette théorie les suppositions suivantes sont faites [2] :

* L'écoulement est totalement axial (l'air ne subit aucun mouvement de rotation).

* L'écoulement est incompressible.

* La vitesse du vent est constante loin du plan du rotor.

* L'écoulement est rotationellement symétrique.

* L'air passe à travers le rotor sans frottement.

En appliquant le principe de la conservation de masse, on obtient:

$$V_0 A_0 = V A = V_1 A_1 = V_i A_i \quad (1)$$

<figure><img src="image_1.png" /><figcaption>Fig. 1: Ecoulement axial de l'air</figcaption></figure>
où: $V_i$ est la vitesse de l'air à la station i ; $A_i$ est l'aire de la section du flux d'air à la station i .

Aussi en appliquant le théorème de variation de la quantité de mouvement, la force axiale du vent sur le rotor est donnée par l'expression suivante :

$$
T = m(V_0 - V_1) = \rho A V(V_0 - V_1) \quad (2)
$$

Appliquant l'équation de Bernoulli en amont du rotor :

$$
p_0 + \frac{1}{2} \rho V_0^2 = p + \frac{1}{2} \rho V^2 \quad (3)
$$

et celle en aval :

$$
p_0 + \frac{1}{2} \rho V_1^2 = p' + \frac{1}{2} \rho V^2 \quad (4)
$$

En combinant (3) et (4), on obtient :

$$
p - p' = \frac{1}{2} \rho (V_0^2 - V_1^2)
$$

Sachant que :

$$
\begin{align*}
T &= A (p - p') \\
T &= \frac{1}{2} A \rho (V_0^2 - V_1^2) \tag{5}
\end{align*}
$$

Les deux équations (2) et (5) permettent d'écrire :

$$
V = \frac{V_0 + V_1}{2} \tag{6}
$$

La vitesse de l'écoulement de l'air à travers le rotor est la moyenne des deux vitesses celle en amont et celle en aval du rotor.

La puissance extraite du vent est :

$$
p = \frac{1}{2} m (V_0^2 - V_1^2) = \frac{1}{2} \rho V A (V_0^2 - V_1^2) \quad (7)
$$

Le coefficient de puissance C<sub>p</sub> est défini comme le rapport entre la puissance disponible du vent et la puissance extraite :

$$
C_p = \frac{P}{\left(\frac{1}{2}\rho V_0^3\right) A} \quad (8)
$$

Introduisant le facteur d'interférence axial a, défini comme la fraction de diminution de la vitesse du vent, entre celle de l'écoulement libre en amont du rotor et celle traversant le plan du rotor :

$$
V = (1 - a) V_0 \tag{9}
$$

En utilisant l'équation (6) :

$$
V_1 = (1 - 2a) V_0 \tag{10}
$$

En substituant $V_1$ et $V$ dans (5), on obtient :

$$
T = \frac{1}{2} \rho A V_0^2 + 4a(1-a) \quad (11)
$$

$$
P = \frac{1}{2} \rho A V_0^3 + 4a(1-a)^2 \quad (12)
$$

L'expression de $C_p$ devient :
$$C_p = 4a(1-a)^2 \qquad (13)$$

La valeur maximale de théorique du coefficient de puissance, connue sous le nom de limite de Betz, peut être obtenue en prenant la dérivée de $C_p$ par rapport à a égal à zéro :

$$\frac{\partial C_p}{\partial a} = 4(1 - 3a^2) \quad \text{et} \quad \frac{\partial C_p}{\partial a} = 0$$

Cela donne : $a = 1/3$

$$C_p = \frac{16}{27} \approx 0,59259$$

Le rendement maximal pour une éolienne idéale est approximativement 59,3 %. Dans la pratique, ce rendement n'est jamais réalisé à cause des effets suivants [2] :

* L'écoulement de l'air a une composante rotative due à la rotation du rotor.

* La force de traînée n'est jamais nulle à cause des frottements.

* L'hélice contient un nombre fini de pales.

Ce modèle simple unidimensionnel ne décrit pas l'écoulement rotatif de l'air, cependant il permet de définir certains concepts fondamentaux dans le fonctionnent des éoliennes.

Pour que l'hélice extrait l'énergie du vent, la vitesse de ce dernier doit diminuer lorsqu'il traverse l'hélice. Une machine idéale doit ralentir la vitesse du vent de 2/3.

D'après l'équation de continuité (1), la section de l'air balayé par le vent diminue de l'amont vers l'aval. Pour une éolienne idéale, l'aire de la section de l'écoulement à l'amont est de 2/3 celui du rotor, et ce dernier est la moitié de celui de la section en aval. Finalement, l'extraction de toute l'énergie disponible du vent est pratiquement impossible, dans la pratique, le rendement des éoliennes ne dépasse pas les 45 %.

## 2.2 Théorie de l'élément de pale

Dans le modèle précédent, on a supposé que l'écoulement ne subit aucun mouvement de rotation, mais en réalité la loi de conservation du moment cinétique impose que l'air doit avoir un mouvement rotatif afin que le rotor puisse extraire un couple utile. Dans ce cas, le sens de rotation de l'écoulement de l'air est opposé à celui du rotor.

L'introduction du mouvement de rotation de l'air permet à ce modèle de mieux approcher la réalité et d'obtenir des résultats plus fiables.

Dans l'élaboration de ce modèle, les suppositions suivantes sont envisagées [3] :

* L'écoulement en amont loin du plan du rotor, est complètement axial.

* Au niveau du plan du rotor, la vitesse angulaire de rotation de l'air est $\omega$, cette vitesse diminue considérablement loin du rotor, en aval, de telle manière que la pression statique à cet endroit peut être considérée égale à la pression atmosphérique.

* Il n'y a pas d'interférence entre les éléments adjacents de la pale.

* L'écoulement de l'air autour d'un élément de la pale est considéré comme bidimensionnel.

Les expressions de la puissance extraite et du couple appliqué sur les pales, peuvent être obtenues en considérant l'écoulement de l'air à travers un élément annulaire de rayon r, d'épaisseur dr et de section dA (Fig. 2) :

$$dA = 2\pi r dr$$

La force axiale élémentaire est déduite de l'équation (11), en remplaçant la surface A du rotor par la section de l'élément annulaire dA, on obtient dans ce cas :

$$dT = 4a(1-a)\rho V_0^2 \pi r dr \qquad (14)$$
<figure><img src="image_2.png" /><figcaption>Fig. 2: Ecoulement de l'air à travers un élément annulaire</figcaption></figure>

La variation de la quantité de mouvement de l'air dans la direction tangentielle, donne la force tangentielle qui s'exerce sur la pale :

$$dF_t = m dV_t = \rho V dA \omega r \quad (15)$$

$$dF_t = 2\pi\rho V \omega r^2 dr \quad (16)$$

Le couple généré dans la section annulaire est :

$$dQ = 2\pi\rho V \omega r^3 dr \quad (17)$$

Puisque la puissance est :

$$P = Q \Omega$$

La puissance élémentaire extraite à travers la section annulaire est :

$$dP = 2\pi\rho V \omega r^3 dr \quad (18)$$

Où : $\omega$ est la vitesse angulaire de rotation de l'écoulement de l'air dans le plan du rotor.

Le couple total et la puissance totale deviennent :

$$Q = 2\pi\rho \int_{0}^{R} V \omega r^3 dr \quad (19)$$

$$P = 2\pi\rho\Omega \int_{0}^{R} V \omega r^3 dr \quad (20)$$

Pour pouvoir calculer la puissance et le couple (P et Q), il faut déterminer $\omega$. Pour cela on introduit le facteur d'interférence tangentiel $a'$, définit comme suit :

$$\omega = a' \Omega \quad (21)$$

Puisque l'écoulement de l'air en amont, loin du rotor est purement axial, on a : $\omega_0 = 0$.

En appliquant la loi de conservation du moment cinétique entre la section 0 et la section 1, située juste avant le plan du rotor, on obtient :

$$\omega_0 r_0^2 = \omega_1 r_1^2 \quad (22)$$

Cela implique que : $\omega_1 = 0$.

On peut conclure que la vitesse angulaire de l'écoulement de l'air subit une discontinuité au niveau du plan du rotor, cela est dû au couple exercé sur le rotor.
La valeur de la vitesse angulaire de l'air au niveau du plan du rotor est la moyenne de celles situées juste avant et derrière le rotor [4] :

$$ \omega = \frac{\omega_1 + \omega_2}{2} \qquad (23) $$

Donc :

$$ \omega_2 = 2 \Omega a' \qquad (24) $$

Appliquant la loi de conservation du moment cinétique entre les sections respectives 2 et 3, on obtient :

$$ \omega_2 r_2^2 = \omega_3 r_3^2 \qquad (25) $$

Etant donné que : $r_1 \approx r_2 \approx r$

L'équation (25) devient :

$$ \omega_2 r^2 = \omega_3 r_3^2 \qquad (26) $$

Donc :

$$ \omega_3 = \omega_2 \left( \frac{r}{r_3} \right)^2 \qquad (27) $$

$$ \omega_3 = 2 \Omega a' \left( \frac{r}{r_3} \right)^2 \qquad (28) $$

Puisque le rayon du sillage croît en s'éloignant du rotor vers l'aval, l'équation (28) montre que la vitesse de rotation de l'air diminue plus rapidement en s'éloignant dans le même sens. Cela confirme la deuxième supposition citée ci-dessus.

<figure><img src="image_3.png" /><figcaption>Fig. 3: Section d'un élément de la pale à un rayon r</figcaption></figure>

Soit W est la vitesse relative de l'air,

$$ dq = \frac{1}{2} \rho W^2 dA = \frac{1}{2} \rho W^2 c \Omega_1 \qquad (29) $$

On a :

$$ C_L = \frac{dL}{dq} \qquad (30) $$

$$ C_D = \frac{dD}{dq} \qquad (31) $$
L : la force de portance; D : la force de traînée; C₁ : le coefficient de portance; C₉ : le coefficient de traînée.

$$C_x = \frac{dF_x}{dq} \qquad (32)$$

$$C_y = \frac{dF_y}{dq} \qquad (33)$$

dFₓ est la force tangentielle élémentaire qui s'exerce sur une seule pale, et dFᵧ la force axiale élémentaire qui s'exerce également sur une seule pale.

D'après la figure 3, on peut écrire :

$$C_y = C_L \cos \phi + C_D \sin \phi \qquad (34)$$

$$C_r = C_L \sin \phi - C_D \cos \phi \qquad (35)$$

On a :

$$dF_t = B \cdot dF_x$$

$$dT = B \cdot dF_y$$

où B est le nombre de pales.

D'après la figure 3, on a :

$$\alpha = \phi - \beta \qquad (36)$$

$$\tan \phi = \frac{(1 - a) V_0}{(1 + a') \Omega r} \qquad (37)$$

où α est l'angle d'attaque (d'incidence) et β est l'angle de vrillage.

La force axiale et le couple pour un élément de la pale sont :

$$dT = B C_y dq = B C_y \frac{1}{2} \rho W^2 c dr \qquad (38)$$

$$dQ = B C_x dq r = B C_x \frac{1}{2} \rho W^2 c r dr \qquad (39)$$

L'égalité de l'équation (14) avec (38) et celle de l'équation (17) avec (39) donne :

$$\frac{a}{1-a} = \frac{c B C_y}{8 \pi r \sin^2 \phi} \qquad (40)$$

$$\frac{a'}{1-a'} = \frac{c B C_x}{8 \pi r \sin \phi \cos \phi} \qquad (41)$$

D'après la figure 3, on a :

$$W = \frac{V}{\sin \phi} = \frac{V_0 (1-a)}{\sin \phi} \qquad (42)$$

Aussi :

$$W = \frac{\Omega r (1 + a')}{\cos \phi} \qquad (43)$$

Le coefficient de plénitude local (local solidity) est défini comme suit [5] :

$$\sigma = \frac{c B}{2 \pi r} \qquad (44)$$

La résolution des deux équations (40) et (41), donne les expressions de a et de a' respectivement :
$$a = \frac{1}{\frac{4 \sin^2 \phi}{\sigma C_y} + 1} \qquad (45)$$

$$a = \frac{1}{\frac{4 \sin \phi \cos \phi}{\sigma C_r} + 1} \qquad (46)$$

### Remarque

Certaines auteurs [6-8] affirment que le coefficient de traînée doit être omis des équations (34) et (35), cela est justifié par le fait que certains paramètres comme la force axiale, ainsi que la puissance ont été calculés par la première théorie (la théorie axiale) sans tenir compte des frictions, l'égalité de ces expressions obtenues par les deux théories oblige à prendre $C_d = 0$.

Par contre, d'autres auteurs [9], continuent à inclure le terme de la traînée. Sachant qu'avec un nombre fini de pales, la supposition que l'écoulement est rotationnellement symétrique ne tient pas évidemment.

Cela a un effet sur la performance énergétique qui engendre des pertes de performance concentrées au niveau de la pointe (le bout libre) de la pale. Ces pertes sont connues sous le nom 'tip losses'.

Ce phénomène est quantifié par le facteur de réduction de la circulation proposé par Prandtl et défini par la formule suivante :

$$F = \frac{B \Gamma}{\Gamma_{\infty}} = \frac{2}{\pi a r \cos(e^{-f})} \qquad (47)$$

Avec:

$$f = \frac{(B/2) (R - r)}{r \sin \phi} \qquad (48)$$

$\Gamma$ est la circulation effective de la pale au rayon $r$. $\Gamma_{\infty}$ est la circulation d'une hélice ayant un nombre infini de pale. Les facteurs d'interférence, axial a et tangentiels a' de la vitesse du vent deviennent :

$$a = \frac{1}{\frac{4 F \sin^2 \phi}{\sigma C_y} + 1} \qquad (49)$$

$$a' = \frac{1}{\frac{4 F \sin \phi \cos \phi}{\sigma C_x} - 1} \qquad (50)$$

Dans la pratique, cette correction a très peu d'influence sur la prévision de la performance de l'éolienne et elle est souvent négligée [3].

## 3. DIMENSIONNEMENT OPTIMAL DES PALES POUR UNE PUISSANCE MAXIMALE

Pour un élément annulaire de l'air de rayon $r$ et d'épaisseur $d r$, on peut appliquer l'équation de Bernoulli pour déduire la différence élémentaire de pression à travers la section annulaire des pales, sachant que lorsqu'on traverse l'hélice la vitesse angulaire relative de l'air croit de $\Omega$ à $\Omega + \omega_2$, alors que la vitesse axiale reste inchangée:

$$p + \frac{1}{2} \rho \Omega^2 r^2 = p' + \frac{1}{2} \rho (\Omega + \omega_2)^2 r^2$$
$$p - p' = \rho \left( \Omega + \frac{1}{2} \omega_2 \right) \omega_2 r^2 \quad (51)$$

La force axiale sur l'élément annulaire du rotor est :

$$dT = \rho \left( \Omega + \frac{1}{2} \omega_2 \right) \omega_2 r^2 2 \pi r dr \quad (52)$$

Introduisant le facteur d'interférence tangentiel $a'$, dont l'expression est donnée par l'équation (24) :

$$a' = \frac{\omega_2}{2\Omega}$$

L'expression de la force axiale sur l'élément annulaire, peut s'écrire comme suit :

$$dT = 4 a' (1 + a') \rho \Omega^2 r^2 \pi r dr \quad (53)$$

La division des deux équations (14) et (53) peut donner :

$$\frac{a(1-a)}{a'(1-a')} = \frac{\Omega^2 r^2}{V_0^2} \quad (54)$$

Prenant la vitesse spécifique (tip speed ratio) comme :

$$\lambda = \frac{\Omega R}{V_0}$$

et la vitesse spécifique locale comme [5] :

$$\lambda_r = \frac{\Omega r}{V_0}$$

L'équation (54) devient :

$$\frac{a(1-a)}{a'(1-a')} = \lambda_r^2 \quad (55)$$

L'expression de la puissance extraite dans l'équation (20) peut s'écrire :

$$P = 2\pi \rho V_0 \Omega^2 \int_0^R (1-a) a' r^3 dr \quad (56)$$

Ou bien :

$$P = \frac{1}{2} \pi \rho R^2 V_0^3 \frac{4 \lambda^2}{R^4} \int_0^R (1-a) a' r^3 dr \quad (57)$$

Le coefficient de puissance $C_p$ est formulé comme suit :

$$C_p = \frac{4 \lambda^2}{R^4} \int_0^R (1-a) a' r^3 dr \quad (58)$$

Le coefficient de puissance dépend du facteur K, tel que :

$$K = a' (1 - a) \quad (59)$$

Ce coefficient atteint sa valeur maximale $C_{pmax}$ si l'expression K est maximale.

Les deux équations (55) et (59) permettent d'exprimer K en fonction de a.

Le maximum de K est atteint si, on prend :
$$ \frac{\partial K}{\partial a} = 0 \qquad (60) $$

Cela implique :

$$ \frac{(1-a)(4a-1)^2}{(1-3a)} = \lambda_r^2 \qquad (61) $$

L'égalité des deux équations (55) et (61) permet d'écrire :

$$ a' = \frac{(1 - 3a)}{(4a - 1)} \qquad (62) $$

Cette dernière formule donne la relation entre $a$ et $a'$ afin que le coefficient de puissance soit maximal $C_{pmax}$, cela implique que la puissance extraite par le rotor dans ce cas soit maximale.

## 3.1 Calcul de la forme optimale de la pale

La conception d'une forme optimale de la pale d'une hélice implique que la dernière relation (62), correspondante à une puissance maximale, soit satisfaite. Si on adopte la thèse de négliger les frottements dans le sens axial, en prenant $C_D = 0$, l'expression de $a$ de l'équation (45) devient :

$$ a = \frac{1}{\frac{4 \sin^2 \phi}{\sigma C_L \cos \phi} + 1} \qquad (63) $$

Aussi, l'expression de $a'$ de l'équation (46) devient :

$$ a' = \frac{1}{\frac{4 \cos \phi}{\sigma C_L} - 1} \qquad (64) $$

où $C_L$ est la valeur optimale du coefficient de portance [10].

Si on substitute l'expression (62) dans cette dernière l'équation, on obtient :

$$ a = \frac{4 \cos \phi}{\sigma C_L + 12 \cos \phi} \qquad (65) $$

L'égalité des deux expressions de $a$ dans les équations (63) et (65) donne l'équation quadratique suivante, dont l'inconnue est $\sigma C_L$ :

$$ (\sigma C_L)^2 + 8 \cos \phi \sigma C_L - 16 \sin^2 \phi = 0 \qquad (66) $$

Cette équation admet deux solutions :

$$ \sigma C_L = -4(\cos \phi + 1) \qquad (67) $$

$$ \sigma C_L = 4(1 - \cos \phi) \qquad (68) $$

La première racine n'est pas acceptable, car $\sigma C_L < 0$. La deuxième racine donnée par (68) donne l'expression optimale de la corde, le long de la pale :

$$ C = \frac{8\pi r}{BC_L} (1 - \cos \phi) \qquad (69) $$

Pour déterminer l'angle optimal $\phi$ et par conséquent l'angle optimal de vrillage $\beta$, des transformations similaires sont à envisager. La formule (37) peut s'écrire comme suit :
$$ \tan \phi = \frac{(1-a)}{(1+a')} \frac{1}{\lambda_r} \qquad (70) $$

Si on substitue l'expression (62) dans cette dernière l'équation, on obtient :

$$ \lambda_r = \frac{(4a-1)(1-a)}{a \tan \phi} \qquad (71) $$

La substitution de l'équation (63) dans (71) donne :

$$ \lambda_r = \frac{(2 \cos \phi - 1)(1 + \cos \phi)}{(1 + 2 \cos \phi) \sin \phi} \qquad (72) $$

Cette équation est simplifiée comme suit :

$$ \lambda_r = \frac{\cos \phi + \cos 2\phi}{\sin \phi + \sin 2\phi} \qquad (73) $$

D'où :

$$ \lambda_r = \frac{1}{\tan \frac{3}{2}\phi} \qquad (74) $$

On aboutit finalement à :

$$ \phi = \frac{2}{3} \tan^{-1} \frac{1}{\lambda_r} \qquad (75) $$

On déduit l'angle optimal de vrillage β de l'équation (36) :

$$ \beta = \phi - \alpha_{\text{opt}} \qquad (76) $$

où $\alpha_{\text{opt}}$ est l'angle d'attaque optimale, qui donne $(C_1 / C_d)_{\text{max}}$.

Les résultats de la distribution optimale de la corde, ainsi que l'angle de vrillage pour deux pales de 5 mètres de diamètre et ayant différents profils sont résumés par les tableaux suivants :

<table><caption>**Tableau 1:** Distribution optimale de la corde et de l'angle de vrillage (Profil NACA 63-421)</caption><thead><tr><td>r/R</td><td>Corde (m)</td><td>Angle de vrillage (°)</td></tr></thead><tbody><tr><td>0,16</td><td>0,60</td><td>19,3</td></tr><tr><td>0,20</td><td>0,53</td><td>15,2</td></tr><tr><td>0,24</td><td>0,47</td><td>12,1</td></tr><tr><td>0,29</td><td>0,42</td><td>9,7</td></tr><tr><td>0,33</td><td>0,37</td><td>7,9</td></tr><tr><td>0,37</td><td>0,34</td><td>6,4</td></tr><tr><td>0,41</td><td>0,31</td><td>5,3</td></tr><tr><td>0,45</td><td>0,28</td><td>4,3</td></tr><tr><td>0,50</td><td>0,26</td><td>3,4</td></tr><tr><td>0,54</td><td>0,24</td><td>2,7</td></tr><tr><td>0,58</td><td>0,23</td><td>2,1</td></tr><tr><td>0,62</td><td>0,21</td><td>1,6</td></tr><tr><td>0,66</td><td>0,20</td><td>1,1</td></tr><tr><td>0,71</td><td>0,19</td><td>0,7</td></tr><tr><td>0,75</td><td>0,18</td><td>0,3</td></tr><tr><td>0,79</td><td>0,17</td><td>0,0</td></tr><tr><td>0,83</td><td>0,16</td><td>-0,3</td></tr><tr><td>0,87</td><td>0,15</td><td>-0,6</td></tr><tr><td>0,92</td><td>0,15</td><td>-0,8</td></tr><tr><td>0,96</td><td>0,14</td><td>-1,0</td></tr><tr><td>1,00</td><td>0,13</td><td>-1,2</td></tr></tbody></table>
<table><caption>**Tableau 2:** Distribution optimale de la corde et de l'angle de vrillage (Profil LS1-0417)</caption><thead><tr><td>r/R</td><td>Corde (m)</td><td>Angle de vrillage (°)</td></tr></thead><tbody><tr><td>0,16</td><td>0,82</td><td>23,3</td></tr><tr><td>0,20</td><td>0,73</td><td>19,2</td></tr><tr><td>0,24</td><td>0,64</td><td>16,1</td></tr><tr><td>0,29</td><td>0,57</td><td>13,7</td></tr><tr><td>0,33</td><td>0,51</td><td>11,9</td></tr><tr><td>0,37</td><td>0,46</td><td>10,4</td></tr><tr><td>0,41</td><td>0,42</td><td>9,3</td></tr><tr><td>0,45</td><td>0,39</td><td>8,3</td></tr><tr><td>0,50</td><td>0,36</td><td>7,4</td></tr><tr><td>0,54</td><td>0,33</td><td>6,7</td></tr><tr><td>0,58</td><td>0,31</td><td>6,1</td></tr><tr><td>0,62</td><td>0,29</td><td>5,6</td></tr><tr><td>0,66</td><td>0,27</td><td>5,1</td></tr><tr><td>0,71</td><td>0,26</td><td>4,7</td></tr><tr><td>0,75</td><td>0,24</td><td>4,3</td></tr><tr><td>0,79</td><td>0,23</td><td>4,0</td></tr><tr><td>0,83</td><td>0,22</td><td>3,7</td></tr><tr><td>0,87</td><td>0,21</td><td>3,4</td></tr><tr><td>0,92</td><td>0,20</td><td>3,2</td></tr><tr><td>0,96</td><td>0,19</td><td>3,0</td></tr><tr><td>1,00</td><td>0,18</td><td>2,8</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 4: Distribution optimale de la corde</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5: Angle optimal de vrillage</figcaption></figure>

La distribution calculée pour les deux paramètres, la corde et l'angle ne sont pas linéaires. Dans la pratique, ces distributions ne sont pas respectées à cause des difficultés de fabrication.

Le compromis serait d'utiliser des distributions linéaires proches de celles calculées sans trop sacrifier la performance énergétique (Fig. 4 et 5).

## 3.2 Calcul des efforts aérodynamiques sur le rotor

Une fois que la forme optimale des pales est déterminée à l'aide des formules de la section précédente, il est possible maintenant, de calculer les efforts et les moments agissant sur les pales ainsi que la puissance extraite par l'hélice.

Afin d'effectuer ce calcul, on suppose que les données suivantes sont disponibles :

* Le rayon de l'hélice R.

* La vitesse spécifique $\lambda$.

* Le nombre de pales B.

* Les caractéristiques du profil $C_L(\alpha)$ et $C_D(\alpha)$.

Les données suivantes sont calculées par la méthode décrite précédemment :

* La distribution de la corde $c(r)$.

* La distribution de l'angle de vrillage $\beta(r)$.
La méthode utilisée consiste à calculer les valeurs des facteurs d'interférence axial et tangentielle pour chaque position r de la pale. Puisqu'il n'existe pas d'expressions analytiques pour ces facteurs, une méthode itérative est utilisée, elle se résume par les étapes suivantes :

1. On donne la valeur de r qui détermine la position (la station) sur la pale par rapport au centre de l'hélice, où ce calcul sera effectué.

2. On suppose des valeurs initiales raisonnables pour a et a' (a = 1/3 et a' = 0 à titre d'exemple).

3. On calcule $\phi$ à l'aide de l'équation (37).

4. On calcule $\alpha$ à l'aide de la formule $\alpha = \phi - \beta$ de l'équation (36).

5. On détermine $C_L(\alpha)$ et $C_D(\alpha)$ pour l'angle calculé $\alpha$ par une méthode d'approximation en utilisant le tableau des caractéristiques du profil utilisé.

6. On calcule a et a' en utilisant les formules (49) et (50).

7. On compare les valeurs de a et a' avec celles précédentes et on répète les étapes 3 à 6 jusqu'à l'obtention de la précision voulue.
La convergence de a et a' veut dire que l'égalité des résultats obtenus par les deux théories est vérifiée.

8. On calcule $dC_p$, $dQ$, $dT$ et $dP$.

9. Si cette approche est répétée pour chaque position r de la pale, cela permettra de calculer la valeur totale de $C_p$, Q, P et T par intégration numérique le long de la pale.

Si les pertes de pression au bout libre de la pale (tip losses) sont prises en considération, les formules adéquates doivent être utilisées. Aussi, pour des valeurs de $\alpha$ où le décrochage de la pale peut avoir lieu, des solutions multiples peuvent en résulter.

L'organigramme de la figure 6 résume cet algorithme de résolution.

<figure><img src="image_6.png" /><figcaption>Fig. 6: Organigramme de la méthode BEM, pour le calcul statique et énergétique du rotor</figcaption></figure>
<table><caption>**Table 3:** Distribution des charges aérodynamiques, vitesse du vent 15 m/s, (profil NACA 63-421)</caption><thead><tr><th>Station (r/R)</th><th>Force axiale</th><th>Force tangentielle</th><th>Moment (N.m)</th></tr></thead><tbody><tr><td>0.16</td><td>86,02</td><td>221,24</td><td>206,30</td></tr><tr><td>0.20</td><td>85,93</td><td>286,49</td><td>263,25</td></tr><tr><td>0.25</td><td>81,92</td><td>351,16</td><td>305,56</td></tr><tr><td>0.29</td><td>73,54</td><td>411,09</td><td>323,48</td></tr><tr><td>0.34</td><td>73,37</td><td>466,19</td><td>372,82</td></tr><tr><td>0.38</td><td>38,41</td><td>535,52</td><td>316,39</td></tr><tr><td>0.43</td><td>57,87</td><td>586,57</td><td>467,49</td></tr><tr><td>0.47</td><td>208,88</td><td>556,61</td><td>1008,82</td></tr><tr><td>0.51</td><td>35,67</td><td>764,83</td><td>724,62</td></tr><tr><td>0.56</td><td>304,36</td><td>649,67</td><td>1690,88</td></tr><tr><td>0.60</td><td>39,33</td><td>908,37</td><td>1120,22</td></tr><tr><td>0.65</td><td>354,32</td><td>759,48</td><td>2370,09</td></tr><tr><td>0.69</td><td>83,60</td><td>998,27</td><td>1686,21</td></tr><tr><td>0.73</td><td>320,40</td><td>901,37</td><td>2780,43</td></tr><tr><td>0.78</td><td>221,50</td><td>1012,80</td><td>2591,90</td></tr><tr><td>0.82</td><td>217,84</td><td>1063,01</td><td>2693,02</td></tr><tr><td>0.87</td><td>211,31</td><td>1109,45</td><td>2746,15</td></tr><tr><td>0.91</td><td>191,81</td><td>1148,48</td><td>2624,72</td></tr><tr><td>0.96</td><td>169,79</td><td>1181,84</td><td>2434,31</td></tr><tr><td>1.00</td><td>140,04</td><td>1206,90</td><td>2100,65</td></tr></tbody></table>

Force axiale totale sur pale = 553,24 N

Force tangentielle totale sur une pale = 4729,56 N

Force axiale totale sur l'hélice = 1659,724 N

Moment = 8218,81 N.m.

La figure 7 donne la distribution de la force axiale, le long de la pale.

<figure><img src="image_7.png" /><figcaption>Fig. 7: Distribution de l'effort axial le long de la pale</figcaption></figure>

## 4. CONCLUSION

Les efforts aérodynamiques obtenus serviront aussi bien pour estimer les performances énergétiques aussi bien que pour analyser la structure de l'hélice.
La méthode itérative utilisée, a pu atteindre la convergence en dépit du volume de calcul impliqué dans la formulation du problème. Plusieurs tentatives étaient nécessaires pour trouver le bon choix des valeurs initiales des facteurs d'interférence, qui assurent la convergence la méthode.

L'un des obstacles qui entrave l'utilisation de la méthode itérative, est la difficulté de déterminer les coefficients de portance et de traînée à cause de la variation de l'angle d'incidence à chaque itération. Ce problème a été surmonté à l'aide d'une approximation mathématique des valeurs (obtenues expérimentalement) des coefficients de portance et de traînée en fonction de l'angle d'incidence.

## NOMENCLATURE

<table><tr><td>A<sub>i</sub></td><td>Aire de la section du flux d'air à la station i</td><td>Q</td><td>Couple aérodynamique</td></tr><tr><td>a</td><td>Facteur d'interférence axial</td><td>T</td><td>Force axiale du vent sur le rotor</td></tr><tr><td>a'</td><td>Facteur d'interférence tangentiel</td><td>V<sub>i</sub></td><td>Vitesse de l'air à la station i</td></tr><tr><td>B</td><td>Nombre de pales</td><td>W</td><td>Vitesse relative de l'air</td></tr><tr><td>C<sub>l</sub></td><td>Coefficient de portance</td><td>α</td><td>Angle d'incidence</td></tr><tr><td>C<sub>L</sub></td><td>Valeur optimale du coefficient de portance</td><td>β</td><td>Angle optimal de vrillage</td></tr><tr><td>C<sub>d</sub></td><td>Coefficient de traînée</td><td>λ</td><td>Vitesse spécifique (Tip speed ratio)</td></tr><tr><td>C<sub>P</sub></td><td>Coefficient de puissance</td><td>λ<sub>r</sub></td><td>Vitesse spécifique locale</td></tr><tr><td>D</td><td>Force de traînée</td><td>Ω</td><td>Vitesse de rotation du rotor</td></tr><tr><td>F</td><td>Facteur de réduction de la circulation</td><td>ω</td><td>Vitesse angulaire de rotation de l'écoulement de l'air</td></tr><tr><td>dF<sub>x</sub></td><td>Force tangentielle élémentaire qui s'exerce sur une seule pale</td><td>Γ</td><td>Circulation effective de la pale au rayon r</td></tr><tr><td>dF<sub>y</sub></td><td>Force axiale élémentaire qui s'exerce sur une seule pale</td><td>Γ<sub>∞</sub></td><td>Circulation d'une hélice ayant un nombre infini de pale</td></tr><tr><td>L</td><td>Force de portance</td><td>φ</td><td>Angle d'écoulement</td></tr><tr><td>P</td><td>Puissance extraite du vent</td><td>σ</td><td>Coefficient de plénitude local(local solidity)</td></tr></table>

## REFERENCES

[1] E. Lysen, 'Introduction to Wind Energy', Netherlands, Amersfort, 2nd Edition, 1983.

[2] J.M. Jonkman, 'Modelling of the UAE Wind Turbine for Refinement of FAST_AD', National Renewable Energy Laboratory, Task N°WER3 2010, NREL/TP-500-34755, December 2003.

[3] D. Wood, 'The Design and Analysis of Small Wind Turbines', University of Newcastle, School of Engineering, NSW 2308 Australia, 1st Edition of Notes, 2002.

[4] S.M. Habalia and I.A. Saleh, 'Local Design, Testing and Manufacturing of Small Mixed Airfoil Wind Turbine Blades of Glass Fiber Reinforced Plastics, Part I: Design of the Blade and Root', Journal of Energy Conversion & Management, Vol. 41, pp. 249-280, 2000.
[5] M.M. Duquette and K.D. Visser, 'Numerical Implications of Solidity and Blade Number on Rotor Performance of Horizontal-Axis Wind Turbines', Journal of Solar Energy Engineering, Vol.125, p. 425, Nov. 2003.

[6] S.M. Habali and I.A. Saleh, 'Design and Testing of Small Mixed Airfoil Wind Turbine Blades', Renewable Energy, Vol. 6, N°2, pp. 161 - 169, 1995.

[7] O. Vries, 'Fluid Dynamic Aspects of Wind Energy Conversion', Agard Publication, AG 243, 1979.

[8] W.M.A. Jansen, 'Horizontal Axis Fast Running Wind Turbines for Developing Countries', Publication SWD 76-3, Netherlands, Amersfoort, June 1976.

[9] R.T. Griffiths and M.G. Woollard, 'Performance of Optimal Wind Turbine', Applied Energy, Vol.4, pp. 261 - 272, 1978.

[10] F. Bertagnolio and N. Sorensen, 'Wind Turbine Airfoil Catalogue', Riso (National Laboratory, Roskilde), Denmark, August 2001.