# Simulation instationnaire de l'écoulement autour d'un rotor éolien à axe horizontal

Mohamed Maïzi¹*, Rabah Dizène² et Ouahiba Guerri¹

¹ Division Energie Eolienne

Centre de Développement des Energies Renouvelables, CDER
B.P. 62, Route de l'Observatoire, Bouzaréah, Alger, Algérie

² Faculté de Génie des Procédés et de Génie Mécanique

Université des Sciences et de la Technologie Houari Boumediène
B.P. 32, El Alia, Bab Ezzouar, Alger, Algérie

(reçu le 29 Septembre 2012 – accepté le 29 Décembre 2012)

**Résumé** - Les pales d'une éolienne sont soumises à des chargements aérodynamiques, génèrent des écoulements extrêmement complexes. Nous nous intéressons à la dynamique des écoulements autour des pales d'éolienne à axe horizontal en rotation. Ce papier présente une méthode de simulation numérique pour la prédiction de l'écoulement autour d'une éolienne à axe horizontal opérant en régime instationnaire, tridimensionnel turbulent utilisant l'approche basée sur les équations de Navier-Stokes moyennées (URANS) en utilisant le code 'CFX'. Un maillage structuré a été adopté en utilisant la technique des interfaces (interface frozen rotor) entre le corps mobile (la pale) et le corps fixe (l'entourage). Le modèle de turbulence SST a été utilisé pour décrire l'écoulement turbulent. Les résultats de simulation ont été comparés avec les résultats expérimentaux, pour les conditions de rotation. En général, des bonnes concordances ont été notées. Une attention spéciale est portée sur les effets instationnaires, conséquents de rotation du rotor qui est mise en évidence par le calcul des pressions fluctuantes sur l'interface entre les parties fixes et mobiles en aval et en amont de la machine. Une analyse fréquentielle de ces signaux en utilisant une transformée de Fourier discrète (FFT) afin d'analyser le comportement instationnaire de la pale.

**Abstract** - In this work the analysis of the flow dynamics around horizontal axis wind turbine blades was carried out using the software CFX with a turbulent three-dimensional flows approach based on the averages Navier-Stokes equations (URANS). A structured mesh was adopted using interfaces technique (frozen rotor interface) between moving parts (blade) and fixed parts, SST turbulence model is used to describe the turbulent flow. The solver ability to modeling the flow field around the wind turbine blades was validated by the experimental results of NREL phase II wind turbine blades. The results show that numerical data have a good agreement with experimental results. A special attention is related to the instationary effects, the rotor rotation consequence which is evidenced by calculating the fluctuating pressure on the interface between the fixed and moving parts of the machine, an analysis of these frequency signals using a discrete Fourier transform (FFT) to analyze the unsteady behavior of the blade.

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

## 1. INTRODUCTION

Parmi les sources potentielles d'énergies renouvelables, l'éolien représente une part importante. L'énergie éolienne, si elle ne peut pas compenser ou remplacer les sources

<footnote>* mohamed.maizi@cder.dz</footnote>
classiques d'énergies fossiles, constitue cependant une voie intéressante de réduction des émissions de CO₂. Le développement de l'énergie éolienne s'est très nettement accéléré depuis une dizaine d'années. Par voie de conséquence, le coût du kilowattheure éolien a progressivement baissé et la fiabilité des machines s'est améliorée.

Les développements réalisés depuis une dizaine d'années portent principalement sur l'accroissement des échelles et des productions unitaires des éoliennes à axe horizontal. Les plus grandes machines atteignent actuellement des puissances unitaires de production de 5 MW.

Avant leur implantation sur site, les éoliennes sont dimensionnées pour un point de fonctionnement nominal, c'est-à-dire en fonction de la puissance attendue pour une vitesse de vent donnée. Dans la plupart des cas, le choix s'est porté sur les éoliennes à haut rendement, c'est-à-dire sur les éoliennes rapides à axe horizontal bi- ou tripales [1].

La prédiction d'écoulement autour d'une éolienne à axe horizontal est un élément clé dans la conception aérodynamique. En règle générale, dans l'ingénierie, la conception d'un rotor est basée sur des approches simplifiées, le concept du disque actuateur et la théorie d'élément de pales. Ces modèles s'avèrent être limités à des écoulements stationnaires, alors que l'écoulement autour du rotor éolien est un écoulement instationnaire en 3-D.

