OpenFOAM

Cette partie s'attache à détailler toutes les fonctionnalités que l'on a utilisées sur OpenFoam tout au long du projet.

Présentation du code


Source : http://hmf.enseeiht.fr/travaux/projnum/book/export/html/898

 

OpenFOAM, pour Open Field Operation and Manipulation, est un code de calculs CFD libre, disponible gratuitement ici. Il s'agit d'un code sans interface graphique, développé en C++.

Son caractère libre permet à ses utilisateurs de modifier autant qu'ils le souhaitent le code puisque les sources sont livrées et accessibles lors du téléchargement du code. Cela donne donc théoriquement la possibilité de résoudre n'importe quel problème physique avec n'importe quelle méthode numérique.

Le point limitant de cette apparente liberté est que ces sources sont parfois confuses dans leur organisation, il faut souvent se référer à la communauté utilisatrice d'OpenFOAM, heureusement très active, pour espérer avoir de l'aide. Finalement, pour des élèves-ingénieurs spécialisés en Mécanique des Fluides et en CFD que nous sommes, l'adaptation des premiers cas tests déjà entièrement codés nous est largement accessible. En revanche, nous serons limités par une modification toute relative des sources, dû à notre manque d'expérience en C++ mais aussi au temps imparti pour ce BEI.

 


Structure d'OpenFOAM ...

Source : http://www.openfoam.org/docs/user/userch1.php

D'un point de vue programmeur, OpenFoam est comme tout code de CFD décomposé en trois parties, le pre-processing, le noyau comportant les solveurs et le post-processing. Toutes ces parties reposent et font appel à une base composée de librairies C++. La partie pre-processing contient principalement les outils de maillage comme blockMesh ou snappyHexMesh. Le noyau du code contient tous les solveurs disponibles, et enfin le post-processing est principalement composé du logiciel de post-traitement Paraview.


D'un point de vue utilisateur, OpenFOAM, une fois installé, peut être symboliquement séparé en deux parties.

La première contient toutes les sources où les calculs sont concrètement réalisés, avec toutes les équations résolvables en CFD selon les hypothèses, modèles et types de discrétisation utilisés.

La deuxième contient ce qui est appelé les "tutoriels". Ce sont des cas tests relativement simples qui sont directement exécutables sans que l'utilisateur ait à modifier quoique ce soit. En fait, sont présents dans ces tutoriels les fichiers qui indiquent la géométrie, les équations résolues, les modèles choisis et les types de schémas de discrétisation désirés. Ils appellent les routines qui sont présentes dans la partie "source" du code où sont décrits les calculs à effectuer. En pratique, lorsqu'on débute un projet avec OpenFOAM, il est extrêmement pertinent de partir d'un de ces tutoriels pour ensuite modifier la géométrie et tous les paramètres numériques initialement présents. C'est d'ailleurs la démarche que nous allons suivre.

Remarquons que ces tutoriels sont répartis entre "solveurs". Chacun de ces solveurs correspond en fait un grand type de problème en CFD, comme par exemple les écoulements incompressibles laminaires (solveur icoFoam), les écoulements incompressibles avec maillage mobile (solveur pimpleDyMFOAM), les écoulements réactifs (reactingFoam) pour ne citer qu'eux. Cela est extrêmement pratique afin de se diriger rapidement vers le tutoriel qui ressemble le plus au problème que l'on souhaite résoudre.

Chacun de ces tutoriels est structuré de la même façon. Un tutoriel est un dossier qui se divise en 3 répertoires :

- le dossier "0" où sont décrites les conditions initiales et conditions limites des différentes grandeurs calculées

- le dossier "constant" où est décrit la géométrie du domaine avec le maillage associé. OpenFOAM utilise principalement des maillages hexaédriques

- le dossier "system" sont décrits les paramètres numériques du calcul, comme les schémas de discrétisation, les pas de temps, la fréquence d'écriture des résultats

Cette structure peut être représentée par la figure suivante :

Source : http://www.openfoam.org/docs/user/case-file-structure.php

Enfin, pour lancer un calcul OpenFOAM, il faut créer le maillage en exécutant la commande blockMesh à la source du tutoriel, puis on lance le solveur en tapant le nom du solveur dans le terminal (pimpleDyMFoam par exemple).

Ceci est une première explication de la structure d'OpenFAOM, nous détaillerons ensuite davantage les sections où notre travail s'est concentré, notamment sur la mise en place de la géométrie, la gestion du maillage mobile, ainsi que l'interpolation des résultats.

