# Modélisation d'un écoulement pulsatif et oscillatoire autour d'une pale d'éolienne HAWT
## Application au décrochage dynamique

Ali Nouioua et Rabah Dizene

Laboratoire de Mécanique Avancée, LMA
Université des Sciences et de la Technologie Houari Boumediene, USTHB
B.P. 32, El Alia, Bab Ezzouar, Algiers, Algeria

(reçu le 10 Août 2014 – accepté le 25 Mars 2015)

**Résumé** – L'étude d'une turbine éolienne à axe horizontal (HAWT) soumise à des conditions de fonctionnement en présence de vent est très complexe, car la machine est soumise à une variation instantanée de la vitesse et de la direction du vent, qui peut conduire inévitablement les aubes de turbine en phase de décrochage dynamique. L'intensité de ce phénomène conduit à une baisse soudaine du coefficient de portance et donc à la dégradation de ses performances aérodynamiques. Ce présent travail est une simulation numérique bidimensionnelle de la variation instantanée de la direction du vent par rapport à la pale de turbine, comparée aux données expérimentales disponibles dans la littérature. L'effet de la variation de l'angle d'incidence d'un profil S809 est donc étudié en utilisant une approche URANS pour prédire le décrochage. La modélisation de la turbulence est assurée à l'aide des modèles k-ε et k-ω SST. Le logiciel Fluent est utilisé pour la résolution numérique des équations URANS. Les résultats numériques sont comparés aux données expérimentales rapportées par la littérature. Les résultats montrent l'importance de prendre en compte le comportement de l'écoulement instationnaire dans les analyses afin d'obtenir une estimation précise des charges aérodynamiques qui agissent sur l'éolienne.

**Abstract** – The study of a horizontal axis wind turbine (HAWT) submitted to operation conditions with wind is very complex, because the machine is subject to instant variation of the speed and the wind direction, which leads necessarily the turbine blades to dynamic pull-out phase. The intensity of this phenomenon generates to a brutal decrease in lift coefficient and so the deterioration of its aerodynamic performances. The present work is a numerical bidimensional simulation of the instant variation of the wind direction according to the turbine blade, compared to available experimental datas in literature. The variation effect od the angle of incidence of an S809 profil is then studied using an URANS approach to predict the off-hook. Model the disturbance is assured with k-ε and k-ω SST models. The software Fluent is used to resolve numerically urans equations. Numerical results are compared with experimental datas reported by the literature. Results show the importance of taking in account the behaviour of the unsteady flow in analyzes in order to obtain precize estimation of aerodynamic charges which affect the wind turbine.

**Mots clés:** Eolienne HAWT - Décrochage dynamique - Profil aérodynamique - CFD - Fluent - Méthodes RANS - URANS.

# 1. INTRODUCTION

