# Simulation de l'écoulement autour d'un profil à l'aide d'une méthode DES

Fatih Bourekba ¹,²*, Ouahiba Guerri ¹ et Nora Nait Bouda ²

¹ Centre de Développement des Energies Renouvelables, CDER
B.P. 62, Route de l'Observatoire, 16340 Bouzaréah, Algiers, Algeria

² Laboratoire de Mécanique des fluides théorique et appliqué, Faculté de Physique
Université des Sciences et de la Technologie Houari Boumediene, USTHB
B.P. 32, El Alia, 16111 Bab Ezzouar, Algiers, Algeria

(reçu le 20 Novembre 2015 – accepté le 28 Décembre 2015)

**Résumé** - Ce travail concerne l'étude du comportement d'un écoulement turbulent autour d'un profil d'aile NACA23012. Nous avons effectué des simulations numériques de l'écoulement. Elle est réalisée par le logiciel CFD qui est basé sur la méthode de volumes finis pour la résolution des équations de Navier-Stokes. Nous avons considéré l'écoulement bidimensionnel, instationnaire et incompressible dans le domaine de vitesses choisies. L'effet de la turbulence a été pris en compte en utilisant le modèle statistique du premier ordre à deux équations SST k-ω et le modèle statistique du second ordre DES (Detached Eddy Simulation).

**Abstract** - This work concerns the study of the behavior of turbulent flow around a wing profile NACA23012. We have, performed numerical flow simulations. The CFD software that is based on the finite volume method for solving the Navier-Stokes equations performs it. We considered two-dimensional flow, unsteady and incompressible in the field chosen speeds. The effect of turbulence was taken into account using the statistical model of the first order equations with two SST k-ω and the statistical model of the second order DES (Detached Eddy Simulation).

**Keywords:** CFD - Wing NACA - Stall - RANS - SST - DES.

## 1. INTRODUCTION

L'aérodynamique, dont l'étymologie évoque l'action de l'air en mouvement, est la science qui étudie les différents aspects de cette action, notamment les forces, pressions et moments qui résultent du déplacement des corps dans l'air. Cette discipline trouve son intérêt dans le domaine de l'industrie des énergies, telles que la conception des pâles d'éoliennes.

La caractérisation des profils d'aile a été un thème d'activités très important dans le domaine de l'aérodynamique. Elle consiste principalement en la détermination des coefficients de portance, de traînée et du moment aérodynamique pour différentes configurations du profil d'aile dans un écoulement de fluide.

Des travaux expérimentaux et numériques sont abordés par plusieurs chercheurs pour montrer l'influence des paramètres caractéristiques des profils sur les grandeurs aérodynamiques globales et locales. Dans ce qui suit, nous présentons quelques travaux antérieurs.

Des systèmes aérodynamiques ont été analysés en fonctionnement et en utilisant une variété de modèles de turbulence basés sur la moyenne de Reynolds dans les équations Navier-Stokes (RANS). Bien que les modèles RANS fournissent une précision raisonnable dans les applications d'ingénierie, des solutions inexactes sont souvent produites à la suite de moyenne dans le temps avec des flux instables, comme dans les

<footnote>* f.bourekba@cder.dz</footnote>
cas avec une grande séparation ou perturbations très fluctuantes. La méthode hybride RANS/LES a été utilisée comme des alternatives efficaces pour résoudre les problèmes de flux instables avec séparation de l'écoulement. Ces procédés traitent la totalité de la couche limite avec un modèle RANS et les régions séparées instables avec le LES, dans lequel le nombre moyen de points de grille pour fournir une précision raisonnable au flux de séparation instable.