Solveurs utilisables pour un maillage mobile

Pour le cas des écoulements incompressibles, OpenFOAM est capable de gérer la résolution des équations avec un maillage mobile grâce à la méthode Arbitrary Lagrangian-Eulerian. Concrètement, le domaine est décrit par un seul maillage dont la vitesse de chaque noeud est imposée et donc connue. Les équations de Navier-Stokes à coder prennent en compte le fait que les volumes de contrôle sont en mouvement et aucune approximation supplémentaire n'est faite par rapport à un calcul avec maillage fixe.

Le solveur d'OpenFOAM permettant de réaliser de tels calculs s'appelle PimpleDyMFoam, pour Pimple Dynamic Mesh Foam.

Pour ce solveur, dans le cas des écoulements incompressibles, quatre tutoriels sont présents dans la version 2.1.0 d'OpenFOAM. Parmi ces quatre tutoriels, deux ont particulièrement retenu notre attention : le cas du movingCone ainsi que le cas du mixerVesselAMI2D.

Le premier traite l'avancée d'un cône solide dans un domaine fermé, un peu à l'idée d'un piston. Il résout une portion angulaire d'un problème axisymétrique. La gestion du maillage mobile est décrite dans le fichier dynamicMeshDict (dans le répertoire constant) et est assurée par le solveur dynamicMotionSolverFvMesh, où l'utilisateur indique que le déplacement des mailles se fait selon l'axe X. Avec ce solveur, il est aussi nécessaire d'introduire un fichier pointMotionUx dans le dossier 0 afin d'indiquer la vitesse de tous les bords du domaine de calcul. Toutes les mailles vont alors se déplacer de manière cohérente. Dans cet exemple, les mailles initialement à gauche du cône vont s'étirer lorsque ce dernier va se déplacer vers la droite, et les mailles à droite vont logiquement se contracter. Remarquons que la qualité du maillage est par conséquent modifiée au cours du calcul. Pour notre étude du mouvement du train, cet aspect sera un point primordial puisqu'il est évidemment impossible de réaliser des calculs avec un maillage de trop mauvaise qualité où certaines mailles auront un rapport hauteur sur longueur éloigné de 1.

 

Voici une animation décrivant l'évolution du maillage au cours du calcul, sans qu'on ait modifié les paramètres de la simulation :

Cas du tutoriel MovingCone

 

Le solveur en pratique


En pratique, l'utilisation de ce solveur nécessite d'effectuer quelques adaptations pour passer d'une configuration quelconque à une configuration fonctionnant avec ce solveur. En effet, il existe quelques contraintes auxquelles nous avons dû nous adapter pour faire fonctionner nos simulations.

Pour noter cas, où nous devions faire translater un objet sur une longue distance, nous avons utilisée une méthode qui consiste à mettre en place trois blocs. Le bloc central sera celui contenant l'objet à déplacer, et aura un maillage constant au cours du temps, au sens où le bloc bougera mais le maillage bougera avec sans se déformer. Il subira une simple translation.

A l'inverse, les deux blocs aux extrémités du domaine auront le rôle de tampon. Ils se compresseront ou s'étireront pour mettre en mouvement le bloc central.

Sur l'animation ci-dessus, on voit le bloc central se déplacer uniformément tandis que les deux blocs autour se compressent d'un côté et s'étirent de l'autre.

Dans toutes nos simulations futures, nous utiliserons cette configuration à trois blocs, avec le bloc central contenant le train, la sphère où tout autre géométrie.

Avec ces trois blocs, en plus des conditions initiales et limites sur la vitesse et la pression, il faut définir le déplacement des points du maillage. Ceci est réalisé dans le fichier pointMotionUx (mouvement suivant x) dans le dossier 0. Ce fichier est organisé de la même manière que les fichiers U ou p. Les parties left et right correspondant aux extrémités du tunnel sont fixées, avec fixedValue et uniform 0. Le train ainsi que le bloc central ont une vitesse de 40m/s définies avec fixedValue uniform 40. Enfin, les deux blocs "tampons" prennent la condition limite slip qui signifie que les points du maillage de ces blocs vont simplement glisser pour venir absorber le mouvement du bloc central.


Moving Reference Frame

