# Prédiction numérique de l'écoulement turbulent autour des pales d'une éolienne à axe horizontal

Rachid Hamel<fnref n="1" /> et Arezki Smaïli<fnref n="2" />*

(reçu le 10 Octobre 2009 – accepté le 20 Novembre 2010)

**Résumé** - Ce papier présente une méthode de simulation numérique pour la prédiction d'écoulement autour d'une éolienne à axe horizontal opérant en régime stationnaire, en utilisant le code 'Fluent'. L'objectif principal est de déterminer les caractéristiques aérodynamiques de la pale en état d'arrêt et en rotation. Les simulations ont été effectuées en considérant les équations de Navier-Stokes moyennées par la procédure de Reynolds. Un maillage non-structuré a été adopté et en utilisant la technique 'moving reference frame'. Le modèle de turbulence Spalart-Allmaras a été utilisé pour décrire l'écoulement turbulent. Les résultats de simulations ont été comparés avec ceux du l'expérimental, pour les conditions à la fois fixe et en rotation. En général, de bonnes concordances ont été notées.

**Abstract** - This work presents a numerical simulation method for prediction of flow around a horizontal axis turbine in steady state condition, using 'Fluent' code. The objective is to determine aerodynamics characteristics for both parked and rotating blade. Simulations are carried out by considering Reynolds averaged Navier-Stokes equations. The rotating reference frame technique was used wherein the blade is fixed in relation to the rotating frame. Turbulent flows with an S-A turbulence model were set as study cases. Computational results compared well with field experimental data, for both parked and rotating blade conditions. In general, good agreements have been noted.

**Mots clés:** Turbine éolienne - Modélisation de la turbulence - CFD - Aérodynamique.

## 1. INTRODUCTION

Le but ultime de ce travail est de développer une méthode de CFD permettant de simuler l'écoulement turbulent autour des pales d'une éolienne à axe horizontal et ce, afin de prédire convenablement les effets tridimensionnels et de la turbulence et par conséquent d'améliorer la précision des prédictions numériques.

Les modèles existants et utilisés jusqu'à maintenant pour le design et la conception des pales sont basées sur des approches simplifiées, c'est-à-dire le concept du disque actuateur et la théorie d'élément de pale [1]. Or, ces modèles s'avèrent être limités à des écoulements stationnaires et uniformes.

Pour tenir compte des effets instationnaires et tridimensionnels, des modèles empiriques sont en outre introduits pour corriger les données aérodynamiques bi dimensionnelles des profils. Le champ d'écoulement autour des pales sera décrit par les équations de Navier-Stokes moyennées par la procédure de Reynolds [2, 3].

<footnote n="1">¹ Laboratoire de Mécanique des Fluides, Unité d'Enseignement & Recherche en Mécanique Appliquée
Ecole Militaire Polytechnique, P.B. 17, Bordj El Bahri, Alger, Algérie</footnote>

<footnote n="2">² Laboratoire de Génie Mécanique et Développement, Ecole Nationale Polytechnique
Avenue Hassen Badi, P.B. 182, El-Harrach, Alger, Algérie</footnote>

<footnote>* arezki.smaili@enp.edu.dz</footnote>
Un modèle de turbulence sera utilisé, pour la fermeture des équations de l'écoulement turbulent moyen. Les équations du modèle mathématique résultant seront implantées dans Fluent. Due à la complexité d'écoulement autour des pales, réaliser une telle étude en utilisant des outils de CFD les plus avancés constitue un grand défi à relever.

En plus, en adoptant une méthode de CFD pour le design et la conception des pales, les retombées économiques et industrielles seraient considérables, puisque les prédictions numériques sont prévues d'être beaucoup plus précises, fiables et robustes.

Dans le présent travail, des simulations autour d'un rotor type NREL phase VI seront présentées. Les résultats de simulations seront comparés avec ceux de la littérature. Nous présenterons enfin nos conclusions et recommandations.

## 2. MODELE MATHEMATIQUE

Pour un écoulement incompressible, les équations régissant la conservation de la masse et la quantité de mouvement sont (moyennées par la procédure de Reynolds) données respectivement par:

Equation de la continuité

$$ \frac{\partial u_i}{\partial x_i} = 0 \qquad (1) $$

Equation de quantité de mouvement

$$ \rho \left( \frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} \right) = - \frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left( \tau_{ij} \right) \qquad (2) $$

$$ \text{Où } \tau_{ij} = \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) $$

### 2.1 Le modèle Spalart Allmaras