Spalart et al., [1] a développé la méthode de simulation DES. Celle-ci est une méthode hybride RANS/LES et l'idée de base est d'utiliser le modèle de turbulence S-A RANS dans la modélisation d'une région et de mur près de l'échelle de la sous grille dans la région du libre de cisaillement, afin d'obtenir des résultats comparables à ceux de l'approche LES. La technique, dépend intrinsèquement de la DES tridimensionnel et de l'heure dépendant fortement des densités de grille et la première formulation de Spalart DES 97 [1] dépend de la taille de la grille d'espacement pour agir comme un modèle LES. Toutefois, le recours à l'espacement de la grille a entraîné des difficultés lors de grilles avec des tailles ambiguës, agi comme LES dans des régions indésirables, provoquant la laminarisation de la couche limite turbulente [2].

Le DDES [2] est une version améliorée du modèle de DES97 d'origine. Il assure une bonne transition de RANS aux LES. Une formation de mélange, similaire à celle utilisée par Menter et al., [4] pour le modèle SST, a été introduite pour limiter l'ampleur de la longueur de DES97. Im et al., [5] a fait une simulation DDES du phénomène de décrochage aérodynamique où DDES a démontré une excellente capacité de simuler le décrochage dynamique pour un débit bloqué de profil NACA0012.

Gritskevich et al., [6] a proposé une stratégie de couplage pour les DDE, ainsi que les résultats numériques qui se basent sur le modèle de turbulence SST k-ω. Des résultats DDES base-SST étaient comparés avec les données expérimentales.

Des applications d'intérêt de l'ingénierie, comme pour les pales du rotor et des aubes de turbine éolienne, puisque leurs caractéristiques aérodynamiques sont affectés par le processus de transition. Des modèles RANS qui ignorent les transferts laminaires/transitions turbulents, peuvent conduire à une prévision de la force de traînée visqueuse [7]. En outre, la transition turbulente est fortement liée à la bulle de séparation laminaire, en particulier pour les faibles nombres de Reynolds [8].

## 2. DEFINITION DU PROBLEME

L'écoulement autour du profil à une incidence donnée crée des forces de surpression sur l'intrados et des forces de dépression sur l'extrados. Le point de pression minimum partage la répartition de pression sur l'extrados en deux zones. Une première région à gradient de pression favorable où l'écoulement est accéléré, suivi d'une région à gradient de pression défavorable où l'écoulement est retardé. Les forces de pression couplées aux effets de la viscosité engendrent sur le profil une force résultante qui peut être décomposée en une force de portance et une force de traînée.

### 2.1 Coefficients aérodynamiques

Les différents coefficients aérodynamiques sont définis à partir de la pression dynamique:

$$q_{\infty} = \frac{1}{2} \rho \times U_{\infty}^{2} \qquad (1)$$

Les coefficients de pression sont définis en tout point de l'aile. Il est exprimé par la relation suivante:
$$C_p = \frac{p_i - p_\infty}{q_\infty} \qquad (2)$$

Pour une aile complète, on définit les coefficients par rapport à la surface S de l'aile. Le coefficient de portance et le coefficient de traînée sont exprimés par ces expressions.

$$C_L = \frac{L}{q_\infty \times S} \qquad (3)$$

$$C_D = \frac{D}{q_\infty \times S} \qquad (4)$$

# 3. APPROCHE NUMERIQUE

## 3.1 Equations du l'écoulement

Les équations régissant l'écoulement sont les équations moyennées de Navier-Stokes. On considère dans cette étude un fluide newtonien, incompressible à propriétés physiques constantes et l'écoulement bidimensionnel en moyenne.

Equation de continuité

$$\frac{\partial U}{\partial x_1} + \frac{\partial V}{\partial y} = 0 \qquad (5)$$

Equation de conservation de la quantité de mouvement

$$\frac{\partial U_i}{\partial x} + U_j \frac{\partial U_i}{\partial y_j} + V_j \frac{\partial V_i}{\partial y} = -\frac{1}{\rho} \frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left( v \frac{\partial U_i}{\partial y_j} - \overline{u_i u_j} \right) \qquad (6)$$