L'étude aérodynamique est une science qui consiste à étudier le comportement d'écoulement des gaz (généralement l'air) autour d'un obstacle. Autrement dit, à étudier l'interaction entre le milieu solide et le milieu fluide en mouvement. Malheureusement, l'étude expérimentale est couteuse et elle n'est pas évidente dans la plupart des cas.

Le CFD 'Computational Fluid Dynamics' est devenue une partie intégrale de la conception de l'ingénierie et un environnement d'analyse pour beaucoup de compagnies
qui ont besoin de l'aptitude de prévoir la performance de nouvelles conceptions ou les processus, avant d'être manufacturé ou mises en application.

Les solutions de CFD à partir des ANSYS sont basées sur la technologie prouvée du logiciel Fluent. Les forces aérodynamiques de la pale responsable de la production de l'énergie doivent être augmentées afin de maximiser la capture de l'énergie. Dans un même temps, les effets de ces charges aérodynamiques sur les composants de l'éolienne doivent être atténués afin de prolonger la durée de vie de la machine.

Ce travail a pour objectif de prédire le phénomène du décrochage dynamique d'un écoulement turbulent autour d'une pale de rotor HAWT. Le profil de la pale choisi pour ces simulations est le NREL (National Renewable energy Laboratory) S809. Nous avons choisi ce profil car il existe un nombre important d'expérimentation dans la littérature le concernant.

La variation du champ de vitesse est décrite par un programme C++ qu'on introduit directement dans le solveur. La modélisation de la turbulence est assurée par l'utilisation de trois modèles de turbulence qui sont le k-ε Standard, k-ε RNG et le k-ω SST.

Dans le but de valider cette démarche numérique, une comparaison des résultats numériques avec des données expérimentales [1, 2], issue de la bibliographie a été faite. L'écoulement considéré est en 2D autour du profil animé d'un mouvement d'oscillation, simulé par un champ de vitesse à incidence variable.

## 2. MODELE MATHEMATIQUE

Pour un écoulement incompressible, les équations régissant la conservation de la masse et la quantité de mouvement sont données respectivement par:

### Equation de la continuité

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

### Equation de quantité de mouvement

$$ \frac{\partial U_i}{\partial t} + U_j \frac{\partial U_i}{\partial x_j} = - \frac{1}{\rho} \frac{\partial P}{\partial x_j} + \nu \frac{\partial^2 U_i}{\partial x_j \partial x_j} \qquad (2) $$

### 2.1 Modèle de turbulence

Les modèles qui se fondent sur l'approche statistique sont les plus répandus. Ils reposent sur la décomposition classique de Reynolds en quantités moyenne et fluctuante.

Parmi les nombreux modèles de turbulence, le code de calcul Fluent propose trois méthodes de fermeture basées sur l'approche statistique:

* Le modèle k-ε et ses variantes;

* Le modèle de turbulence k-ω;

* Le modèle des contraintes de Reynolds (RSM) et ses variantes.

On s'intéresse à ces trois modèles de turbulence, puisque Fluent propose d'autres modèles, tel que Spalart-Allmaras et aussi l'approche LES.

### 2.2 Modèle k-ε [3]

Ce modèle propose deux équations pour déterminer la viscosité turbulente, ce dernier fait donc partie des modèles dits à deux équations. Ce modèle est le plus connu
et le plus utilisé dans l'industrie jusqu'à présent. Il se base sur les équations de l'énergie cinétique $k$ et de la dissipation turbulente $\varepsilon$.

Ce modèle suppose que le régime de turbulence est pleinement établi dans tout le domaine et que les effets de la viscosité moléculaire sont négligeables par rapport à ceux de la viscosité turbulente (loin des parois). Il est basé sur l'hypothèse de Boussinesq.

Dans ce travail, on a opté pour les variantes de $k-\varepsilon$ suivantes:

* $k-\varepsilon$ standard.

* $k-\varepsilon$ RNG (Reynolds Normalization Group).

## 2.3 Modèle $k-\omega$ SST [4]

Le modèle SST ou transport des contraintes de cisaillement turbulentes été développé efficacement par Menter (1996), il est dérivé du modèle $k-\omega$ Standard. Ce modèle combine la robustesse et la précision de la formulation du modèle $k-\omega$ dans la région proche de la paroi avec le modèle $k-\varepsilon$ et toutes ses caractéristiques citées plus haut pour les écoulements libres loin de la paroi.

Ce modèle est principalement conseillé dans le cas d'applications, telles que les fluides subissant des changements soudains de contrainte, s'écoulant dans des surfaces courbées ou dans les cas de séparation de la couche limite.

Nous avons choisi d'adopter ce modèle à notre phénomène en vue des différents avantages qu'il offre pour ce type de configurations.

# 3. MODELISATION NUMERIQUE

L'étude de notre problème implique la résolution d'un système d'équations mathématiques, aux dérivées partielles de type elliptique. Cela nécessite une approche par une voie numérique. Nous avons opté dans cette étude pour l'utilisation du code de calcul 'Fluent' comme solveur, et le 'Gambit' pour la description géométrique et le maillage.

'Gambit' est un logiciel de génération tridimensionnelle de maillage. Les maillages générés peuvent être utilisés par de nombreux logiciels d'analyse que ce soit pour la résistance des matériaux par éléments finis ou pour la simulation des écoulements fluides dont fait parti 'Fluent'.

## 3.1 Méthode des volumes finis

Plusieurs méthodes de discrétisation des équations différentielles aux dérivées partielles sont utilisées actuellement, telles que: la méthode des volumes finis, des différences finies et des éléments finis.

Le code Fluent utilise la méthode des volumes finis, car elle apporte une généralité beaucoup plus importante par son formalisme, son traitement conservatif, ainsi que par son adaptation au problème physique, et par sa simplicité dans la linéarisation des termes. Ces critères ont donné à cette méthode, une stabilité numérique et une convergence plus efficace.

La méthode des volumes finis est caractérisée par son avantage à satisfaire la conservation de masse, de quantité de mouvement et d'énergie dans tous les volumes finis ainsi dans tout le domaine de calcul. Elle facilite la linéarisation des termes non linéaires dans les équations de conservation, telle que le terme source par exemple.
La méthode consiste à partager le domaine de calcul en plusieurs volumes, où chaque volume entoure un nœud. En utilisant différents schémas d'approximations on peut intégrer les termes des équations différentielles modélisant sur chaque volume de contrôle, où les valeurs et les quantités sont stockées aux nœuds du volume de contrôle.

Ces équations algébriques produites expriment la conservation des quantités pour le volume de contrôle et pour tout le domaine de calcul.

## 3.2 Domaine et Maillage

Les dimensions du domaine de calculs sont les suivants: Une hauteur de 25 cordes, une longueur de 32.5 corde dont 12.5 en amont du bord de fuite et 20 en aval.

Le nombre de Reynolds amont utilisé pour ces simulations est de $10^6$, le maillage est structuré, et ce maillage est présenté sur la figure 1, possédant 12451 nœuds.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Maillage du domaine d'étude</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Fig. 2: Maillage autour du profil</figcaption></figure>

## 3.3 Equation de champ de vitesse

Pour cette étape très importante dans notre étude, on utilisera l'UDF (User Defined Function) qui est une composante du logiciel Fluent, à l'aide d'un programme 'C++'. Pour qu'on puisse varier l'angle d'incidence ou le champ de vitesse en fonction du temps, l'équation qui gouverne la variation du champ de vitesse est décrite comme suit:

$$ \alpha = \alpha_0 + \alpha_m \times \cos(\Omega t) \quad (3) $$

Avec, $\alpha_0$, l'angle d'incidence moyen, $\alpha_m$, l'amplitude maximale des oscillations et $\Omega$, la fréquence d'oscillation.
<figure><img src="image_3.png" /><figcaption>Fig. 3: Variation du champ de vitesse</figcaption></figure>

# 4. RESULTATS

Les résultats de la simulation numérique présentés sont obtenus à l'aide de la résolution des équations de Navier-Stockes pour un écoulement 2D incompressible, avec un nombre de Reynolds $Re = 10^6$.

Une discussion qualitative est faite grâce à la présentation des valeurs isovitesse instationnaire respectivement avec le modèle $k-\epsilon$ standard, $k-\epsilon$ RNG et le $k-\omega$ SST.

La confirmation de la description précédente est réalisée grâce à la présentation de la distribution du coefficient de pression dans les mêmes conditions et respectivement pour les trois modèles de turbulence cités juste avant.

On termine enfin par montrer un résultat quantitatif important qui concerne le coefficient de portance obtenu avec les trois modèles et compare avec un résultat expérimental.

## 4.1 Distribution de vitesse

La figure 4 montre les isovitesses de l'écoulement autour du profil pour les trois modèles, pour $\alpha = 15^\circ$.

On observe à $\alpha = 15^\circ$ déjà un écoulement qui commence à décoller, à $\alpha = 20^\circ$, l'écoulement est complètement décollé sur l'extrados.

Par ailleurs, on observe que le décrochage tourbillonnaire prend naissance juste derrière le bord d'attaque comme le montrent les modèles $k-\epsilon$ RNG et le $k-\omega$ SST.

<figure><img src="image_4.png" /><figcaption>Fig. 4: Isovaleurs de vitesse pour $\alpha = 15^\circ$</figcaption></figure>

## 4.2 Distribution de la pression statique

La figure 5 montre les distributions de la pression statique en fonction de la position $x/c$ respectivement pour les trois modèles pour un angle d'attaque de $\alpha = 15^\circ$, la
représentation de la pression est quasi semblable entre les trois modèles sur la surface entière du profil. Les distributions sur les faces intrados et extrados présentent la même tendance que celle observée généralement sur les résultats expérimentaux. On respecte l’existence d’un bull de décollement qui se situe aux environs de (x/c = 0.45), ce qui est conforme aux résultats trouvés dans la littérature (Reuss Ramsay, M 1995) [5].

Pour les différents angles notables sur ces distributions entre les trois modèles, on note cependant des pressions fortement négatives à $\alpha =15°$ et qui s’accentuent à $\alpha =20°$ sur la face extrados du profil. Ce pic de dépression localisé juste après le bord d’attaque, diminue au fur et à mesure que l’écoulement avance vers le bord de fuite.

Certains auteurs associent l’ampleur de ce pic de dépression à sa dépendance de la pression dynamique le long de l’envergure. Ces gradients de pression négatifs d’envergure et de corde fournissent un mécanisme par lequel l’on va se diriger vers des zones fortement décollée. Cette analyse et ces constatations laissent prédire que le décrochage dynamique risque d’avoir lieu aux alentours de $\alpha = 20°$.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Coefficients de pression en fonction de la position x/c pour $\alpha =15°$</figcaption></figure>

La comparaison du phénomène entre les trois modèles de turbulence fait ressortir que le modèle k-ε standard se distingue des deux autres, dont la valeur $\alpha = 25°$ où les distributions de pressions obtenues restent relativement constantes sur la surface de l’extrados et vers le bord d’attaque.

## 4.3 Coefficient de portance

La prédiction de l’écoulement dans toute sa complexité est importante non seulement pour déterminer la puissance moyenne récupérée, mais aussi pour anticiper sur les efforts cycliques qui induisent la fatigue des matériaux. En effet, ces sollicitations affectent aussi la structure qui assure la fixation des machines.

Le coefficient de portance $C_1$ est défini par l’équation suivante:

$$L = 1/2 \times C_1 \times \rho \times V^2 \times S \quad (4)$$

Avec, L, la portance, (N); $C_1$, le coefficient de portance, (-); $\rho$, la masse volumique de l’air, (kg/m$^3$); V, la vitesse, (m/s); S, la surface à l’air, (m$^2$).
<figure><img src="image_6.png" /><figcaption>Fig. 6: Coefficients de portance pour les trois modèles de turbulence pour Re=10⁶</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7: Résultats expérimentaux de coefficient de portance en fonction de l'angle d'attaque selon I. Dobrev [1] (à gauche), et P. Walter [2] (à droite)</figcaption></figure>

La figure 6 présente l'évolution du coefficient de portance en fonction de l'angle de l'incidence $\alpha$ pendant un intervalle de temps $t = 2.75$ secondes et pour un Reynolds de $10^6$.

Comparés aux modèles RSM et $k-\epsilon$, ainsi qu'à la courbe expérimentale de la référence [1, 2], nos calculs reproduisent de manière très satisfaisantes le comportement du coefficient $C_l$ autour du profil S809.

Les valeurs maximales de $C_l$ sont bien prédites dans la plage des angles d'attaque de $\alpha=15°$ à $\alpha=22°$, alors que la courbe expérimentale montre un coefficient de portance $C_l$ maximale se trouvant entre $\alpha=22°$ à $\alpha=18°$. Cependant, les valeurs de ce présent calcul sont au dessus des valeurs calculées et mesurées dans les références [1, 2].

Les résultats des modèles $k-\epsilon$ standard et $k-\epsilon$ RNG sont moins précis que ceux du modèle $k-\omega$ SST de [1], comparant aux valeurs expérimentales. Des insuffisances numériques essentiellement dues au maillage et au type de maillage utilisé sont à suspectés. Cependant, la prospection des résultats expérimentaux de [1, 2] montre que les mesures effectuées ne prenaient pas en compte la variation instationnaire de l'incidence simultanément avec l'instationnarité, telle que considéré dans cette étude. Il se pourrait donc que ces écarts proviennent de ces conditions de mesures.
# 5. CONCLUSION

Une modélisation d'un champ de vitesse variable est proposée dans cette étude afin de simuler les variations de la vitesse du vent auquel est soumis un rotor d'éolienne d'axe horizontal (HAWT) de profil NREL S809. Les résultats de calcul instationnaire et turbulent sont basés sur la discussion des champs de vitesse, de pression statique, ainsi que l'évolution du coefficient de portance.

Les résultats de calcul de l'évolution globale du champ d'écoulement autour du profil montrent un comportement identique dans son ensemble aux résultats trouvés dans la littérature [6], qui montrent que le déclenchement du décrochage se fait à partir du bord d'attaque et à $\alpha =15°$.

Les résultats des coefficients de pression statique reproduisent de façon satisfaisante la zone de transition suspectée être localisée à une distance de 40% de corde à partir du bord d'attaque, qui est aussi le siège du début du décollement de la couche. Cependant les courbes expérimentales montrent des perturbations de la couche limite, dès que la position $x=20\%$ de corde est atteinte.

Les coefficients de portance $C_l$ de notre calcul présentent une prédiction jugée très satisfaisante, comparés à ceux de l'expérience [2], puisque le décrochage dynamique est suspecté être aux environ de $\alpha =20°$.

Comparés avec l'expérience, les modèles $k-\omega$ SST et RSM semblent les mieux adaptés pour la simulation du phénomène du décrochage.

Un travail plus approfondi sur les effets de variation de l'angle d'attaque est à entreprendre comme perspective à donner à ce présent travail, car il existe plusieurs mécanismes à maîtriser à l'aide d'autres modélisations quelles que soient empiriques ou périodiques.

# REFERENCES

[1] I. Dobrev, 'Modèle Hybride de Surface Active pour l'Analyse du Comportement Aérodynamique des Rotors Eoliens à Pales Rigides ou Déformables', Thèse de Doctorat, Ecole Nationale Supérieure d'Arts et Métiers, Paris, 2009.

[2] W.P. Wolfe and S.S. Ochs, 'CFD Calculations of S809 Aerodynamic Characteristics', AIAA, Aerospace Sciences Mtg, 1997.

[3] D.C. Wilcox, 'Turbulence Modeling for CFD', DCW Industries, Inc.; La Canada, California 91011, November 1994.

[4] A. Johansson, 'Engineering Turbulence Models and Their Development, with Emphasis on Explicit Algebraic Reynolds Stress Models', Department of Mechanics, KTH, SE-100 44 Stockholm, Sweden.

[5] R.R. Ramsay, J.M. Janiszewska and G.M. Gregorek, 'Wind Tunnel Testing of Three S809 Aileron Configurations for use on Horizontal Axis Wind Turbines', The Ohio State University Columbus, Ohio.

[6] N.N. Sorensen, J.A. Michelsen and S. Schreck, 'Navier-Stokes Predictions of the NREL Phase-IV Rotor in the NASA Ames 80-By-120 Wind Tunnel'. AIAA-2002-0031, ASME Wind Energy Symposium, pp. 94 - 105, 2002.

[7] R. Dizene, A. Nouioua et L. Bourane, 'Modélisation d'un Écoulement autour d'une Pale d'Eolienne, Champs de Vitesse Oscillant', Master Energies Renouvelables, USTHB, 2010.