Le modèle de Spalart-Allmaras est un modèle relativement simple à une équation. Il s'agit de l'équation de transport de la viscosité cinématique turbulente.

Ceci représente une classe relativement nouvelle de modèles à une équation dans lesquels il n'est pas nécessaire de calculer l'échelle de longueur liée à l'épaisseur locale de la couche de cisaillement.

Ce modèle a été conçu spécifiquement pour les applications de l'aérospatiale comportant des écoulements confinés.

### 2.2 Equation de transport pour le modèle de Spalart Allmaras

La variable transportée, dans le modèle de Spalart-Allmaras, $\tilde{v}$ est identique à la viscosité cinématique turbulente, sauf aux voisinages de la paroi (couches affectées par la viscosité).

L'équation de transport pour $\tilde{v}$ est:
$$ \frac{\partial}{\partial t}(\rho \tilde{v}) + \frac{\partial}{\partial x_i}(\rho \tilde{v} u_i) = G_v + \frac{1}{\sigma_{\tilde{v}}} \left[ \frac{\partial}{\partial x_i} \left\{ (\mu + \rho \tilde{v}) \frac{\partial \tilde{v}}{\partial x_j} \right\} + C_{b2} \rho \left( \frac{\partial \tilde{v}}{\partial x_j} \right)^2 \right] - Y_v + S_v \quad (3) $$

$G_v$ et $Y_v$ représentent respectivement la production et la destruction de la viscosité turbulente produite dans la région proche des parois. $\sigma_{\tilde{v}}$ et $C_{b2}$ sont des constantes et $\mu$ est la viscosité dynamique. $S_{\tilde{v}}$ est un terme source défini par l'utilisateur.

# 3. METHODE NUMERIQUE

## 3.1 Discrétisation de système des équations

La discrétisation des équations de Navier-Stokes et de conservation de la masse est faite sur des éléments triangulaires à cause de la complexité de la géométrie. Elle conduit à un système d'équations non-linéaires couplées. Ce système est résolu de manière itérative.

L'algorithme SIMPLE est employé pour résoudre le couplage vitesse-pression. Le code utilisé pour les simulations de notre travail est FLUENT. Il est basé sur la méthode des volumes finis pour résoudre les équations régissant les écoulements des fluides. GAMBIT a été employé pour la génération du maillage.

## 3.2 Domaine de calcul et maillage

En négligeant les effets de la tour, de la nacelle et du sol, nous avons considéré deux configurations du domaine de calcul.

### 3.2.1 Configuration libre

Le domaine de calcul se compose d'un sous domaine qui englobe la pale du rotor (Fig. 1) où le maillage est plus raffiné et un grand domaine sous forme d'une hémisphérique (Fig. 2). Le rayon de ce dernier est de 12 R (R représente la longueur de la pale) et il est utilisé pour la simulation de la pale fixe et en rotation. Le maillage généré se compose de 1.518.306 éléments.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Sous domaine englobant la pale</figcaption></figure>
<figure><img src="image_2.png" /><figcaption>Fig. 2: Domaine de calcul</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3: Domaine de configuration libre avec le maillage</figcaption></figure>

### 3.2.2 Configuration tuyère

Le domaine de calcul est représenté par un demi-cylindre qui a la forme d'une soufflerie de rayon égal à 16.86 m.

<figure><img src="image_4.png" /><figcaption>Fig. 4: Domaine de configuration tuyère</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 5: Domaine de configuration tuyère avec le maillage</figcaption></figure>

## 4. RESULTATS

Le profil utilisé pour la réalisation de la pale est le S809, utilisé par le NREL dans le développement de leur rotor à deux pales NREL Phase VI.

<figure><img src="image_6.png" /><figcaption>Fig. 6: Pale dans la position 12 heures et la position 47 % de R</figcaption></figure>

### 4.1 Pale fixe

La figure 6 montre la pale d'éolienne dans la position 12 heures. On voit également la position 47 % de R, que l'on va prendre comme référence pour l'angle d'attaque ($\alpha_{47}$).

A $\alpha_{47} = 8.2°$ (Fig. 7), où l'angle d'attaque varie de $\alpha = -1,4°$ à $30°$ du rayon de la pale à $14,4°$ à $95°$ du rayon, l'accord est presque bon pour toutes les positions radiales avec une surestimation pour le coefficient normal $C_N$ et une sous-estimation pour le coefficient tangentiel $C_T$.