$\overline{u_i u_j}$, est un tenseur du second ordre dont les éléments sont les tensions de Reynolds. La fermeture de ce système d'équations (5, 6) est assurée par la modélisation des tensions de Reynolds.

## 3.2 Modèle mathématique

Nous avons dans ce contexte le modèle de turbulence du premier ordre à deux équations SST $k-\omega$ qui prédit convenablement les décollements. C'est un modèle hybride qui conjugue les avantages respectifs de deux modèles à bas-Reynolds, $k-\omega$ en proche paroi, et haut-Reynolds $k-\epsilon$ loin de la paroi.

Les 2 équations à résoudre sont: l'équation de l'énergie cinétique turbulente k et le taux de dissipation spécifique $\omega$ exprimés ci-dessous.

$$\rho \frac{\partial k}{\partial t} + \frac{\partial}{\partial x_i} (\rho k U_i) = \frac{\partial}{\partial x_j} \left( (\mu + \sigma_k \mu_t) \frac{\partial k}{\partial x_j} \right) + P_k - \beta * \rho k \qquad (7)$$

$$\rho \frac{\partial \omega}{\partial t} + \frac{\partial}{\partial x_i} (\rho \omega U_i) = \frac{\partial}{\partial x_j} \left( (\mu + \mu_t) \frac{\partial \omega}{\partial x_j} \right) + \alpha \frac{\omega}{k} P_k - \beta \rho \omega^2 \qquad (8)$$

$$\omega = \varepsilon / k \qquad (9)$$
$$P_k = \tau_{ij} \frac{\partial U_i}{\partial x_j} \qquad (10)$$

$$\tau_{ij} = \mu_t \times \left( 2S_{ij} - \frac{2}{3} \frac{\partial \mu_k}{\partial x_k} + \delta_{ij} \right) - \frac{2}{3} \rho k S_{ij} \qquad (11)$$

$$S_{ij} = \frac{1}{2} \left( \frac{\partial \mu_i}{\partial x_j} + \frac{\partial \mu_j}{\partial x_i} \right) \qquad (12)$$

Dans ce modèle, la viscosité turbulente est définie par la relation suivante:

$$\mu_t = \rho \times k / \omega \qquad (13)$$

Comme précédemment, les constantes sont celles des couches limites atmosphériques neutres.

$$\alpha = 0.3706, \beta = 0.0275, \beta^* = 0.033, \sigma_w = 0.5, \alpha_k = 0.5,$$

Le modèle DES dans Fluent réalise la transition entre RANS et LES par une comparaison de l'échelle de longueur turbulente avec l'espacement de la grille. Le modèle sélectionne les limites de deux modes et passe ainsi entre le RANS et LES en remplaçant $\epsilon$ dont l'équation de k :

$$\epsilon = \frac{k^{3/2}}{L_{k-\omega}} \rightarrow \epsilon = \frac{k^{3/2}}{\min(L_{k-\omega}, C_{DES}\Delta_{\max})} \qquad (14)$$

Une fois que le modèle a sélectionné l'espacement de la grille que le minimum, le modèle fonctionne en modeLES. L'espacement de la grille entre explicitement dans le modèle DES. Cela peut affecter la solution RANS dans les régions, où la grille est comprise entre RANS et LES résolution ou lorsque l'instabilité de l'écoulement n'est pas assez fort pour générer des structures. Un autre problème à considérer avec DES est le problème de la grille séparation induite.

$$C_{DES} = (1 - F_1) C_{DES}^{k-\epsilon} + F_1 C_{DES}^{k-\omega} \qquad (15)$$

$$C_{DES}^{k-\epsilon} = 0.61 C_{DES}^{k-\omega} = 0.78$$

Il se produit si la grille pour un écoulement de couche limite de la paroi mitoyenne est raffinée jusqu'à un point où le limiteur DES devient actif et la solution affecte RANS. Cependant, dans de telles situations, l'instabilité de l'écoulement n'est pas assez forte pour équilibrer le contenu RANS réduite par la turbulence résolu. Ce sera typiquement lieu à une séparation de l'écoulement artificiel à l'emplacement de la grille de raffinement.