Le deuxième concerne l'étude d'un mixeur en rotation dans un domaine fluide. Son intérêt est probablement de mélanger un fluide contenant plusieurs espèces afin d'homogénéiser leur répartition dans le fluide. Comme dans le premier exemple, le solveur qui gérera le maillage mobile est indiqué dans la fichier dynamicMeshDict et on peut alors remarquer qu'il s'appelle cette fois-ci solidBodyMotionFvMesh. Avec ce solveur, il faut indiquer certaines zones du domaine (zones MRF, pour Moving Reference Frame) de calcul et le déplacement de celles-ci sont indiquées dans le fichier MRFZones dans le dossier constant. La particularité supplémentaire par rapport au premier exemple est que, en plus de la méthode ALE pour gérer le fait que certains points de calculs sont mobiles, il est aussi nécessaire de faire des interpolations entre les zones MRF et les zones fixes. En effet, leur interface est fixe, mais les points de calculs de la zone MRF sur cette interface est mobile alors que ceux de la zone fixe sont immobiles. L'avantage de cette méthode est qu'il ne modifie pas la forme des mailles et donc si le maillage initial est de qualité correcte, il le sera tout au long du calcul.

 

Cas du tutoriel MixerVesselAMI2D

Pour deux raisons, nous préférerons donc partir sur le premier cas test. Le premier argument est géométrique. Il vient simplement du fait que nous souhaitons étudier la translation d'un train dans un domaine rectiligne, et ce cas du movingCone est relativement proche de ce que nous souhaitons finalement obtenir.

Le deuxième argument est davantage une anticipation de difficultés liées au solveur indiquant ces zones MRF. Outre le fait que nous aurions dans tous les cas du indiquer à certaines mailles de se contracter et de s'étirer (à hauteur du train), certaines zones auraient pu tout à fait être fixes (pour des hauteurs différentes de celles du train). Il aurait fallu interpoler les valeurs à l'interface entre ces deux zones du maillage. Cela n'étant pas indispensable, il ne nous apparaissait pas vraiment utile de se surcharger d'une contrainte supplémentaire avec cet outil d'interpolation, source probable d'erreurs.

Maillage d'une géométrie "complexe"

Le tutoriel dont on s'inspire traite initialement du mouvement d'un cône dans un cylindre. Ce type de géométrie peut être qualifié de "simple" puisqu'elle peut être décrite rapidement grâce au fichier blockMeshDict. En effet, avec seulement 5 blocs et 24 points, tout le domaine pouvait être défini.

Dans notre projet, il s'agit de décrire la géométrie d'un train, géométrie que l'on peut qualifier de "complexe". En effet, le cône est remplacé par un train, avec plusieurs wagons, ce qui engendre d'introduire beaucoup plus de points, et donc de blocs, avec le fichier blockMeshDict.

Nous avons également pu s'initier et utiliser le mailleur automatique snappyHexMesh qui permet de réaliser le maillage autour d'une géométrie qui, elle, est définie dans un fichier au format bien spécifique.

Nous allons ici vous présenter les deux méthodes possibles pour définir une géométrie plus complexe, que ce soit avec un fichier classique blockMeshDict ou avec l'utilisation du mailleur snappyHexMesh.

Avec BlockMesh

La première idée, a priori la plus simple, est de décrire manuellement les nombreux points composant la géométrie du train au sein même du fichier blockMeshDict. En effet, avec cette méthode, nous serons assurés d'avoir des résultats puisque l'on devra définir la vitesse des points des bords du domaine. Or, comme nous définissons nous même le domaine avec cette méthode, il suffira d'indiquer la vitesse adéquate à chacune des parois du domaine. Toutes les mailles à l'intérieur du domaine auront alors une vitesse en cohérence avec la vitesse des bords du domaine.

Voici pour rappel la géométrie du train qui nous a été fournie par Mr Tekam :

Dans un premier temps, nous allons négliger les congés qui sont censés améliorer l'aérodynamique globale du train. En effet, cela complique la création d'un maillage avec un blockMeshDict puisqu'il faut créer des arcs de cercle avec deux points de base et un troisième point par lequel passera l'arc de cercle. Suite à cette concession, voici à quoi devra ressembler notre géométrie :

Il est par conséquent nécessaire de décrire tout les points qui vont servir à décrire l'ensemble des blocs d'une telle géométrie.

Au lieu des 24 points définis pour définir le rectangle, il est nécessaire d'en introduire 32 supplémentaires pour un total de 56 points. De plus, au lieu des 5 blocs initialement présents, il faudra en définir ici 17. Voici ci-dessous un schéma permettant de définir tous ces blocs.

