Etude d'OpenFoam






Index

    1. Présentation du logiciel

    2. Cas test : Profil NACA0012

    3. Architecture du logiciel

    4. Solveurs utilisés

    5. Evaluation du modèle Spalart-Allmaras

            a) Conditions initiales et conditions limites
            b) Influence du maillage
            c) Comparaison du 1er et 2nd ordre spatial
            d) Effet du schéma temporel

    6. Evaluation du modèle k-Epsilon

            a) Calculs stationnaires
            b) Calculs instationnaires
            c) Comparaison des polaires
            d) Conclusion

    7. Conclusion


1. Présentation du logiciel


    OpenFoam (Open Field Operation and Manipulation) est un des premiers outils, écrit dans le langage C++, permettant de résoudre des problèmes de mécanique des milieux continus et incluant des modèles de mécanique des fluides numériques (CFD).

    Son développement a débuté à l’Impérial Collège de Londres dans les années 1980. Le but était de trouver une plateforme de simulation générale plus puissante et plus modulable que celles déjà existantes à l’époque et qui utilisaient le langage Fortran. Il est ensuite devenu Open-source sous licence GNU General Public License en 2004.

    OpenFoam rivalise avec la majorité des logiciels commerciaux de CFD (Computational Fluid Dynamics) et est applicable à une large gamme de problèmes. Il contient des solveurs standards permettant de simuler des problèmes d’écoulements compressibles, incompressibles ou multiphasiques mais aussi des problèmes de combustion, d’échange de chaleur, ou d’électromagnétisme par exemple. L’un des avantages de ce logiciel est la relative facilité offerte à l’utilisateur de créer de nouveaux solveurs adaptés à son problème.

    De plus, l’utilisateur a le choix entre une large gamme de logiciels compatibles avec OpenFoam pour le pre ou post-processing des données. Ainsi, OpenFoam est capable d’utiliser un maillage créé sous Gambit par exemple. Dans notre cas, nous utiliserons un maillage créer pour Fluent et que l'on exportera grâce à une fonction spécifique (Fluent3DMeshToFoam). Le post-processing et la réalisation des animations sera réalisé à l'aide de Paraview.

Le site officiel est à l’adresse suivante : http://www.openfoam.com/download/ Le logiciel peut y télécharger librement ainsi que différents tutoriels et le manuel d’utilisation.

Un autre avantage de ce logiciel est le dynamisme de sa communauté. Le principal forum http://www.cfd-online.com/Forums/openfoam/ est entièrement consacré à ce logiciel et peut s’avérer d’une aide précieuse.



2. Cas test : Profil NACA0012

    Les maillages utilisés ici ont été créé sous Fluent, puis convertis au format OpenFoam grâce à la fonction fluent3DMeshToFoam. Ce sont des maillages conçus pour réaliser des simulations 2D, cependant OpenFoam ne travaillant qu'avec des maillages 3D ils comportent tous deux une cellule dans la troisième direction.

    Les deux maillages seront nommés respectivement maillage original et maillage extended. Le maillage extended reprend le maillage original en étendant la distance entre le profil et la frontière afin de vérifier que les effets induits par la frontière n'influencent pas les résultats. Le maillage extended compte 50 cordes alors que la maillage initial n'en compte que 15.


fig 1. Maillage original


fig 2. Maillage sur le profil


    Dans ce projet, nous allons étudier l'impact d'un écoulement sur un profil NACA0012 en incidence

        Nombre de Reynolds : Re = 100.000

        Corde du profil : 2.5 10-5 m




3. Architecture du logiciel

    Quelque soit le solveur utilisé, l'architecture du répertoire est redondante. Par la suite, on se servira de ce schéma afin d'expliquer les différentes étapes du travail effectué sur le logiciel.


fig 3. Architecture d'un répertoire de calcul
Sous OpenFOAM, le choix de l’application se fait par le biais du fichier system/controlDict. Ce fichier permet également de spécifier le pas de temps, l’instant de départ et de fin et l’intervalle d’écriture des fichiers de sortie.
Dans le fichier controlDict, le paragraphe functions() permet d’activer des fonctions complémentaires appelées function objects, notamment le calcul des efforts sur une frontière (forceCoeffs).