Par suite de l'évolution rapide de la puissance de calcul dans ces deux dernières décennies, la signification des simulations numériques pour la conception des machines tournantes et dans les éoliennes en particulier, s'est développée jusqu'à un degré considérable [2].

<figure><img src="image_1.png" /><figcaption>Fig. 1: Géométrie de la pale éolienne phase II en 3D</figcaption></figure>

Les modèles numériques de dynamique des fluides sont aujourd'hui très performants. La vitesse élevée des processeurs permet de réduire les temps de calcul pour obtenir des solutions numériques complexes et précises. Les simulations numériques se comparent maintenant à de véritables essais expérimentaux virtuels.

Les modèles physiques utilisés peuvent être en 2D ou 3D et de complexités différentes. Les modèles 3D les plus avancés comme ceux utilisés dans Ansys CFX utilisent la théorie des éléments finis ou des volumes finis pour résoudre les équations de Navier-Stokes.

L'industrie éolienne est actuellement confrontée à de nombreuses difficultés dans la construction des machines éoliennes causées par l'incapacité de prédire la charge structurale et de la puissance de sortie.
Le but ultime de ce travail est de développer une méthode CFD permettant d'analyser le comportement instationnaire de l'écoulement turbulent autour des pales d'une éolienne à axe horizontal et ce, afin de prédire convenablement les effets tridimensionnels et instationnaire, de la turbulence et par conséquent d'améliorer la précision des prédictions numériques. Le champ d'écoulement autour des pales sera décrit par les équations de Navier-Stokes moyennées (URANS) [3].

Dans le présent travail, des simulations instationnaires autour d'un rotor éolien type NREL phase II seront présentées. L'éolienne étudiée est une éolienne tripale de 10 m de diamètre ayant un profil S809 (Fig. 1). Ce choix a été fait en raison de la présence de données expérimentales sur cette éolienne qui constitue actuellement une référence d'étude.

Un modèle de turbulence $\kappa - \omega$ SST 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 Ansys CFX. 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 qui régissent 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 \times \left( \frac{\partial u_i}{\partial t} + u_j \times \frac{\partial u_i}{\partial x_j} \right) = - \frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} (\tau_{ij}) \qquad (2) $$

$$ \tau_{ij} = \mu \times \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \qquad (3) $$

Avec $u$, la vitesse du fluide, $\rho$, la masse volumique et $\mu$, le coefficient de viscosité dynamique

### 2.1 Modèle de Turbulence

Le modèle $\kappa - \omega$ SST est une évolution du modèle $\kappa - \omega$ BSL (Baseline) qui est lui-même une combinaison du modèle $\kappa - \epsilon$ et du modèle $\kappa - \omega$ de Wilcox. Menter [4] a combiné les deux modèles $\kappa - \epsilon$ et $\kappa - \omega$ d'une manière à utiliser le deuxième modèle près de la paroi, en passant à travers une fonction du modèle $\kappa - \epsilon$, loin de la paroi, sachant que ces équations ont été transformées à un format du modèle $\kappa - \omega$.

Cette méthode intègre la suite 'suggestion Bradshaw' que la contrainte de Reynolds de cisaillement doit être considérée comme étant proportionnelle à l'énergie cinétique turbulente. Le modèle $\kappa - \omega$ SST représente le transport des contraintes de cisaillement
turbulent et donne des prédictions très précises sur le commencement et l'évolution des séparations du flux sous les contraintes du gradient de pression. Les équations finales sont les suivantes [5]:

$$ \frac{\partial(\rho \times \kappa)}{\partial t} + \overline{U_j} \times \frac{\partial(\rho \times \kappa)}{\partial x_j} = \tau_{ij} \times \frac{\partial \overline{U}}{\partial x_j} - \beta^* \times \rho \times \omega \times \kappa + \frac{\partial}{\partial x_j} \left[ (\mu + \sigma_{\kappa} \times \mu_1) \times \frac{\partial \kappa}{\partial x_j} \right] \quad (4) $$

$\kappa$ représente l'énergie cinétique du mouvement turbulent.

$$ \frac{\partial(\rho \times \omega)}{\partial t} + \overline{U_j} \times \frac{\partial(\rho \times \omega)}{\partial x_j} = \frac{\mu}{\mu_t} \tau_{ij} \times \frac{\partial \overline{U}}{\partial x_j} - \beta \times \rho \times \omega^2 + \frac{\partial}{\partial x_j} \left[ (\mu + \sigma_{\omega} \times \mu_1) \times \frac{\partial \omega}{\partial x_j} \right] + 2(1 - F_1) \times \rho \times \sigma_{\omega}^2 \times \frac{1}{\omega} \times \frac{\partial \kappa}{\partial x_i} \times \frac{\partial \omega}{\partial x_j} \quad (5) $$