On note ici que l'écoulement est attaché à la paroi, qui explique le faible écart entre les coefficients calculés et mesurés à faible angle d'attaque.
<figure><img src="image_7.png" /><figcaption>Fig. 7: Coefficients de force normale et tangentielle pour la pale fixe</figcaption></figure>

Les distributions de pression

<figure><img src="image_8.png" /></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8: Comparaison de la distribution de pression avec l'expérimental pour la pale fixe</figcaption></figure>

Les résultats de la comparaison des répartitions de pression entre la simulation et de l'expérimental sont présentés dans la figure 8 qui montre que l'écoulement est le plus souvent attaché à la paroi et dans ce cas de la pale fixe, nous avons une bonne prédiction des distributions de pression, qui est justifié par les courbes C<sub>N</sub> et C<sub>T</sub>.

## 4.2 Pale tournante

En examinant les distributions des coefficients des forces normales et tangentielles dans le sens de l'envergure pour la pale tournante, nous remarquons un bon accord, comme le montre la figure 9 avec des écarts remarquables pour les coefficients de force tangentielle.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Coefficients de force normale et tangentielle pour la pale tournante</figcaption></figure>

## Les distributions de pression

En examinant la figure 10, il est évident que les mesures montrent une séparation de l'écoulement, alors que les calculs de préserver un maximum d'aspiration.
Un examen plus attentif à la figure 10 montre que les résultats de la configuration du libre sont les plus proche des mesures.

## 5. COURBES DE PUISSANCE

Les courbes puissance-vitesse du vent des éoliennes indiquent la puissance fournie en fonction de la vitesse du vent.

Pour concevoir et évaluer le rendement du rotor d'éolienne, un code informatique a été développé, basé sur la théorie de l'élément de pale, et il a été mis en œuvre après validation.

Nous avons obtenu les résultats présentés dans le graphhe illustré sur la figure 11.

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

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

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

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 10: Comparaison de la distribution de pression avec l'expérimental pour $U_∞ = 10$ m/s</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 11: Courbes de la puissance</figcaption></figure>

On remarque un bon accord entre les résultats de calcul et de l'expérimental pour les faibles vitesses de vent où nous trouvons une faible erreur entre eux, surtout pour la configuration libre, pour $U_∞ = 7, 10, 13$ m/s. Mais l'écart augmente et devient plus important en fonction de l'augmentation de la vitesse du vent.

# 6. CONCLUSION

Dans ce papier, nous avons présenté deux travaux.

Le premier consiste en l'étude de l'écoulement autour de la pale fixe, et le deuxième consiste en l'étude de l'écoulement autour d'une pale en rotation.

Pour ce faire, nous avons proposé une méthode numérique capable d'effectuer des simulations avec précision acceptable et un temps de calcul adéquat. Les résultats
montrent qu'il y a de bonnes corrélations entre les valeurs calculées et les données expérimentales de la NREL phase VI notamment pour les faibles vitesses de vent.

Avec l'augmentation de la vitesse, les puissances données par le modèle BEM sont moins prédictives. Cette divergence est probablement due à l'entrée dans une zone de fonctionnement en décrochage.

## REFERENCES

[1] J.F. Manwell, J.G. McGowan and A.L. Rogers, ‘Wind Energy Explained: Theory, Design and Application’, John Wiley & Sons Ltd, 590 p., 2002.

[2] J. Johansen, N.N. Sørensen and J.A. Michelsen and S. Schreck, ‘Detached-Eddy Simulation of Flow Around the NREL Phase VI Blade’, Wind Energy, Vol. 5, N°2-3, pp. 185 - 197, 2002.

[3] N.N. Sørensen, J.A. Michelsen and S. Schreck, ‘Navier-Stokes Predictions of the NREL Phase VI Rotor in the NASA Ames 80ft x 120ft Wind Tunnel’, Wind Energy, Vol. 5, N°2-3, pp. 151 - 169, 2002.

[4] A. Le Pape and J. Lecanu, ‘3D Navier-Stokes Computations of a Stall-Regulated Wind Turbine, The Science of Making Torque from Wind’, ONERA, pp. 78 - 88, 2004; Wind Energy Journal, Vol. 7, N°4, pp. 309 - 324, 2004.

[5] C. Thumthae and T. Chitsomboon, ‘Optimal Pitch for Untwisted Blade Horizontal Axis Wind Turbine’, The 2nd Joint International Conference on ‘Sustainable Energy and Environment (SEE 2006)’, B-034, November 21-23, Bangkok, Thailand, 2006.