Les propriétés du fluide sont spécifiées dans le fichier constant/transportProperties.

Le fichier de conditionnement permettant le choix des schémas de discretisation, est le fichier system/fvSchemes. Ce fichier est décomposé en paragraphes, correspondant à un opérateur mathématique (divergence, gradients...), avec à l'intérieur de chacun, la possibilité du choix de discrétisation de chaque terme des équations de conservation.


4. Solveurs utilisés


    1. SimpleFoam
   
SimpleFoam est un solveur de turbulence stationnaire incompressible, qui peut utiliser plusieurs modèles de turbulence. Ce solveur est utilisé pour les angles d'incidences faibles (< 10-12°), pour lesquels le décollement n'est pas important.

    2. PisoFoam

PisoFoam est un solveur turbulent instationnaire incompressible. Ce solveur sera utilisé pour les grands angles d'incidence impliquant un décollement important, et donc de forts effets instationnaires.

    3. Modèles étudiés

Dans ce projet nous effectuerons différentes simulations avec les modèles suivants :


Spalart-Allmaras est un modèle de transport à une équation :

 

La viscosité turbulente est définie par :
Relations auxiliaires liées au modèle :

Valeur des coefficients utilisés :

Condition limite à la paroi : on impose une viscosité turbulente nulle à la paroi.

Contrairement au précédent, celui-ci est un modèle à deux équations de transport.
La viscosité turbulente est donnée localement par:
Modélisation des équations de transport :

Coefficients standard :

Le modèle Launder-Sharma est un modèle k-Epsilon adapté aux bas Re.

Les fonctions d'amortissements fµ, f1 et f2 ainsi que les termes D et E sont actifs uniquement près des parois solides et permettent de résoudre k et ε jusqu'à la sous-couche visqueuse.

Remarque :

2 autres options peuvent également être envisagé pour fµ,
Conditions limites :


5. Evaluation du modèle Spalart - Allmaras

Cette partie est consacrée aux différents tests effectués sur le modèle Spalart-Allmaras, avec les solvers simpleFoam et pisoFoam, afin de constater son potentiel et de comparer les résultats obtenus à des expériences.

        a) Conditions initiales et conditions limites

Rappel du contenu des fichiers :

Conditions initiales et limites : fichier 0
Propriétés du fluide : fichier constant
Schémas de discrétisation : fichier system
Voir le contenu des fichiers pour : Cas stationnaire : simpleFoam

Remarque :
Le cas instationnaire, traité avec pisoFoam, ne différe que par légèrement du fait de l'utilisation d'une disrétisation au niveau temporel (cf. constant/fvSchemes)

        b) Influence du maillage

Dans un premier temps, nous avons étudié l'effet du maillage sur les résultats des simulations numériques. Pour cela, quelques points ont été calculé avec des schéma de discrétisation au 1er et au 2nd ordre spatial, à partir des deux maillages mis à notre disposition : le maillage original et extended.

Résultats des calculs au 1er ordre :


On constate donc que le maillage n'a quasiment aucun effet sur un calcul au premier ordre. Ceci s'explique probablement par la très forte dissipation de ce genre de modèle.
Seul un léger écart dans l'amplitude des courbes du coefficient de portance Cl est constaté aux grandes incidences, mais ce phénomènes est totalement négligeable.

Résultats des calculs au second ordre :


On remarque qu'au 2nd ordre, la courbe correspondant au coefficient de portance est lissée, cependant les valeurs obtenues avec le maillage extended ne diffèrent que peu des valeurs du maillage original.
Pour ce qui est du coefficient de traînée, on constate un changement effectif uniquement pour une incidence de 15°, ce qui ne permet toujours pas de marquer suffisament l'augmentation de traînée due au décrochage par rapport à l'expérience. On peut donc négliger cette influence.

Conclusion :

Même si quelques changements sont observés au 2nd ordre, ils ne sont pas significatifs, et le gain apporté par le maillage extended est jugé négligeable par rapport à l'augmentation du temps de calcul qu'il occasionne du fait d'un plus grand nombre de cellules.
La suite de l'étude sera donc menée avec le maillage original, jugé suffisant pour évaluer les possibilités du code dans notre configuration.