Les constantes dans le modèle $\kappa - \omega$ SST sont calculées comme suit: Si $\varphi$ est une constante dans le modèle $\kappa - \omega$ SST, $\varphi_1$ et $\varphi_2$ sont des constantes dans les modèles $\kappa - \omega$ et $\kappa - \varepsilon$, respectivement [5], alors

$$ \varphi = F_1 \times \varphi_1 + (1 - F_1) \times \varphi_2 \quad (6) $$

Les constantes pour $\kappa - \omega$ sont: $\sigma_{\kappa 1} = 0.85$; $\sigma_{\omega 1} = 0.5$; $\beta_1 = 0.075$; $\beta^* = 0.09$; $a_1 = 0.31$; $\kappa = 0.41$.

$$ \gamma_1 = \frac{\beta_1}{\beta^*} - \frac{\sigma_{\omega_1 \kappa^2}}{\sqrt{\beta^*}} \quad (7) $$

Les constantes pour $\kappa - \varepsilon$ sont: $\sigma_{\kappa 2} = 1.0$; $\sigma_{\omega 2} = 0.856$; $\beta_2 = 0.0828$; $\beta^* = 0.09$; $\beta_{\kappa 1} = 0.075$.

$$ \gamma_1 = \frac{\beta_2}{\beta^*} - \frac{\sigma_{\omega_2 \kappa^2}}{\sqrt{\beta^*}} \quad (8) $$

$\tau_{ij}$ est donnée par l'expression suivante:

$$ \tau_{ij} - \rho \times \overline{U_i U_j^*} = -\frac{2}{3} \times \rho \times \kappa \times \delta_{ij} + \mu_i \times \left( \frac{\overline{\partial U_i}}{\partial x_j} + \frac{\overline{\partial U_j}}{\partial x_i} \right) \quad (9) $$

### 3. DOMAINE DE CALCUL ET MAILLAGE

En négligeant les effets de la tour, de la nacelle et du sol, nous avons considéré une configuration du domaine de calcul. Le domaine de calcul pour la pale d'éolienne est illustré sur la figure 3. Une grille à deux blocs est utilisée. Puisque l'éolienne dispose de trois pales, le domaine est choisi sous forme d'un tiers de cylindre pour appliquer les conditions de périodicité.

Le domaine de calcul se compose d'un sous-domaine intérieur DI (mobile) qui englobe la pale du rotor (Fig. 2) où le maillage est plus raffiné et un grand domaine extérieur DE sous forme d'un tiers de cylindre (fixe) (Fig. 3).
Il est utilisé pour la simulation de la pale en rotation. Le maillage généré se compose de $3 \times 10^6$ nœuds (Fig. 4).

<figure><img src="image_2.png" /><figcaption>Fig. 2: Sous domaine (mobile)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 3: Domaine de calcul</figcaption></figure>

Le maillage utilisé a été généré avec le préprocesseur ICEM. Les frontières extérieures du domaine DE ont été fixées respectivement à 3 et 10 fois la longueur de la pale du centre du domaine. L'interface DE-DI permet de mettre en œuvre la technique du maillage glissant disponible dans le logiciel (interface frozen rotor), l'information passe du domaine fixe au domaine mobile via une interpolation.

### Condition aux limites

<table><thead><tr><th>En amont</th><th>En aval</th></tr></thead><tbody><tr><td>U<sub>x</sub> = U<sub>∞</sub></td><td>P = P<sub>∞</sub> = const</td></tr><tr><td>U<sub>y</sub> = 0</td><td>Vitesse de rotation de la pale- ω<sub>Z</sub> = 71.68 tours/min</td></tr><tr><td>U<sub>Z</sub> = 0</td><td>Avec- U<sub>∞</sub> : Vitesse de vent en amont</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 4: Domaine de calcul avec maillage</figcaption></figure>

## 4. RESULTATS

Le code commercial ANSYS CFX est utilisé pour la simulation présentée. Les figures 5 et 6 présentent le coefficient de pression le long du profil pour une vitesse 7.2 m/s. Tout d'abord, on remarque que les distributions de pression très similaires avec l'expérimental pour l'intrados et l'extrados surtout pour r/R = 0.80.