## 3.3 Méthode de résolution numérique

Nous avons recensé dans la littérature deux types de maillage pour la géométrie d'un profil, le maillage en 'O' et le maillage en 'C'. Nous avons opté pour le deuxième représenté sur la figure 1.

C'est un maillage structuré non uniforme qui comporte près de 90400 mailles en prenant le soin de le raffiner près de la paroi. La distance à la paroi de la première cellule est de l'ordre du $10^{-4}$.
<figure><img src="image_1.png" /><figcaption>Fig. 1: Maillage appliqué</figcaption></figure>

Le calcul est mené par le code CFD Fluent qui utilise la méthode des volumes finis. L'algorithme SIMPLE est appliqué pour résoudre les équations couplées pression-vitesse. Dans Fluent, toutes les variables sont calculées aux centres des mailles. Le calcul des flux sur les faces nécessitent toutefois de connaître la valeur de la variable sur cette face. Celle-ci est déterminée par interpolation entre les valeurs voisines connues. La discrétisation spatiale des termes des équations de transport est basée sur un schéma Upwind de second ordre. Les critères de convergences des résidus normalisés des variables dépendantes sont fixés à $10^{-5}$.

# 4. RESULTATS

De la confrontation des résultats numériques, les paramètres de calcul, maillage, nombre de mailles, taille des mailles et modèles de fermeture [9] ont pu être réajustés afin d'améliorer la prédiction des coefficients aérodynamiques. Nous présentons dans ce travail les résultats obtenus pour le nombre de Reynolds de $6.85 \times 10^4$ comparés aux résultats expérimentaux, déjà effectués précédemment [10].

Les essais expérimentaux ont été réalisés en soufflerie subsonique $EA \ 10^3$ à circuit ouvert aux laboratoires LMFTA USTHB. La veine d'essai de section carrée $30 \times 30$ cm, et d'une longueur de 1 m, a été réalisée en altuglas transparent.

## 4.1 Evolution du coefficient de pression

L'évolution du coefficient de pression, sur le profil NACA 23012, obtenu par les deux approches est représentée sur la figure 2 pour le nombre de Reynolds $6.85 \times 10^4$.

<figure><img src="image_2.png" /></figure>
De cette figure, on trouve un résultat classique, relatif à une augmentation du coefficient de pression, quand l'angle d'incidence augmente jusqu'à l'angle de décrochage. Par ailleurs, la comparaison avec les données expérimentales montre que l'approche DES est supérieure à l'approche SST et présente généralement un bon accord avec les données expérimentales.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Evolution du coefficient de pression autour du profil à différents angles d'incidences (5°, 10°, 17.5°, 20°)</figcaption></figure>

## 4.2 Coefficient de portance

L'évolution de la portance en fonction de l'angle d'incidence du profil, à Re = 6.85 × 10⁴, est représentée sur la figure 3. Nous avons regroupé les résultats obtenus par DES comparés à ceux prévus par SST. La variation de portance avec l'angle d'attaque α est linéaire sur un intervalle des angles d'incidence allant approximativement de -5° jusqu'à 10°, le décrochage se manifestant par une diminution du coefficient à un angle avoisinant 12.5°.

La différence entre les résultats présentés est notable après le décrochage. Aussi, les solutions qui ne présentent pas de SST. La tendance des courbes calculée $C_L$ de DES sont en bon accord avec les données expérimentales. Cependant, DES prédit un point de décrochage un peu plus retardé que celle observée dans l'expérience et il ne parvient pas à capturer la crête de coefficient de portance pour le cas étudiés.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Evolution du coefficient de portance en fonction de l'angle d'incidence</figcaption></figure>
## 4.3 Coefficient de trainée