Grâce à ces premières simulations, on peut déjà constater une différence considérable entre les résultats obtenus au premier et au second ordre.
La partie suivante traite de cette remarque.
 

        c) Comparaison du mdèle Spalart-Allmaras pour des schémas spatiaux du 1er et 2nd

Pour étudier le phénomène mis en évidence précédemment, ces calculs ont été approfondis et comparé aux résultats expérimentaux. Cette partie explicite les résultats obtenus.

Etude d'un cas stationnaire

On étudie ici un cas à moyenne incidence : 10°. On relève plusieurs paramètres, visualisés à l'aide du logiciel Paraview, afin de vérifier la cohérence des résultats. Ceci n'est qu'une première approche afin de regarder que les résultats obtenus ne sont pas aberrants. Cette approche ne peut en aucun cas nous permettre de porter un jugement quant à la validité du code. Ce travail sera mené un peu plus loin.

Tracé de la vitesse :
Tracé de la pression :

Tracé de la viscosité turbulente :
Première analyse :
  • Vitesse : l'angle d'incidence 10° est théoriquement assez proche de l'angle de décrochage, il est donc logique d'observer une zone de recirculation sur l'extrados, qui doit cependant rester faible. On peut également observé un point d'arrêt au bord d'attaque et une vitesse supérieure au niveau de  l'extrados par rapport à l'intrados. 
  • Pression : pour 10° d'incidence la portance doit être relativement importante, on retrouve bien une pression supérieur à l'intrados par rapport à l'extrados. Ceci assure donc bien une portance positive de l'aile.
  • Viscosité : la viscosité turbulente est bien supérieure dans les zones de recirculation où la turbulence est plus importante.

Cette première analyse nous permet donc de voir que ces résultats ne sont pas aberrant, et nous autorise à approfondir l'étude de ce modèle au 1er ordre. Il est cependant possible que certains changements aient besoin d'être effectués pour améliorer les résultats notamment pour les fortes incidences . Le phénomène de décrochage étant instationnaire, les résultats stationnaires restent limités dans leur représentation de la réalité.



Etude d'un cas instationnaire

Etude d'un cas à incidence élevée : 20°.

Tracé de la vitesse :
Tracé de la pression :

Tracé de la viscosité turbulente :
Première analyse :
  • Vitesse : l'angle 20° d'incidence correspond théoriquement à un cas décroché, on s'attend donc à trouver une zone de recirculation importante, que l'on observe bien ici.
  • Pression : la pression se concentre au point d'arrêt, et le gradient de pression paraît faible, ce qui correspondrait bien à une chute de portance. Cependant ce graphique ne permet en aucun cas de conclure sur ce point.
  • Viscosité : même remarque que précédemment, cette fois l'augmentation de la viscosité turbulente avec la vitesse est évidente.

Comme pour le cas stationnaire, cette première analyse montre des résultats qui semblent cohérent avec la théorie.
Ces premiers cas testés nous permettent de nous lancer dans les calculs avec un a priori positif.


Comparaison des modèles au 1er et 2nd ordre Spatial

Remarque :
Pour les courbes obtenues ci-dessous, le solver stationnaire simpleFoam a été utilisé jusqu'à un angle d'incidence de 10°, ensuite, le phénomène instationnaire étant prédominant, le solveur pisoFoam a pris le relais.


    Sur la 1er courbe représentant le coefficient de portance (Cl), on remarque que le 1er ordre ne capte pas du tout le phénomène de chute de portance.
    Le 2nd ordre capte mieux le phénomène de chute de portance, cependant le décrochage est plus tardif par rapport au résultats expérimentaux. L'amplitude de la chute de portance constatée est également moindre.
    Enfin, on remarque que le point situé à une incidence de 20°, si on tient compte du décalage entre les courbe, semble mal capter le palier situé entre 15 et 17° sur la courbe expérimentale, ce qui se traduit par une réaugmentation de la portance.
    Aux grands angles d'incidence, le Cl discrétisé en espace au 2ème ordre, se rapproche de l'expérimental, tandis que le Cl discrétisé au 1er ordre semble poursuivre une progression linéaire.

    La courbe 2, représentant le coefficient de traînée (Cd), ne permet de conclure quant à la validité d'un des modèles. En effet le phénomènes d'augmentation de traînée due au décrochage, et bien visible sur les résultats expérimentaux, n'apparaît pas clairement sur les modélisations au 1er ou 2nd ordre.