On peut remarquer que 4 blocs sont des triangles. C'est ici une particularité acceptée par OpenFOAM. En fait, au lieu de décrire une face d'un bloc par 4 points distincts, on peut décrire une face d'un bloc par 4 points dont 2 sont identiques, on aura alors un triangle. Lors de l'étape du maillage, nous aurons des mailles triangulaires dans l'un des sommets de ce triangle.

Concernant l'étape du maillage réalisée dans ce fichier blockMeshDict, il faut simplement faire attention à toujours définir un nombre de mailles identique pour les deux blocs voisins. Avec suffisamment de rigueur, cela est tout à fait faisable.

Voici ce que l'on obtient après avoir exécuté blockMesh dans le terminal et lancé paraFoam pour visualiser le domaine maillé :

Avant de commenter la géométrie finalement obtenue, on peut remarquer une particularité probablement gênante de ce maillage. Cela concerne les mailles triangulaires. Pour assurer un nombre minimum satisfaisant de mailles dans le domaine, il est en effet obligatoire d'imposer un nombre de mailles dans ces blocs triangulaires. Cela a pour conséquence d'avoir des mailles triangulaires extrêmement petites par rapport à ses voisines rectangulaires. Ceci est certainement source de divergence et rend plus contraignant le choix du pas de temps du calcul pour assurer sa convergence.

L'idée est alors d'éviter ces mailles triangulaires en modifiant légèrement la définition des blocs. Au lieu d'introduire des blocs triangulaires, nous pouvons définir 2 points supplémentaires pour chaque côté du train afin d'obtenir un bloc trapézoïdal. Plus ces deux points sont éloignés de la paroi du train, plus la qualité des mailles initialement triangulaires (qui sont évidemment devenues des trapèzes) sera acceptable. Cela se fait par contre au détriment de la qualité des mailles initialement rectangulaires des blocs voisins qui deviennent des trapèzes. Par conséquent, on induit quelques erreurs de calculs puisque ces mailles ne sont plus orthogonales.

Voici finalement ce qu'on peut obtenir en visualisant le maillage après avoir exécuté paraFoam.

Cette méthode permet donc d'obtenir de manière assez simple un maillage correct, mais reste valable uniquement pour des géométries très simples composées de carrés, de triangles ou d'autres polygones réguliers.

Avec SnappyHexMesh

SnappyHexMesh est un mailleur automatique intégré à OpenFOAM. Il consiste, à partir d'un maillage hexaédrique et d'une géométrie superposée dessus, à créer un nouveau maillage hexaédrique prenant en compte cette géométrie.

L'outil SnappyHexMesh

Pour utiliser, en tapant la commande snappyHexMesh -overwrite, il est nécessaire d'avoir préalablement créé une géométrie au format .stl ou .obj (description des points et surfaces de l'objet) ainsi qu'un maillage de fond.

Création du fichier .obj ...


SnappyHexMesh peut lire deux types de format pour définir l'objet que l'on souhaite introduire sur le maillage de fond. En effet, cet outil est capable de lire les fichiers .stl ou .obj. Dans notre étude, nous sommes passés par des fichiers de type .obj.

Pour cela, il faut utiliser des logiciels de CAO pour définir une surface fermée par des segments, des courbes ou des points. Le logiciel choisi se dénomme FreeCAD, qui est un logiciel disponible librement et gratuitement ici

Cette étape se résumait à ​définir les premiers points définissant grossièrement la géométrie. Ensuite, pour réaliser les congés à l'avant et l'arrière du train, de nomreux points ont été ajoutés afin de correspondre au mieux à la géométrie désirée. Enfin, il a fallu réaliser une extrusion de cette géométrie 2D afin de lui donner une épaisseur et d'avoir une surface 3D fermée.

Voici ce que l'on peut obtenir au final :

 

Création du premier maillage de "background" ...


Ceci se fait se façon assez classique avec un fichier blockMeshDict. On peut créer le maillage de fond avec la commande blockMesh.

 

Lancement de snappyHexMesh ...


Une fois ce fichier .obj obtenu et bien placé, et une fois le maillage de background créé, on peut lancer snappyHexMesh avec la commande snappyHexMesh -overwrite.

​Concrètement, lancer ce mailleur consiste à réaliser les étapes suivantes : 

Étape 1 : Intersection des mailles avec la géométrie

Repérer les intersections du maillage de fond avec la géométrie et affiner les mailles au niveau de ces intersections. Le niveau de raffinement de ces mailles dépend de ce que l'utilisateur a indiqué dans le fichier snappyHexMeshDict. Le niveau 0 correspond au maillage de fond, le niveau 1 correspond à des mailles deux fois plus petites et ainsi de suite...