Les distributions de pression proposées sont presque superposées aux valeurs expérimentales. Les valeurs du coefficient de pression proposées par la simulation pour
l'intrados sont légèrement supérieures aux valeurs expérimentales jusqu'à y/c = 0.45 pour r/R =0.63, puis suivent parfaitement les valeurs expérimentales.

Ensuite on suppose qu'une prédiction trop importante, malgré la sous-estimation du coefficient de pression avant y/c = 0.45 pour r/R = 0.63.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Distribution de coefficient de pression Cp - Uₓ=0.72 m/s (71.68 tour/min, r/R = 0.80)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: Distribution de coefficient de pression Cp - Uₓ=0.72 m/s (71.68 tour/min, r/R=0.63)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7: Vue en 3D des lignes de courant</figcaption></figure>

Il est clair que près de la racine (en pied de la pale), une importante séparation de flux sur intrados, un décollement apparaît, lorsque l'angle incidence est important, c'est-à-dire que la couche limite décolle sans rattachement de la couche limite sur le profil

Cette séparation influence largement sur la structure de l'écoulement derrière le rotor. Les gradients de pression adverses sont forts et ne sont plus complètement balancés par l'accélération convective (et les forces de viscosité).

On remarque que la limitation des viscosités tourbillonnaires, permet au modèle $κ - ω$ SST de modéliser le transport des contraintes de cisaillement turbulent, et donc lui permet de modéliser le déplacement des turbulences. Ce qui lui permet notamment de réaliser de bonnes prédictions pour les structures tourbillonnaires.
# 5. SIMULATION INSTATIONNAIRE

Les calculs stationnaires basés sur la technique du maillage glissant débutent après convergence de la solution stationnaire utilisant l'approche 'Frozen-Rotor'. Ceci étant après la stabilisation des variations des grandeurs de l'écoulement issues du calcul URANS.

Afin de localiser les instationnarités et les zones susceptibles d'être le siège des perturbations aérodynamiques dans l'écoulement autour des pales du rotor, les pressions statiques fluctuantes ont été calculées sur différents points au niveau de l'interface en amont et en aval de la pale (Fig. 13), où la variation de pression statique en fonction du temps correspond à 3 tours du rotor.

<figure><img src="image_8.png" /><figcaption>Fig. 8: Positions des points de calcul</figcaption></figure>

Le signal récupéré est un signal de pression statique mettant en évidence les phénomènes physiques. L'avantage de cette grandeur est d'être également invariant par changement de repère (passage du repère tournant au repère fixé). L'analyse spectrale de ces signaux pour ces points a fourni des signaux fréquentiels correspondants en utilisant une transformée de Fourier discrète (FFT).

Afin d'analyser le comportement instationnaire de l'écoulement autour de la pale, la simulation numérique instationnaire a été réalisée pour un cas d'un fonctionnement nominal.

<figure><img src="image_9.png" /><figcaption>Fig. 9: Fluctuations de pression en fonction du temps</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10: Fluctuations des signaux fréquentiels</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 11: Fluctuations de pression en fonction du temps</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 12: Fluctuations des signaux fréquentiels</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 13: Fluctuations de pression en fonction du temps</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 14: Fluctuations des signaux fréquentiels</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 15: Fluctuations de pression en fonction du temps</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 16: Fluctuations des signaux fréquentiels</figcaption></figure>

Les figures 9, 11, 13 et 15 montrent que le signal de pression issue des simulations est pratiquement périodique. On constate clairement qu'ils font apparaître de manière dominante 08 évolutions correspondantes au passage des pales du rotor (3 pales).
La figure 9 montre que la pression du point 1 est mise en oscillation périodique relativement grande, de l'ordre de 15 Pa par rapport aux autres points. La figure 10 montre la présence d'un pic dominant sur le signal fréquentiel (FPP). Ce pic correspond à la fréquence fondamentale de passage de la pale et autres harmoniques qui correspondent aux multiples de la fréquence fondamentale.

L'origine des fluctuations de forces périodiques sur les pales est due au fait que la vitesse moyenne dans le repère fixe (vitesse absolue) à l'entrée du rotor n'est pas parfaitement uniforme, en particulier sur le radial à cause de l'absence de vrillage de la pale. Ce qui entraîne une variation périodique de la résultante des forces aérodynamiques sur la pale à des multiples de la fréquence de rotation.