L'évolution du coefficient de trainée $C_D$ sur le profil en fonction de l'angle d'inclinaison est présentée sur la figure 4. On remarque que pour les faibles angles d'attaque de l'écoulement, le coefficient de trainée est faible; mais augmente lorsque $\alpha$ augmente. De l'effet de trainée de forme est rajoutée la trainée de frottement qui s'amplifie avec l'épaississement de la couche limite.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Courbe de coefficient de traine pour différents angles d'incidence et pour différents calculs</figcaption></figure>

## 4.4 Visualisation de l'écoulement

La figure 5 représente un aperçu général de l'écoulement autour de l'aile à des incidences différentes. On constate qu'avec l'augmentation de l'angle de l'incidence, la zone de recirculation augmente, le décollement est marqué par une zone de recirculation.

Pour un angle d'incidence égale à 10°, il y a apparition d'une zone de séparation dans la DES, cette zone prend de l'ampleur avec l'augmentation de l'angle d'incidence.

D'autre part, nous remarquons l'apparition, au-delà de 15°, d'un deuxième tourbillon prenant naissance au bord de fuite. Ces deux tourbillons sont en sens opposés.

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

a- STT AOA=10°

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

b- DES AOA=10°
<figure><img src="image_8.png" /><figcaption>Fig. 5: Lignes de courant et pour les grands angles d'incidence pour les deux approches -a- SST, -b- DES</figcaption></figure>

## 5. CONCLUSION

Cet article étudie numériquement l'écoulement turbulent séparé instable autour d'un profil d'aile NACA. L'écoulement est dans le régime relativement faible du nombre de Reynolds sur la base de la longueur de corde de profil. Le URANS et les approches les plus avancées DES sont utilisés. RANS est couplé avec un modèle de turbulence, à savoir la pointe du modèle de transition SST $k-\omega$ et le DES.

Les simulations DES ont été appliquées avec succès pour évaluer l'écoulement externe entièrement turbulent. La DES exploite les avantages de la modélisation RANS à l'intérieur de la couche limite et le LES en dehors de la couche limite turbulente. Les résultats ont généralement montré que d'excellente solutions apportées par la simulation DES pour l'écoulement de réveil instable, mais non pas pour les écoulements transitoires parce que la DES est fondamentalement basée sur une hypothèse complètement turbulente.

Le RANS est toujours le cheval de travail du logiciel CFD pour les applications de flux industriels où il est intuitif à utiliser des modèles de turbulence simples, tels que le modèle SST $k-\omega$ et il est numériquement moins cher que les méthodes plus complexes, comme la simulation de DES et la simulation LES. Toutefois, la précision des prédictions numériques utilisant des modèles de turbulence simples doit être quantifiée. Le cas étudié concerne le décrochage, les charges aérodynamiques
instantanées exercées sur le profil et la structure de l'écoulement autour de profil
associée.

En outre, les solutions DES révèlent un modèle d'écoulement périodique, comme le
RANS effectue pour les faibles angles d'attaque. En conclusion, cette étude sera suivie
par la caractérisation du phénomène de décrochage dynamique, un phénomène que l'on
rencontre particulièrement dans les éoliennes à axe vertical.

REFERENCES

[1] P.R. Spalart, W.H. Jou, M. Strelets and S.R. Allmaras, ‘Comments on the Feasibility of LES for Wings, and on a Hybrid RANS/LES Approach’, Advances in DNS/LES, 1st AFOSR, International Conference on DNS/LES, H. Columbus, Aug. 4–8, 1997.

[2] P.R. Spalart, S. Deck, M. Shur, K. Squires, M. Strelets and A. Travin, ‘A New Version of Detached-Eddy Simulation, Resistant to Ambiguous Grid Densities’, Theoretical and Computational Fluid Dynamics, Vol. 20, N°3, pp. 181–195, 2006.