Source : http://www.openfoam.org/docs/user/snappyHexMesh.php

Étape 2 : Suppression des mailles solides

Pour l'instant, tout le domaine reste maillé. Il faut pouvoir indiquer à snappyHexMesh quelle est la région qui correspond au solide (à l'intérieur ou à l'extérieur de la géométrie). Cela se fait dans le fichier blockMeshDict.

Voici ci-dessous une illustration du domaine obtenu :


Source : http://www.openfoam.org/docs/user/snappyHexMesh.php 

On voit alors que toutes les mailles sont carrées. Il serait donc utile de pouvoir éviter d'avoir ces mailles carrées au niveau de la surface de l'objet. SnappyHexMesh est capable de gérer cette problématique, comme on peut le voir lors de l'étape 3.

Étape 3 : Amélioration des mailles à la surface de l'objet

Par l'intermédiaire d'une légère déformation des mailles proches de l'objet, on peut réussir à parfaitement suivre les points de la surface de l'objet, comme on peut le voir sur l'image suivante.


Source : http://www.openfoam.org/docs/user/snappyHexMesh.php

                                                  

Cet outil, au premier abord difficile à adopter et à utiliser, reste toutefois extrêmement performant. Son principal inconvénient reste qu'il maille automatiquement le domaine. Par conséquent, l'utilisateur n'a pas un contrôle total sur toutes les mailles, notamment leur qualité.

 

Le fichier snappyHexMeshDict

Nous allons maintenant nous intéresser à la structure du fichier snappyHexMeshDict et détailler les principales parties qui composent ce fichier. Les paramètres modifiés seront également expliqués. Toutefois, les options qui n'ont pas été modifiées durant le projet ne seront pas explicitées. Pour plus d'informations concernant snappyHexMesh, le User Guide constitue une bonne base. Le rapport de Projet de Fin d'Etudes réalisé par la promotion 2011 de l'IPSA et intutilé "Présentation, essai et validation du logiciel Open-Source OpenFOAM" détaille également de manière approfondie chaque partie du fichier snappyHexMeshDict.


Le fichier snappyHexMeshDict est situé dans le dossier system/. Il commence comme tout autre fichier OpenFoam, par une en-tête contenant les informations sur le fichier.

Il s'ensuit alors trois lignes pour définir quelles options l'utilisateur souhaite activer. Ces options s'activent en mettant true ou false à la suite du mot clé.

La première option castellatedMesh intervient pour raffiner le maillage autour de l'objet considéré. L'option snap concerne la découpe des mailles et enfin, l'option addLayers désigne l'ajout ou non d'un couche de prisme le long du solide afin par exemple d'obtenir de meilleurs résultats dans la couche limite.

Une fois les options choisies, l'utilisateur doit alors désigner le nom du fichier .obj ou .stl à mailler. Il peut également changer la façon dont est désigné l'objet dans la suite du fichier avec la ligne name.

La suite du fichier s'organise globalement de la manière suivante :

La première partie concerne les paramètres de raffinement de l'option castellatedMesh. Parmi ces options, une est particulièrement intéressante puisqu'il s'agit de définir les régions de raffinement, qui permettent d'obtenir différents niveaux de raffinements du maillage. Il existe plusieurs manières de définir ces zones. Nous les avons définies par rapport à leur distance de la surface de l'objet. La syntaxe est la suivante :

mode distance ;

levels ((distance1 niveau1)(distance2 niveau2)(distance3 niveau3)...);

On définit ainsi les différentes zones de raffinement.

Les autres paramètres de raffinement contiennent divers paramètres comme par exemple le nombre de cellules pour passer d'un niveau de raffinement à un autre, le nombre de cellules maximum etc...

La seconde partie appelée snapControls contrôle la découpe du maillage autour de la géométrie. Elle contient des paramètres de type tolerance où nombre d'itérations maximum.

Enfin, la dernière partie concerne l'ajout de couche de prismes le long de la surface de l'objet. Elle contient de nombreux paramètres sur la qualité des cellules, l'épaisseur des mailles, le nombre de cellules de la couche de prismes etc...

Globalement, seule les paramètres de la première partie ont été modifié pour parvenir au résultat final. Voici par exemple un des maillages que nous avons pu obtenir avec snappyHexMesh :

On remarque bien les différents niveaux de raffinements autour de l'objet, définis par leur distance à ce dernier.