D'après l'analyse fréquentielle des signaux de ces différents points, les figures montrent que l'amplitude de l'harmonique fondamentale (la fréquence de passage de pale FPP 3.5 Hz) est plus importante pour le point 1 situé à 0.1 de la hauteur de la pale en aval du rotor, puis l'amplitude de fréquence pour le point 3 situé à 0.1 de la hauteur de la pale en amont du rotor et l'amplitude la plus faible correspond au point 4 situé à 1.1 de la hauteur de la pale en amont du rotor.

On constate que la fréquence de passage des pales est toujours la structure dominante du point de vue temporel et fréquentiel pour les points 1 et 3. La zone la plus perturbée est la zone qui correspond à la plus importante amplitude, ou bien qui correspond à la fluctuation de pression la plus importante, c'est la zone du point 2 où apparaît l'influence de séparation de flux sur intrados (Fig. 7), où le facteur dominant est le plus perturbant derrière la pale du rotor.

La zone en aval de la pale subit alors de fortes fluctuations de pression résultantes de l'interaction du sillage avec l'écoulement moyen, ce qui se traduit par une forte émission de fréquence au passage des pales.

Le mouvement relatif des pales du rotor constitue la principale source d'instationnarité. On constate que la fréquence de passage des pales est toujours la structure dominante du point de vue temporel et fréquentielle.

L'analyse fréquentielle des fluctuations de pression a montré l'émergence d'un pic à la fréquence de passage des pales (3.5 Hz). Pour les points 2 et 4 qui sont proches à l'extrémité de la pale où l'amplitude de FFP est faible par rapport à 2FPP et 3FPP.

On peut voir clairement que l'amplitude de fréquence 10.5 Hz est le facteur dominant pour les deux points 2 et 4 où l'amplitude du point 2 est plus importante que l'amplitude du point 4 de la même fréquence, puisque le point 2 est situé dans la zone qui subit une interaction entre les tourbillons marginaux à l'extrémité de la pale et l'écoulement moyen.

## 6. CONCLUSION

Dans ce papier, nous avons présenté un travail qui consiste en l'étude de l'écoulement autour d'un rotor éolien en rotation. Pour le 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 II. Une comparaison des répartitions du coefficient de pression entre la simulation et l'expérimental a été présentée. Nous avons une bonne prédiction des distributions de pression.
La modélisation numérique basée sur la résolution des équations URANS a été réalisée. L'analyse des champs de pression a montré qu'il existe de fortes fluctuations de pression notamment en aval du rotor, caractérisées par la présence d'un fort gradient de grandeurs locales de l'écoulement.

L'analyse fréquentielle des fluctuations de pression a montré l'émergence d'un pic à la fréquence de passage des pales. L'augmentation de la fluctuation de pression qui se produit à cause du mouvement des pales est assimilée à des forces de pression qui se propagent vers l'amont et vers l'aval du rotor sous forme d'onde de pression.

Le mouvement relatif des pales du rotor constitue la principale source d'instationnarité, et on constate que la fréquence de passage des pales est toujours la structure dominante du point de vue temporel et fréquentielle.

## REFERENCES

[1] J.L. Menet et N. Bourabaa, 'Prévision des Performances Aérodynamiques d'un Nouveau Type d'Eolienne à Axe Vertical: Le Rotor à Contrevoiles', 18ème Congrès Français de Mécanique, Grenoble, 27-31 Août 2007.

[2] M. Maïzi, 'Interaction Stator-Rotor: Cas d'une Soufflante Axiale', Premières Journées sur la Physique, la Sûreté et la Technologie des Réacteurs Nucléaires, Centre de Recherche Nucléaire de Draria, 2008.

[3] J.L. Harion, 'Simulations Numériques des Performances d'une Nouvelle Géométrie d'Eolienne à Axe Vertical', 18ème Congrès Français de Mécanique, Grenoble, 27-31 Août 2007.

[4] F.R. Menter, 'Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications', American Institute of Aeronautics and Astronautics Journal, 'AIAA Journal', Vol. 32, N°8, pp. 1598 - 1605, 1994.

[5] R. Dizène and I. Yous, 'Numerical Investigations on Rotor Dynamic Stall of Horizontal Axis Wind Turbine', Revue des Energies Renouvelables, CICME'08, Sousse, pp. 119 – 126, 2008.