Conclusion :

    Le 1er ordre semble beaucoup trop dissiper, ce qui lisse les courbes obtenues. Les résultats donnés par ce modèle sont très éloignés de la courbe expérimentale et le modèle au 1er ordre spatial n'est à retenir. Ceci était connu mais cette étude nous a servie de base à notre travail.
    Le 2nd ordre donne de meilleur résultats, et capte la physique de la chute de traînée. Cependant, la courbe obtenue comporte un décalage angulaire, et quelques problèmes persistent au niveau de certains angles.
    Au vu des résultats, le modèle de 2nd ordre spatial a donc été retenu.


        d) Effet du modèle temporel

Le modèle temporel ayant été comparé en même temps que le modèle spatial, il a été décidé de se reprendre la configuration initiale avec des schémas spaciaux du 1er ordre afin de ne changer qu'un seul paramètre à la fois.

Voici les résultats obtenus :


En comparant les deux courbes obtenues, il semble qu'il n'y ait aucune différence notable entre les différents schémas temporels.
On peut seulement noter une différence au niveau du point à 20° d'incidence sur le coefficient de traînée, qui paraît plus proche des résultats expérimentaux. Cependant ce résultats doit être fortement nuancé par le fait que le schéma de Crank-Nicholson a beaucoup de mal à converger, et donc ce point doit être mis entre parenthèse par rapport aux problème rencontrés avec ce schéma. (Il faut de plus savoir que ce schéma a été testé pour une incidence de 15° mais il a totalement divergé).

On peut donc conclure que les différents schémas temporels au 2nd ordre n'apporte rien par rapport au schéma d'Euler.
Au vue des résultats précédent, il est cependant nécessaire de noter qu'une analyse utilisant les schémas spatiaux de 2ème ordre est fortement souhaitable voire nécessaire pour conclure à ce sujet. Le manque de temps ne nous a pas permis de la réaliser.


6. Evaluation du modèle k-Epsilon Launder-Sharma

    L'objectif étant de tester les capacités du logiciel OpenFoam, nous avons également testé un modèle à deux équations. De nombreux modèles sont disponibles sous OpenFoam, cependant notre maillage n'étant pas un maillage loi de paroi, nous avons choisi le seul modèle Bas Reynolds disponible, en l'occurrence le modèle k-epsilon Launder-Sharma.

Calculs stationnaires