[3] F.R. Menter, ‘Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications’, AIAA-Journal, Vol. 32, N°8, pp. 269-289, 1994.

[4] F.R. Menter and M. Kuntz, ‘Adaptation of Eddy-Viscosity Turbulence Models To Unsteady Separated Flow Behind Vehicles, The Aerodynamics of Heavy Vehicles: Trucks, Buses, and Trains’, Springer Berlin Heidelberg, pp. 339–352, 2004.

[5] H.S. Im and G.C. Zha, ‘Delayed Detached Eddy Simulation of a Stall Flow over NACA0012 Airfoil using High Order Schemes’, AIAA Paper, 2011-1297, 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, Florida, 4-7 January 2011.

[6] M.S. Gritskevich, A.V. Garbaruk, J. Schütze and F.R. Menter, ‘Development of DDES and IDDES Formulations for the k-ω Shear Stress Transport Model’, Flow, Turbulence Combustion, Vol. 88, N°3, pp. 431–449, 2012.

[7] J. Johansen, ‘Prediction of Laminar/Turbulent Transition in Airfoil Flows’, Risø-R-987(EN), Denmark, 1997.

[8] J. Windte, U. Scholz and R. Radespiel, ‘Validation of the RANS-Simulation of Laminar Separation Bubbles on Airfoils’, Aerospace Science and Technology, Vol. 10, pp. 484–494, 2006.

[9] P. Wernert, W. Geissler, M. Raffel and J. Kompenhans, ‘Experimental and Numerical Investigations of Dynamic Stall on a Pitching Airfoil’, AIAA Journal, Vol. 34, N°5, pp. 982–989, 1996.

[10] F. Bourekba, ‘Etude d'un Ecoulement Turbulent autour d'une Aile NACA', Mémoire de Magister, Université des Sciences et de la technologie Houari Boumediène, 2012.

[11] W. McCroskey, K. McAlister, L. Carr and S. Pucci, ‘An Experimental Study of Dynamic Stall on Advanced Airfoil Sections', Vol. 1: Summary of the Experiment. NASA TM, 1982.

[12] T. Lee and P. Gerontakos, ‘Investigation of Flow over an Oscillating Airfoil', Journal of Fluid Mechanics, Vol. 512, pp. 313–341, 2004.
[13] J.G. Leishman, 'Dynamic Stall Experiments on the NACA 23012 Aerofoil', Experiments in Fluids, Vol. 9, N°1, pp. 49 - 58, 1990.

[14] S. Wang, D.B. Ingham, L. Ma, M. Pourkashanian and Z. Tao, 'Numerical Investigations on Dynamic Stall of Low Reynolds Number Flow Around Oscillating Airfoils', Computers & Fluids, Vol. 39, pp. 1529 – 1541, 2010.

[15] S. Wang, D.B. Ingham, L. Ma, M. Pourkashanian and Z. Tao, 'Turbulence Modeling of Deep Dynamic Stall at Relatively Low Reynolds Number', Journal of Fluids and Structures, Vol. 33, pp. 191 – 209, 2012.

[16] A. zanotti, R. Nilifard, G. Gibertini, A. Guardone and G. Quaranta, 'Assessment of 2D/3D Numerical Modeling for Deep Dynamic Stall Experiments', Journal of Fluids and Structures, Vol. 51, pp. 97 - 115, 2014.

[17] Dong-Ha Kim, Jo-Won Chang; 'Low-Reynolds-Number Effect on the Aerodynamic Characteristics of a Pitching NACA0012 Airfoil', Aerospace Science and Technology, Vol. 32, pp. 162 - 168, 2014.

[18] Jeong Hwan Sa, Soo Hyung Park, Chang Joo Kim and Jung Keun Park, 'Low-Reynolds Number Flow Computation for Eppler 387 Wing Using Hybrid DES/Transition Model', Journal of Mechanical Science and Technology, Vol. 29, N°5, pp. 1837 - 1847, 2015.