Comme pour le modèle précédent, l'étude a d'abord été menée en stationnaire (Toujours à l'aide du solveur SIMPLEFOAM).


On trace l'énergie cinétique et la dissipation pour une incidence de 10° :
Tracé de l'énergie cinétique turbulente, k :
Tracé de la dissipation, Epsilon :

Le tracé de l'énergie cinétique nous permet de constater que l'on a bien une élévation de l'énergie dans les zones de forte turbulence.

La dissipation semble bien se concentrer au niveau des forts gradients.


Les courbes de traînée et de portance ont été tracées pour des incidence allant de 0 à 30 degrés.
Cependant comme on pouvait s'y attendre, ces calculs stationnaires rencontrent des problèmes de convergences à partir de l'angle de décrochage. On retrouve ci dessous les évolutions des coefficients de traînée et de portance au cours de la simulation pour un angle de 10° et de 11°.



Le critère de convergence de ces valeurs a été fixé à 10-2 .

Il est donc clair qu'au delà d'un angle de 10° d'incidence, ces calculs ne satisfont pas cette exigence. Les calculs instationnaires sont donc inévitables et ont aussi été réalisé avec le solveur instationnaire PISOFOAM.


Calculs instationnaires

Tracé de la vitesse, pour un angle de 23° :



On sait que le nombre de strouhal de cet écoulement est de 0,12.
Par définition :f : fréquence du phénomène,
C : la longueur de la corde du profil (0,025cm),
U : Vitesse de l'écoulement (64,36m/s)

On déduit ainsi la fréquence du phénomène, f = 310 Hz

Cette fréquence doit se retrouver sur les simulations instationnaires. Pour cela, nous avons réalisé une transformation de Fourier à l'aide de Matlab. La fréquence d'échantillonage a été choisie à 10000 hz, respectant ainsi le critère de Shannon mais aussi un échantillonage de l'ordre de 30 points par période pour assurer une bonne analyse du signal.

Les résultats suivants ont été obtenus pour un angle d'incidence de 16°. On retrouve bien l'ordre de grandeur de la fréquence calculée précédemment avec une fréquence de 260 Hz. Ceci est donc un premier élément positif pour valider l'utilisation de ce solveur.



Les résultats sont présentés ci-dessous :

Comme on peut le voir ci-dessus, les simulations n'ont pas été menés pour des angles inférieur à 8° pour limiter le nombre de calcul. Cependant, elles auraient pu l'être et on devrait retrouver les même résultats, car des calculs instationnaires sont aussi adaptés à des calculs stationnaires mais la réciproque n'est pas vrai.


Résultats :

Une première remarque reside dans le faire que les calculs instationnaires sont à présent convergés même si souvent une moyenne statistique a été nécessaire, On remarque que ce modèle donne une sous-estimation de la traînée encore plus forte qu'avec des calculs stationnaires. On retrouve aussi une meilleure prédiction du profil de portance. L'angle de décrochage se rapproche de l'angle observé expérimentalement et la chute de portance est plus importante.

Comparaison des polaires

Pour juger au mieux des résultats et de la validité de ce modèle, la courbe suivante combine les résultats stationnaires jusqu'à une incidence de 10° et les simulations instationnaires pour des angles supérieurs.


Conclusion

En comparaison avec le modèle précédent, les résultats semblent plus proche de la réalité :
  • La chute de portance est bien transcrite et l'angle de décrochage est proche de l'expérimental.
  • Pour les petits angles d'incidence, la traînée est sur-estimée alors que dans les grandes incidences, elle est sous-estimé. On obtient donc une courbe de traînée lissée. On sait que lors du décrochage, la traînée augmente, alors que dans ce cas, sans les résultats expérimentaux, la courbe ne présente aucune rupture de pente permettant de déterminer un angle de décrochage.
Ces résultats peuvent cependant paraître étrange, car un schéma du premier ordre même avec un modèle à 2 équations ne doit pas être plus précis que le Spalart-Allmaras au second ordre. Il aurait donc fallut faire tourner ce modèle au second ordre pour une meilleure comparaison mais les nombreux tests, tant pour les conditions limites que les schéma n'ont pas permis de réussir à faire tourner ces calculs (quelque soit l'angle).


7. Conclusion

Les différents modèles testés montrent que le logiciel OpenFoam démontre un potentiel certain quant à la modélisation numérique de phénomènes aérodynamiques simples. Cependant il est nécessaire d'effectuer certains changements afin de s'adapter au cas étudié ainsi qu'au maillage utilisé.
Au vue de nos résultats, la configuration la plus adaptée semble un modèle k-epsilon launder-sharma discrétisé au premier ordre. Cependant nous sommes surpris de constater que ces derniers sont meilleurs qu'un modèle Spalart-Allmaras deuxième ordre. Il serait donc nécessaire de faire converger les calculs du modèle k-epsilon deuxième ordre pour pouvoir tirer de meilleures conclusions.
En effet, le maillage en O utilisé ici a démontré certaines limites, en particulier sur le choix de certaines conditions limites inadaptées à ce maillage. Un axe de recherche pour améliorer les résultats serait de modifier le maillage afin de pouvoir imposer des conditions en entrée et en sortie distinctes( condition de type InletOutlet qui s'adapte selon les besoins de l'écoulement en imposant une condition soit  freestream, soit zerogradient) . Ici le maillage en O rassemble ces deux conditions du fait de sa géométrie et ne permet donc pas d'approfondir l'étude de manière fine.
Enfin, un nouveau solveur AeroFoam, disponible dans la version 1.6, semble adapté à la résolution d'écoulement aérodynamique, ainsi une des prochaines étapes de validation de ce code pourrait être de tester ce solveur.