Simulation du mouvement d’un train avec le logiciel OpenFOAM - Mise en place du maillage mobile


Simulation du mouvement d’un train avec le logiciel OpenFOAM 

Mise en place du maillage mobile


 

 


 

 

 


Idriss Benkour, Nicolas Brett & Thibaut Moy

Encadré par Dominique Legendre & Annaïg Pedrono (INP - IMFT)

Contact industriel : Sydney Tekam (Techam-Ingénierie)


Ce sujet présente la mise en place d'une simulation du mouvement d'un train avec le logiciel libre OpenFOAM afin d'étudier les phénomènes physiques ayant lieu dans un espace confiné. Plus particulièrement, sera détaillé ici la mise en place du maillage mobile avec un tel code de calculs de CFD.

Tout le travail a été réalisé en collaboration avec l'entreprise Techam-Ingénierie et son fondateur Mr Sydney Tekam ainsi qu'avec un autre groupe d'étudiants se penchant davantage sur les modèles de turbulence pertinents à utiliser dans de telles géométries.

L'objectif est alors de pouvoir simuler de passage d'un train dans un tunnel en implémentant différents modèles de turbulence. Ce travail s'inscrit dans une étude des différentes contraintes que les parois du tunnel subissent dues à la présence du train.

Contexte général

Le sujet ...

La traversée d'un train dans un tunnel est un phénomène physique qui n'est pas parfaitement connu par le monde de la Mécanique des Fluides. Il rallie en effet plusieurs grandes problématiques de la CFD comme un domaine mobile, des écoulements probablement compressibles, un nombre de Reynolds important impliquant une modélisation turbulente du problème. Du point de vue de l'aérodynamique, là encore, la problématique est plutôt rarement étudiée vu qu'il s'agit d'un écoulement externe, autour du train, mais ne pouvant pas se rapprocher à l'aérodynamique classique vu que l'espace est confiné, ce qui rend tout point de comparaison à des phénomènes vus en aérodynamique en domaine ouvert impossible. Les distances mises en jeu constitue également un problème du point de vue du temps de calcul, ce qui contraint les industriels à adopter des solutions telles que des modèles 1D. Toutefois, la CFD reste un outil puissant quand il s'agit de vérifier et de valider des modèles.

D'une part, une connaissance plus précise des phénomènes mis en jeu dans ce type de mécanisme intéresse les constructeurs de trains dans l'idée d'améliorer les performances aérodynamiques de leur produit et le confort des passagers. D'autre part, cela intéresse aussi les constructeurs de tunnels qui souhaitent connaître les efforts que subissent les parois du tunnel dans le but de mieux décrire la fatigue mécanique éprouvée au cours du temps par ces structures. Cela leur permettrait de renforcer de manière plus efficace et plus pertinente les zones du tunnel sensibles au passage du train afin d'éviter tout écroulement du tunnel sans avoir à prendre trop de précautions grâce à une relative confiance en les résultats d'un calcul CFD.


Le projet ...

C'est dans ce cadre que l'entreprise Techam-Ingénierie, société de conseils en ingénierie, a été sollicitée par la société d'ingénierie et de conseils SETEC, spécialisée notamment dans le BTP, afin de réaliser une étude sur le mouvement d'un train dans un tunnel.


Source : http://www.lerm.fr/le-lerm-et-le-groupe-setec.html                                   Source : http://techam-ingenierie.fr/

Mr. Tekam a alors accepté cette mission et a décidé, en accord avec ses clients, d'utiliser le code de calculs libre OpenFOAM pour une première approche de ce problème. C'est dans ce cadre que Mr. Tekam s'est renseigné auprès des encadrants du BEI à l'ENSEEIHT afin de réaliser ces travaux de type "étude de faisabilité". Ce sujet nous paraissant intéressant, nous nous sommes tous portés volontaires pour participer à ce projet et nous avons eu le plaisir d'y être acceptés.


Les objectifs ...

Plus concrètement, deux groupes participent à ce sujet. L'un de ces groupes réalisera une étude des modèles de turbulence dans le cas d'aérodynamique en espace confiné (que vous trouverez ici). Un maillage précis d'un train dans un tunnel leur a été fourni et ils mettront en valeur les défauts et les atouts de quelques modèles judicieusement choisis.

Le deuxième groupe, nous concernant, essaiera quant à lui de mettre en place les procédés de maillage mobile proposés par OpenFOAM et de résoudre les problèmes géométriques apparaissant lorsque le train parcours le domaine d'un bout à l'autre.

Une phase finale aura ensuite lieu pour un travail commun entre ces deux groupes qui consistera à, une fois le maillage mobile maîtrisé et les modèles de turbulence mieux connus dans ce type de problèmes physiques, étudier les modèles de turbulence avec le maillage mobile puisque c'est finalement l'objectif final de ce sujet de BEI.


Les hypothèses ...

Ce sujet étant fortement axé sur un travail de CFD, nous devons faire une première hypothèse qui, il est vrai, n'est pas forcément respectée. Nous nous sommes, en accord avec Mr Tekam, tout de suite tourné vers un écoulement de type incompressible. En effet, il nous paraissait dangereux, dans le temps imparti, de prendre compte des problèmes numériques rencontrés avec les écoulements compressibles comme la capture des chocs ou le traitement aux bords du domaine pour éviter toute réflexion d'onde acoustique ou isentropique. De plus, le cahier des charges stipule que la vitesse du train doit être de l'ordre de $40 m/s$,  ce qui correspond à un nombre de Mach de $Ma=0.11$ et la limite pour les écoulements incompressibles est autour de $Ma=0.3$, ce qui valide toutefois cette hypothèse.

L'étude a également été menée en deux dimensions, pour des raisons de temps de calcul principalement, étant donné les moyens de calcul dont nous disposons pour ce projet. En effet, par la suite, le nombre de mailles utilisé en 2D est de l'ordre de 100 000. Pour avoir un maillage correct dans la troisième direction, il aurait fallu au moins une vingtaine de mailles aux endroits les moins raffinés, ce qui aurait impliqué plusieurs millions de mailles, et ce n'était clairement pas envisageable pour ce projet.


Plan de l'étude

L'étude s'est alors organisée de la façon suivante :

    1. Problème physique et cahier des charges

    2. OpenFoam

    3. Validation des outils utilisés

    4. Mise en place de la simulation numérique

    5. Ajout de la turbulence

    6. Résultats : analyse et comparaisons

    7. Perspectives

    8. Remerciements

    9. Références

 

Problème physique & Cahier des Charges

Le problème à étudier dans ce sujet de BEI est le parcours d'un train dans un tunnel.

L'objectif donné par les responsables industriels est principalement de mettre en place une méthodologie permettant de réaliser cette étude du mouvement du train dans le tunnel. Les résultats, en soi, n'auront que peu d'importance. Nous devrons bien sûr exploiter quelques résultats afin de vérifier la cohérence des simulations obtenues, mais ce que souhaite avant tout Mr Tekam est d'avoir une méthodologie permettant cette analyse. Il se chargera lui-même de fixer les valeurs des différents paramètres du problème (géométrie, modèles de turbulence, pas de temps, modèles de discrétisation, etc ...).

Néanmoins, pour rester au plus proche de la réalité, avec Mr Tekam, nous nous sommes fixés une certaine géométrie de train, une longueur de tunnel, ainsi qu'une vitesse de déplacement du train.


La géométrie du train et du tunnel ...

Mr Tekam nous a fourni une géométrie du train par l'intermédiaire d'un maillage Fluent déjà utilisé pour des études similaires.

Les dimensions du train, qui fait cinq mètres de hauteur, ainsi que celles du tunnel sont données :


Le déplacement du train ...

Concernant la vitesse de déplacement du train, elle est fixée à 40 m/s, ce qui est tout à fait cohérent avec ce qu'il se passe en réalité. Cela justifie évidemment la prise en compte d'un modèle turbulent. En effet, pour de l'air, et avec un train de cette hauteur, nous avons un nombre de Reynolds de l'ordre de :

 

$Re=\frac { U \cdot L}{\nu} = \frac {40  \cdot 5}{10^{-5}}=20 \cdot 10^6$

 

Les écoulements étudiés seront donc turbulents. Cependant, dans notre démarche, nous allons dans un premier temps réaliser des simulations en régime laminaire, puis dans un deuxième temps nous augmenterons la vitesse en incorporant un modèle de turbulence. Ce modèle, ainsi que la valeur de ses paramètres, sera celui que le groupe 24 nous aura conseillé car ils auront préalablement effectué des tests de différents modèles, en maillage fixe, sur la géométrie du train.


Les paramètres de la simulation ...

Ensuite, l'idée sera de simuler ce problème pour une durée suffisamment grande pour que l'influence de la condition initiale soit nulle. Concrètement, d'après nos responsables industriels, si l'on simule ce problème pour une durée de 5 à 10 secondes, tous les phénomènes physiques auront le temps de s'établir. Avec une telle vitesse et ce train donné, cela signifie que l'on doit faire déplacer le train sur environ 5 à 10 fois sa taille. On peut alors en déduire une longueur du domaine représentant le tunnel, qui sera alors fixée à 470m.

Ces paramètres pourront bien sûr être modifiés à la guise de l'utilisateur. Ces données servent uniquement à prendre connaissance de tous les ordres de grandeurs d'un cas concret classique.

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.

Mise en place de la simulation numérique

Une fois le choix du solveur effectué, et l'étape de maillage bien comprise, il faut maintenant parvenir à raccorder tous les outils. C'est généralement la difficulté avec OpenFOAM, qui est une excellente boîte à outils regroupant énormément d'utilitaires, de solveurs etc... mais il est parfois très difficile de faire concorder tous ces outils.

Pour cela, la partie suivante traite de la stratégie que nous avons mis en place afin de parvenir à lancer une simulation complète du train se déplaçant sur une grande distance. Les problèmes de déformation de mailles posés par le solveur doivent être résolus sinon ils mèneront à la divergence du calcul.

Stratégie utilisée

Le cas initial ...

Lorsqu'on reprend le tutoriel choisi, à savoir l'étude d'un cône se déplaçant dans un domaine fermé, nous avons pu voir qu'il y avait trois sections distinctes.

Le bloc central, contenant le cône ainsi que les mailles situées entre lui et la paroi supérieure du cône, se translatait de la gauche vers la droite. Il est important de noter que l'on pourra incorporer la géométrie du train dans ce bloc central.
De part et d'autre de cette partie centrale du domaine, chacun des blocs se dilate (celui de gauche) ou se contracte (celui de droite). Cela permet bien évidemment de conserver un domaine aux dimensions fixes, mais un problème majeur important se met en valeur : celui de la déformation des mailles.

Pour mettre en place notre raisonnement, rappelons que nous conserverons toujours cette idée de séparer le domaine en 3 blocs. le bloc central contiendra la géométrie complexe maillée avec snappyHexMesh :


Le problème rencontré ...

 

En effet, à l'instant initial, le maillage semble un peu trop contracté à gauche et un peu trop étiré à droite. Cependant, au fur et à mesure que le calcul avance, le bloc central se translate vers la droite. Ainsi, le bloc de gauche, qui avait initialement des mailles compactes, voit ses volumes de contrôle s'étirer de manière conséquente. A l'inverse, les mailles du bloc de gauche de resserrent de façon non négligeable.

On peut illustrer cette déformation avec ces deux images. la première est l'instant initial du cas de base fourni par OpenFOAM. La seconde correspond à l'instant final.

Globalement, le maillage nous paraît être de mauvaise qualité si l'on discute en terme de ration "hauteur par rapport à longueur" lors des instants initial et final. Concernant notre étude, ce domaine sera assimilé à un tunnel. Le domaine sera donc encore plus long que celui du tutoriel. Les blocs de gauche et de droite seront encore plus déformés. Dès lors, il nous paraissait inconcevable d'utiliser cette technique sans rien modifier tant la pauvreté du maillage rendrait notre confiance en ces résultats assez faible.

 


La solution envisagée ...

L'idée est de résoudre cette problématique tout en conservant ce solveur qui sépare le domaine de calcul en 3 blocs distincts.

Il faudrait alors, à l'instant initial, mailler le domaine de manière à avoir un maillage à la qualité satisfaisante. On exécute OpenFOAM et le solveur pimpleDyMFoam afin de réaliser les premières itérations. Le maillage est alors légèrement déformé au cours de ces itérations.

Toutefois, ce premier calcul n'ira pas au terme de la simulation. En effet, nous considérons que au bout d'un certain temps $t_0$, le maillage aura été trop déformé, on s'arrête à cet instant là. Il faut alors retenir la solution de chaque champ des différentes variables calculées (U, p, k, $\epsilon$, etc ...). Une fois retenue cette solution à $t=t_0$, il faut alors remailler le domaine de manière à avoir une qualité de maillage de nouveau satisfaisante. A partir de ce nouveau maillage, on peut relancer le solveur pimpleDyMFoam. On réalise encore quelques itérations jusqu'à ce qu'on arrive à une date $t=t_1$ où l'on considérera que la qualité du maillage n'est plus assez satisfaisante. On retiendra les différents champs des grandeurs calculées pour appliquer cette solution à un nouveau troisième maillage pour relancer une troisième sous-simulation. On procède ainsi jusqu'à ce qu'on ait atteint le temps final de simulation désiré.

Concrètement, afin de choisir ces dates $t_i$ où il faudra remailler et donc changer de sous-simulation, nous diviserons la simulation complète en N sous simulations. Ainsi, pour une durée de simulation de T secondes, nous remaillerons le domaine toutes les $\frac {T}{N}$ secondes. Il faut ensuite choisir ce paramètre N de manière avoir suffisamment de sous-simulations pour toujours conserver un maillage ayant une qualité recevable.

Enfin, nous verrons dans la paragraphe suivant que OpenFOAM a un outil extrêmement pratique permettant d'appliquer une condition initiale particulière à un calcul. Cette condition initiale peut tout à fait venir d'un calcul déjà effectué. Si le maillage est différent, OpenFOAM pourra interpoler les champs de toutes les grandeurs de l'ancien maillage sur le nouveau maillage.


Par la figure ci-dessous, nous pouvons donc schématiquement représenter notre vision du problème et notre organisation pour traiter le parcours complet du train dans un tunnel. La page suivante permettra alors de retranscrire cet idéogramme avec un script en bash.

 

Automatisation du calcul

Nous nous sommes rendus compte dans le paragraphe précédent qu'il fallait diviser la simulation complète en plusieurs sous-simulations. Entre chacune de ces sous-simulations, il faut remailler le domaine et interpoler les différents champs calculés grâce à l'ancien maillage sur le nouveau maillage.

Vu la longueur du tunnel par rapport à sa hauteur, nous savons qu'il va falloir réaliser un nombre relativement conséquent de sous-simulations, de l'ordre de 10 à 20 a priori. Il serait donc extrêmement utile de créer un script bash permettant de tout automatiser. L'objectif serait de uniquement créer "manuellement" un dossier type OpenFoam classique comme si la simulation ne se faisait qu'en une seule fois, avec tous les paramètres désirés (maillage initial, pas de temps, temps de simulation, schémas, etc ...). Puis à partir d'un script bash, les étapes de remaillage et d'interpolation séparant les sous-simulations seraient automatisées.

 

Vous trouverez ci-après des explications sur l'utilisation d'OpenFOAM pour mener à bien le passage d'une sous-simulation à la suivante ainsi qu'un autre lien menant à la page détaillant le script avec plus de précisions.

Gérer avec OpenFOAM le passage d'une sous-simulation à la suivante

 

Avant d'expliquer plus en détails la structure et les commandes introduites dans le script, revenons à la façon dont OpenFOAM est capable d'appliquer une solution (connue, calculée avec un maillage donné) sur un autre maillage en interpolant toutes les valeurs sur les nouveaux points de calcul.

Cette procédure est très bien détaillée dans un tutoriel disponible sur le site d'OpenFOAM (ici, section 3.1.5.3), mais expliquons comment cela fonctionne.

Pour cela, il faut disposer d'une ancienne simulation déjà réalisée, contenant tous les résultats dont la solution finale correspondant à la date $t=t_0$. Dans le fichier controlDict de cette ancienne simulation, l'utilisateur avait indiqué que le temps final de la simulation était bien $t_0$.

Il faut également créer un deuxième dossier, comme s'il s'agissait d'une nouvelle simulation. Ce dossier contient alors classiquement les dossiers constant, system, et également un dossier de condition initiale qui doit être dénommé $t_0$ au lieu de 0.

Dans ce deuxième dossier, il faut également introduire un fichier mapFieldsDict. Ce fichier indique quelles zones de l'ancien maillage doivent être interpolées sur quelles zones du nouveau maillage. Cela se fait par les noms de patch du domaine. Voici, dans notre cas, notre fichier mapFieldsDict :

La section cuttingPatches est utile dans le cas où les domaines des deux simulations sont différents. Les zones absentes sur le nouveau maillage, par le nom de certains patches, sont indiquées dans ce paragraphe ; ce qui n'est pas le cas pour notre étude, les deux domaines sont parfaitement identiques.

Une fois ce fichier créé dans le dossier system, il faut d'abord générer le maillage sur la nouvelle simulation, avec la commande blockMesh suivi de snappyHexMesh.
Puis on indique que la condition initiale de cette nouvelle sous-simulation est la solution finale de l'ancienne grâce à la commande mapFields ../simulation_précédente.
Enfin, on lance l'exécutable avec la commande classique pimpleDyMFoam si on veut un calcul en séquentiel.

 

Le script en détails

Le script peut être décomposé en plusieurs parties, comme le montre schématiquement la figure présentée dans le chapitre précédent qui illustre la stratégie employée.

Avant toute chose, pendant la première étape, l'exécutable a besoin de deux paramètres que sont le nombre de sous-simulations ainsi que le nombre de processeurs pour réaliser le calcul. En effet, pour gagner du temps, le calcul se fait en parallèle. Cette entête du script contient également quelques informations supplémentaires. Il est en effet nécessaire de repérer où le script est exécuté, d'indiquer le nom du dossier initial qui contient tous les paramètres choisis par l'utilisateur, et aussi d'importer des formats de nombre et autres fonctions spécifiques.
Cette partie prépare également le dossier Simulations où seront répertoriés tous les dossiers de sous-simulation.


La deuxième étape est de définir tous les paramètres de simulation. Parmi ces paramètres il a fallu :

  • Lire le premier argument du script afin de connaître le nombre de sous-simulations
     
  • Lire le deuxième paramètre qui indique le nombre de processeurs désiré. Il faut ensuite modifier de fichier OpenFOAM decomposePartDict qui gère cette parallélisation du calcul. Cette modification se fait grâce à la commande sed -ri, qui va soit rechercher la ligne contenant le mot clé indiqué après 's/, soit le numéro de ligne indiqué par exemple par " '21 " pour la ligne 21. Après avoir repéré cette ligne, cela supprimer la ligne en question et on remplace toute la ligne par ce qui est indiqué dans la suite de la commande, avec le vocabulaire spécifique aux commandes bash relatives à l'écriture dans des fichiers.
     
  • Repérer, dans le dossier de base créé par l'utilisateur, les différents paramètres de la simulation comme le temps final, la vitesse du train, le pas de temps, la fréquence de sauvegarde des résultats, le nombre de mailles initiales à gauche et à droite du train, etc ...
    Pour chacun d'eux, il a fallu utiliser astucieusement la commande echo pour faire afficher une ligne isolée par la combinaison de head et de tail. Par exemple, pour isoler la ligne 26, on affiche les 26 premières lignes du fichier qui contient l'information voulue. Puis sur ce résultat, avec un pipe et la commande tail, on ne retient que la dernière ligne.
    Une fois cette ligne isolée, il faut la diviser en plusieurs champs grâce à la commande cut qui va séparer la ligne en plusieurs parties en indiquant le séparateur.
    Une fois le paramètre enregistré dans le script, il est parfois pratique d'en tirer d'autres informations, comme la longueur du domaine par exemple. Réaliser des opérations en bash peut se faire avec la commande echo "opération" |  bc.
     
  • Une dernière section permet de calculer le nombre de mailles idéal qu'il faudra retirer à gauche du train, et ajouter à droite. On calcule ce nombre optimal en imposant le fait de conserver la même taille caractéristique de maille entre le maillage initial et le maillage final. Cela est faisable vu que l'on peut connaître la position du bloc central à l'instant final. En effet, on a à notre disposition la vitesse du train, sa position initiale et le temps de simulation.
    D'ailleurs, ayant ces informations, on peut d'ores et déjà savoir si ce bloc "train" va dépasser du domaine. Si c'est le cas, la commande exit permet de stopper le calcul.

 

 


La troisième étape consiste en la création des différents dossiers de sous-simulations. En effet, comme nous l'avons vu avec la commande mapFields d'OpenFOAM, pour chaque sous simulation, il faut un dossier distinct. Chacun de ces dossiers est à l'origine une copie du dossier de base créé par l'utilisateur. On modifie ensuite le fichier controlDict en changeant le temps initial et temps final. Cela se fait avec des commandes déjà expliquées plus haut : le calcul des temps se fait avec echo "Opération" | bc, puis on modifie la ligne désirée du fichier controlDict avec sed -ri.

 


La quatrième étape est la boucle à proprement parler, que l'on a déjà pu approcher dans la page "Stratégie utilisée".

  • Avant de lancer la boucle pour enchaîner les simulations, pour une question d'ordre d'enchaînement des commandes OpenFOAM, il faut lancer à part la création du maillage de la première sous-simulation avec les commandes runApplication blockMesh, snappyHexMesh et extrudeMesh.
     
  • Cette génération de maillage crée des fichiers non désirés. Par conséquent, en calculant les temps initial et final de la sous-simulation en cours de traitement, on peut se déplacer dans le dossier adéquat pour supprimer ces fichiers.
     
  • Puis on peut lancer le calcul en parallèle avec les commandes runApplication decomposePar, `getApplication et reconstructPar. S'il s'agissait de la dernière simulation, on sort de la boucle grâce à un test et la commande break si ce test s'est avéré positif.
     
  • Une fois le calcul terminé, l'idée est de préparer le remaillage de la sous-simulation suivante. En effet, nous savons déjà combien de mailles il faudra supprimer et ajouter aux blocs, respectivement de gauche et de droite. Par contre, nous ne savons pas où il faut positionner ces blocs, tout comme la géométrie du train décrite dans le fichier .obj.
    Il faut donc être capable de lire la position de l'interface séparant le bloc de gauche et le bloc central. Pour cela, nous avons ajouté une section au solveur afin qu'il indique en sortie, dans le fichier log.pimpleDyMFoamModified la position de cet edge. Cette démarche est expliquée à la page "Modifications du solveur".
    Ainsi, en ouvrant ce fichier texte log.pimpleDyMFoamModified avec cat, puis en recherchant la ligne contenant le mot clef "Left" avec grep, on peut repérer la position de cet edge séparant le bloc central du bloc de gauche à la fin de la simulation qui vient de se terminer.

 

  • Une fois cette position connue, on peut concrètement rentrer dans le dossier de la sous-simulation suivante. Il faut modifier le maillage en spécifiant la nouvelle position des edges en retouchant les lignes correspondantes du fichier blockMeshDict. Il faut également modifier le nombre de mailles par bloc. Tout ceci est réalisé grâce à la commande sed -ri que l'on a déjà expliquée un peu plus haut. De plus, il est nécessaire de translater le fichier .obj qui permet de décrire le train. Cela peut se faire avec la commande surfaceTransformPoints -translate ' "vecteur" ' train.obj.
     
  • Enfin, une fois que les fichiers définissant le maillage de la sous-simulation sont correctement définis, on peut lancer le processus de maillage avec les commandes runApplication blockMesh, snappyHexMesh,  extrudeMesh et mapFields pour interpoler les valeurs des variables sur ce nouveau maillage.
    Suite à cette interpolation de maillage, il faut juste renommer le fichier pointMotionUx qui définit la vitesse des points de chaque noeud du maillage.
     
  • On peut passer à la sous-simulation suivante, et ainsi de suite ...


La cinquième et dernière étape permet de rassembler tous les résultats, pour l'instant séparés dans des dossiers de simulations différents, dans un seul et même dossier. Cela est bien plus pratique pour le post-traitement, comme par exemple avoir toute la simulation en une seule animation (une seule ouverture de paraview).
Ce dossier contenant l'ensemble des résultats est appelé data. On le crée avec la commande mkdir. Puis on parcoure tous les dossiers de sous-simulation et on exécute la commande foamToVTK. Cela crée un dossier VTK rassemblant les résultats de la sous-simulation sous le format .vtk. On peut alors les visualiser avec la commande paraview, et non paraFoam.
Enfin, tous ces fichiers .vtk sont déplacés dans le dossier data. Il est juste nécessaire d'incrémenter, à chaque pas de sauvegarde successif, l'indice à la fin du nom du fichier. Connaissant l'instant auquel correspond chacun des fichiers, on peut attribuer correctement un indice permettant de ranger dans l'ordre croissant tous ces nombreux fichiers.

Modifications du solveur

Préparation pour la compilation du nouveau solveur


Afin de mener à bien la simulation finale, celle-ci doit être décomposée en plusieurs sous-simulations. Le fait de diviser ainsi la simulation ajoute un certain nombre de contraintes qui doivent être résolues. Une de ces contraintes est qu'il faut repérer entre chaque simulation la position finale du train de la simulation $n-1$ pour pouvoir communiquer à la simulation $n$. En effet, c'est cette information qui est nécessaire pour pouvoir remailler la géométrie translatée. Cependant, dans le solveur utilisé pimpleDyMFoam, cette information n'est présente nulle part de manière explicite. Il faut donc l'implémenter dans le solveur. Cette étape nécessite de bien comprendre comment s'organise un solveur, et d'être capable de coder quelques lignes en C++, même si de nombreuses aides sont disponibles ce qui facilite grandement la tâche.

Il faut tout d'abord créer une copie du solveur que l'on veut modifier, pimpleDyMFoam. Pour cela, on crée un dossier qui contiendra les solveurs modifiés par l'utilisateur, et on vient copier le solveur original dans ce dossier.

On possède alors une copie du solveur placé dans le dossier $WM_PROJECT_USER_DIR et renommée pimpleDyMFoamModified. Plusieurs fichiers sont présents dans le dossier mais celui qui nous intéresse ici est le fichier principal du solveur qui est pour le moment pimpleDyMFoam.C, que l'on renomme en pimpleDyMFoamModified.C, et on supprime les dépendances et autres fichiers binaires créés lors de la compilation du solveur d'origine.

Une dernière modification est nécessaire avant de pouvoir passer à l'implémentation à proprement parler. Il faut aller modifier le fichier Make/files dans le dossier du solveur qui indique quels fichiers seront à compiler. Le dossier Make est l'équivalent d'un Makefile. On remplace les lignes du fichier par les suivantes :

On indique ici que le fichier à compiler sera pimpleDyMFoamModified.C et que l'exécutable créé se nommera pimpleDyMFoamModified et sera placé dans le bin de l'utilisateur, c'est-à-dire là où sont rangés tous les exécutables créés par l'utilisateur.

 

Modification du solveur


Une fois toute cette préparation préliminaire effectuée, on peut alors passer à l'ajout du calcul de la position du train à la fin d'une simulation. Pour cela, on se place à la fin du fichier principal du solveur et on ajoute les lignes suivantes :

La méthode utilisée est ici assez simple. Notre domaine étant divisé en trois blocs distincts, on vient chercher la position du point le plus à droite du patch farFieldLeft (qui correspond à la condition limite paroi du tunnel sur le bloc de gauche) à la fin de la simulation. Connaissant sa position initiale et sa position finale, on peut alors calculer le déplacement du train.
Les lignes codées ci-dessus reprennent exactement cette démarche. La ligne 113 permet de trouver le patch que l'on cherche. On vient alors récupérer tous les points de cette partie du maillage. Pour cela, on trouve d'abord la partie du maillage correspondant au patch trouvé, à la ligne 115 puis les points correspondants. On déclare ensuite un nombre de type float nommé "max" que l'on initialise à -10000. Ce nombre va en fait nous permettre de trouver la position extrême du patch, en faisant une boucle sur tous les points identifiés comme appartenant au patch. Dès qu'un nombre est plus grand que "max", ce dernier prend sa valeur. Finalement, on affiche la position du segment de gauche du bloc du milieu, qui vaut "max". On est alors en mesure de connaître le déplacement du train durant une simulation.

 

Compilation et test du solveur


Il reste enfin à compiler ce nouveau solveur. Pour cela, on se place dans le dossier du solveur, et on exécute la commande wmake, ce qui aura pour résultat de créer l'exécutable du solveur, et puisque $FOAM_USER_APPBIN est ajouté au PATH, la commande pimpleDyMFoamModified sera disponible depuis le terminal. On peut alors tester notre nouveau solveur sur une simulation simple. Voici ce qui s'affiche à la fin d'une simulation :

Pour récupérer cette valeur, il suffit d'écrire les logs de la simulation dans un fichier, et de venir récupérer cette valeur avec les commandes grep, cut, etc ...

Une fois cette étape terminée, il faut notamment déterminer une procédure pour lancer les simulations et l'intégrer au sein d'un script en bash, qui permet d'automatiser ceci.

Validation des outils utilisés

Validation avec le cas du déplacement d'une sphère dans une conduite


Afin de valider les outils utilisés, que sont le solveur pimpleDyMFoam ainsi que le mailleur snappyHexMesh, nous allons comparer les résultats obtenus avec notre méthode par rapport à un écoulement type dont les résultats sont connus.

Pour cela, nous avons recherché dans la littérature des expériences, numériques ou non, qui proposaient des résultats et où l'utilisation d'un maillage mobile dans le cas d'écoulements confinés était appropriée. C'est dans ce cadre là que nous avons trouvé un document fort intéressant traitant du mouvement d'une sphère dans une conduite cylindrique.

Ce document, intitulé "Couplage fluide-solide appliqué à l'étude de mouvement d'une sphère libre dans un tube vertical", est le rapport de thèse effectuée par Mr Thibaut Deloze à l'Université de Strasbourg. Nous nous sommes notamment référés au chapitre 7 : Etude d'un écoulement autour d'une sphère en translation uniforme le long d'un tube. Il est disponible en téléchargement libre sur le lien suivant : scd-theses.u-strasbg.fr/2203/

L'aspect intéressant de cette étude est que la vitesse de la sphère est fixée par l'utilisateur, ce qui est proche de notre étude finale qui est l'étude du déplacement du train à une vitesse imposée.

Nous allons ici comparer les résultats obtenus avec notre code par rapport à ceux de cette thèse dans le cadre de l'influence de la distance de la sphère par rapport à la paroi ; en effet, la sphère n'est forcément au centre de la conduite. Remarquons également que le maillage utilisé dans cette thèse contient 2.6 millions de mailles pour une étude en 3 dimensions.

Source : scd-theses.u-strasbg.fr/2203/

Les résultats de cette thèse donnent la valeur des coefficients de traînée et de portance en fonction de la position de la sphère dans la conduite, et ce pour 3 nombre de Reynolds différents : 50, 200 et 250. Le paramètre adimensionnel permettant de rendre compte de cette distance à la paroi est le paramètre L/d sachant que D/d est fixé à 3.3. Ce paramètre varie de 1 à 1.6.

 

 

Implémentation du problème physique au sein du code

Pour lancer ces simulations, nous allons ainsi utiliser le maillage mobile avec le solveur pimpleDyMFoam. Afin de prendre compte de la sphère solide qui se déplace, nous avons utilisé le mailleur automatique snappyHexMesh en ayant préalablement défini la conduite circulaire.

En revanche, pour des raisons de ressources de calculs, nous devrons nous restreindre à un cas bidimensionnel. En effet, réaliser une étude paramétrique avec un maillage mobile de 2.6 millions de mailles nous paraît irréalisable dans le temps imparti pour ce BEI. Nous devons donc bien nous rendre compte que le cas étudié n'est pas vraiment le même et ceci peut être source de non concordance entre nos résultats et ceux de la thèse.


Définition de la sphère ...

Pour prendre en compte la sphère, nous allons utiliser le mailleur snappyHexMesh. Ainsi, il nous faut un fichier .obj. Ceci peut être obtenu, en 3 dimensions, avec le logiciel Paraview. Il faut en effet, après avoir lancé le logiciel, cliquer sur "Source" puis "Sphere". Une sphère s'affiche alors à l'écran. On exporte ensuite cette sphere au format x3d qui peut être finalement converti en .obj.

Cette sphère fait donc, comme dans la thèse, 1 mètre de diamètre et aura une vitesse imposée de 1m/s, et est placé selon les cas plus ou moins proche de la paroi du domaine.

Voici finalement le maillage, autour de la sphère, que l'on obtient :


Définition du domaine ...

La stratégie employée est toujours celle décrite précédemment avec une séparation du domaine en 3 blocs. Le bloc central, en translation, contient la sphère et et les deux autres se contractent ou s'étirent.

Voici donc schématiquement comment ce domaine peut être représenté :

Au final, le maillage complet comporte 28 000 cellules. Nous aurons donc une longueur caractéristique de maille du même ordre de grandeur que le maillage utilisé pour la thèse sachant que ce dernier contient 2.6 millions de mailles mais réalise un calcul en trois dimensions.

Pour ce cas de la sphère, nous avons choisi de diviser la simulation en 10 sous-simulations afin de conserver un maillage de qualité correcte. Par conséquent, voici le type de simulation que l'on peut obtenir :


Obtenir le coefficient de traînée avec OpenFOAM ...

Ce paragraphe traite de l'évolution du coefficient de traînée Cx de la sphère dans l'écoulement. Il faut ainsi expliquer comment nous avons pu avoir accès à ce coefficient. Ce dernier est issu de la force de traînée qui est présente suite aux effets visqueux de l'écoulement ainsi que les efforts de pression de l'écoulement sur la surface de la sphère. On peut relier Cx à cette force par la relation suivante :

$ F_x=\frac {1} {2} \cdot \rho \cdot S \cdot C_x \cdot V² $

Grâce aux données de la vitesse et de pression autour de la sphère, il aurait été tout à fait possible de remonter à la force de traînée en réalisant une intégrale sur la surface de cette sphère.

Cependant, directement dans le fichier controlDict et par OpenFOAM, il est possible d'importer la librairie libforces.so qui permet ensuite d'obtenir comme résultat un coefficient de traînée. L'utilisateur doit indiquer une vitesse de référence (la vitesse de la sphère dans notre cas), une longueur de référence (le diamètre), une densité de référence ainsi qu'une surface de référence. Cette surface de référence, dans notre cas 2D, est calculée à partir du diamètre que l'on multiplie par l'épaisseur de la maille dans la troisième direction. Cette surface correspond en effet à la surface apparente de la sphère, vu par l'écoulement.

Voici par conséquent ce qu'il a fallu ajouter à la fin du fichier controlDict :

 

En ayant inséré ce paragraphe dans le fichier controlDict avant de lancer les calculs, OpenFOAM créera alors un dossier intitulé forceCoeffs où sera présent un fichier comportant le coefficient de traînée.

 


Les paramètres numériques ...

 

Enfin, nous avons utilisés les paramètres numériques suivants pour la résolution du problème sous OpenFOAM :

Solveur Discrétisation temporelle Discrétisation spatiale Pas de temps (s) Temps Final (s)
PimpleDyMFoam

Euler implicite,
Ordre 2

Gauss, Ordre 2 10-3 10

 

Comparaison des Résultats : influence de la distance à la paroi pour différents Reynolds

Etude de Convergence en maillage ...

Nous avons précédemment fait remarquer que le maillage obtenu était constitué de 28 000 cellules pour un calcul 2D. Il faut quand même pouvoir justifier un tel maillage pour le domaine décrit lors de la section précédente.

Pour cela, nous nous fixons à un nombre de Reynolds donné, égal à 200, et pour une distance à la paroi de manière à avoir un rapport $ \frac {L}{d}$ égal à 1. Pour différents maillages, chacun associé à un nombre de maille, nous avons utilisé OpenFOAM pour calculer le coefficient de traînée.

Ainsi, on peut tirer la valeur de Cx pour les 5 maillages testés, sachant que le résultat de la thèse, pour cette configuration, est $C_x = 1.01$

 

Nombre de mailles 3200 7200 12800 28000 64000
Cx 2.0678 2.0674 2.0598 2.055 2.0521

 

Ce tableau, sous forme de courbe, est représenté ci-dessous :

On peut alors remarquer que, quelque soit le maillage utilisé, la valeur du coefficient de traînée est constante et ainsi l'écart par rapport au résultat de la thèse reste identique peu importe la précision du maillage.

Nous avons ainsi choisi de conserver le maillage avec 28000 volumes de contrôles puisque ce dernier, avec le pas de temps choisi de 10-3 s, nécessite environ 1H de calcul. Un tel temps de calcul est convenable pour plusieurs raisons. Il assure un minimum de précision dans les résultats, avec un précision équivalente à celle de la thèse et permet également d'enchaîner les simulations de l'étude paramétrique dans de bonnes conditions.


Comparaison des résultats avec la thèse ...

Après avoir justifié ce choix de maillage, nous pouvons maintenant réaliser l'étude de l'influence du nombre de Reynolds et de la distance paroi-sphère. Notons simplement que nous faisons varier le Reynolds par l'intermédiaire de la viscosité du fluide, la vitesse de déplacement du solide étant toujours fixée à 1 m/s.

Ainsi, avec un post-traitement via Matlab, voilà la courbe que nous obtenons :

Ces courbes mettent en valeur un décalage certain entre les deux sources de résultats. Effectivement, les calculs 2D sous OpenFOAM ont tendance à surestimer ce coefficient de traînée par rapport aux résultats 3D de la thèse. Peu importe le Reynolds et la distance à la paroi, il faudrait retirer 1 à nos résultats pour retrouver environ ceux de la thèse.

La différence entre 2D et 3D peut sans doute être source de ce décalage. En effet, le problème 2D n'est pas parfaitement équivalent au cas 3D. Des phénomènes fluides ayant lieu lorsqu'on considère la 3ème direction ne peuvent être retranscrit en 2D et peuvent être à l'origine d'une augmentation de la traînée dans un cas bidimensionnel.

Malgré tout, les résultats semblent tout à fait corrects dans la mesure où les tendances sont respectées et les ordres de grandeur équivalents. Plus le nombre de Reynolds augmente, plus le coefficient de traînée diminue et l'éloignement de la sphère a une légère tendance à faire diminuer ce coefficient.

Ajout de la turbulence

Pour obtenir une simulation réaliste, il est évidemment nécessaire de prendre en compte la turbulence dans le calcul, puisque le $Re$ est de $20 \cdot 10⁶$, ce qui correspond à un régime pleinement turbulent.

Nous nous sommes entièrement basés sur les conclusions du travail du groupe n°24 (Voir ici). Selon eux, deux modèles étaient à privilégier, le modèle $k-\epsilon$ et le modèle $k-\omega$. Nous nous attacherons ici à détailler brièvement ces deux modèles et les équations qui les composent. Les valeurs utilisées ainsi que les conditions limites appliquées pour les nouveaux paramètres de la turbulence seront explicités.

Rappel sur les modèles utilisés


Les modèles de turbulence sont basés sur des équations de transport. Il existe des modèles à une équation, à deux équations ($k-\epsilon$ et $k-\omega$ par exemple) et plus encore (Reynolds Stress à 7 équations etc...). Les modèles à deux équations sont basés sur l'hypothèse de Boussinesq qui suppose que le tenseur de Reynolds qui apparaît lorsque l'on réécrit les équations de Navier-Stokes en utilisant une décomposition de Reynolds des grandeurs et en les moyennant, est proportionnel au tenseur des vitesses de déformations.

Le modèle $k-\epsilon$


Le modèle $k-\epsilon$ est le modèle le plus commun et le plus utilisé en CFD, il est un peu le "couteau-suisse" des modèles de turbulence. C'est un modèle à deux équations de transport, une pour l'énergie cinétique turbulente $k$ et une pour le taux de dissipation $\epsilon$.

Dans ce modèle, la viscosité turbulente s'écrit ainsi :

$\nu_t=\rho C_{\mu} \frac{k²}{\epsilon}$

Ces deux équations sont les suivantes :

$\frac{\partial k}{\partial t}+ U_k \frac {\partial k}{\partial x_k}=\frac {\partial}{\partial x_k}\left[ \left( \nu+ \frac {C_{\mu} k²}{\sigma_k \epsilon}\right)\frac {\partial k}{\partial x_k}\right]+\frac {C_{\mu} k²}{\sigma_k \epsilon}\left(\frac {\partial U_i}{\partial x_k}  +\frac {\partial U_k}{\partial x_i}\right)\frac {\partial U_i}{\partial x_k}  -\epsilon$

$\frac{\partial \epsilon}{\partial t}+ U_k \frac {\partial \epsilon}{\partial x_k}=\frac {\partial}{\partial x_k}\left[ \left( \nu+ \frac {C_{\mu} k²}{\sigma_{\epsilon} \epsilon}\right)\frac {\partial \epsilon}{\partial x_k}\right]+C_{\epsilon_1}C_{\mu}k\left(\frac {\partial U_i}{\partial x_k}  +\frac {\partial U_k}{\partial x_i}\right)\frac {\partial U_i}{\partial x_k}  -C_{\epsilon_2} \frac{\epsilon²}{k}$

Les paramètres prennent les valeurs suivantes :

$C_{\mu}$ $\sigma_k$ $\sigma_{\epsilon}$ $C_{\epsilon_1}$ $C_{\epsilon_2}$
0.09 1.0 1.22 1.44 1.92

Ce modèle est connu pour avoir des comportements pathologiques dans certains types d'écoulements, notamment ceux où de forts gradients de pression adverses sont présents. Comme tous les modèles de type RANS, la qualité de la modélisation dépend aussi fortement des valeurs de l'énergie cinétique turbulente et du taux de dissipation injectées dans l'écoulement. En effet, il est incapable de générer de la turbulence, et par conséquent de mauvaises valeurs de $k$ ou d'$\epsilon$ vont soit dissiper trop de turbulence, soit en générer beaucoup trop. Il faut donc être très prudent avec les résultats et les paramètres utilisés.

Le modèle  $k-\omega$ SST


Le modèle $k-\omega$ SST est également un modèle à deux équations. Il s'agit toujours de deux équations de transport, une sur $k$, l'énergie cinétique turbulente, et une autre sur $\omega$ qui peut être vu comme une fréquence caractéristique de la turbulence. Il s'écrit aussi $\omega=\frac{\epsilon}{k}$.

La viscosité cinématique turbulente s'écrit :

$\nu _T = \frac{a_1 k}{max(a_1 \omega, S F_2)}  $

Et les deux équations de transport sont les suivantes :

$\frac{\partial k}{\partial t} + U_j \frac{\partial k}{\partial x_j } = P_k - \beta ^* k\omega + \frac{\partial}{\partial x_j } \left[ {\left( {\nu + \sigma_k \nu _T } \right){{\partial k} \over {\partial x_j }}} \right]$

${{\partial \omega } \over {\partial t}} + U_j {{\partial \omega } \over {\partial x_j }} = \alpha S^2 - \beta \omega ^2 + {\partial \over {\partial x_j }}\left[ {\left( {\nu + \sigma_{\omega} \nu _T } \right){{\partial \omega } \over {\partial x_j }}} \right] + 2( 1 - F_1 ) \sigma_{\omega 2} {1 \over \omega} {{\partial k } \over {\partial x_i}} {{\partial \omega } \over {\partial x_i}}$

Les différents termes introduits dans ces deux équations sont les suivants :

$F_2=\mbox{tanh} \left[ \left[ \mbox{max} \left( { 2 \sqrt{k} \over \beta^* \omega y } , { 500 \nu \over y^2 \omega } \right) \right]^2 \right]$

$P_k=\mbox{min} \left(\tau _{ij} {{\partial U_i } \over {\partial x_j }} , 10\beta^* k \omega \right)$

$F_1=\mbox{tanh} \left\{ \left\{ \mbox{min} \left[ \mbox{max} \left( {\sqrt{k} \over \beta ^* \omega y}, {500 \nu \over y^2 \omega} \right) , {4 \sigma_{\omega 2} k \over CD_{k\omega} y^2} \right] \right\} ^4 \right\}$

$CD_{k\omega}=\mbox{max} \left( 2\rho\sigma_{\omega 2} {1 \over \omega} {{\partial k} \over {\partial x_i}} {{\partial \omega} \over {\partial x_i}}, 10 ^{-10} \right )$

$\phi = \phi_1 F_1 + \phi_2 (1 - F_1)$

Les constantes du modèle sont définies dans le tableau suivant :

$\alpha_1$ $\alpha_2$ $\beta_1$ $\beta_2$ $\beta^*$ $\sigma_{k1}$ $\sigma_{k2}$ $\sigma_{\omega_1}$ $\sigma_{\omega_2}$
$0.55$ $0.44$ $0.075$ $0.0828$ $0.09$ $0.85$ $1$ $0.5$ 0.856

Ce modèle a été développé par Menter et cherche à bénéficier de l'avantage du modèle $k-\omega$ en proche-paroi et du modèle $k-\epsilon$ dans la partie principale de l'écoulement, évitant ainsi les comportements indésirables des deux modèles dans l'écoulement.


Paramètres de la simulation

La simulation Fluent réalisée par SETEC a été lancée avec une intensité turbulente de $1\%$. A partir de là, on peut en déduire les valeurs pour $k$ et $\epsilon$ pour obtenir une intensité équivalente :

$k=\frac{3}{2}(U_\infty I)^2$

$\epsilon=\frac{C_{\mu}^{\frac{3}{4}}k^{\frac{3}{2}}}{l}$

 

On obtient alors les valeurs suivantes pour $k$ et $\epsilon$ : $k=0.06~m^2.s^{-2}$ et $\epsilon=0.35~m^2.s^{-3}$.

Ces valeurs permettent d'initialiser les champs de $k$ et d'$\epsilon$ au début de la simulation. La vitesse initiale dans le domaine est initialisée à $20~m/s$ afin de simuler un courant d'air dans le tunnel.

Pour la simulation $k-\omega$ SST, on réalise les mêmes calculs, et on trouve $k=0.06~m^2.s^{-2}$ et $\omega=64.8s^{-1}$.

Conditions limites

On rappelle que l'intensité turbulente vaut $1\%$.

La figure suivante répertorie les différentes conditions limites utilisées dans cette simulation.

Pour la simulation avec $k-\epsilon$, les conditions limites sont les suivantes :

Champs Right Left farFieldLeft farFieldRight farFieldMoving Hull
k fixedValue 0.06 fixedValue 0.06 kqRWallFunction kqRWallFunction kqRWallFunction kqRWallFunction
$\epsilon$ fixedValue 0.35 fixedValue 0.35 epsilonWallFunction epsilonWallFunction epsilonWallFunction epsilonWallFunction
$\nu_t$ calculated calculated nutkWallFunction nutkWallFunction nutkWallFunction nutkWallFunction

Pour la simulation avec $k-\omega$ SST :

Champs Right Left farFieldLeft farFieldRight farFieldMoving Hull
k fixedValue 0.06 fixedValue 0.06 kqRWallFunction kqRWallFunction kqRWallFunction kqRWallFunction
$\omega$ fixedValue 64.8 fixedValue 64.8 omegaWallFunction omegaWallFunction omegaWallFunction omegaWallFunction
$\nu_t$

calculated

calculated nutLowReWallF* nutLowReWallF* nutLowReWallF*

nutLowReWallF*

nutLowReWallF* : nutLowReWallFunction

Notre $y+$ proche des parois étant toujours  élevé, des modèles de loi de paroi doivent être utilisés, d'où les conditions limites pour $k$,$\omega$,$\epsilon$ et $\nu_t$ aux parois.


Finalement, une fois la turbulence implantée, et le script prêt, il ne reste plus qu'à lancer les simulations finales.

Résultats : Analyse et Comparaisons

Analyse et comparaisons des résultats obtenus sur OpenFoam avec l'étude numérique réalisée sur Ansys Fluent


Dans le cadre de ces simulations du mouvement d'un train dans un tunnel, nous avons décidé de faire déplacer le train à une vitesse de 40 m/s, comme indiqué par nos responsables industriels. De plus, la simulation dure 8 secondes, ce qui permet au train de se déplacer sur une distance de 320 m.

Simulation OpenFOAM avec le modèle k-epsilon

Pour ce modèle $k-\epsilon$, nous avons suivi les conseils de nos collègues ayant réalisé une étude approfondie sur les différents modèles de turbulence. Nous rappelons par conséquent que nous avons imposé une intensité turbulente égale à 1% dans tout le domaine initial, ainsi qu'en tant que condition limite sur tous les bords du tunnel.

La simulation peut donc être lancée, pour une durée physique égale à 8 secondes et un train se déplaçant à une vitesse de 40 m/s.

Les autres paramètres numériques de la simulation sont récapitulés ci-dessous :

Maillage : 105525 cellules

Modèle : $k-\epsilon$

Intensité turbulente : 1%

$k=0.06 m²/s²$

$\epsilon=0.35 m²/s²$

$ Schéma temporel : Euler implicit

Schéma pour les termes de Laplacien/Gradient : Gauss ($2^{ème}$ ordre)

Schéma pour les termes de divergence : Upwind

$\Delta t=0.0001$s

$t_f=8s$

 

Sur ce champ de vitesse, on peut observer l'établissement d'une traînée ce qui est tout à fait caractéristique d'un déplacement d'un objet solide dans un fluide.

Concernant le champ de pression au cours du temps, on peut voir l'apparition d'une surpression à l'avant du train. Là encore, cela paraît a priori rassurant dans la mesure où c'est bien caractéristique de ce type d'écoulement.

On peut donc, après avoir post-traité les résultats sur paraview, faire afficher les champs d'énergie cinétique turbulente à l'instant final de la simulation, c'est-à-dire pour $t=8s$. Cette énergie cinétique turbulente, représentée par la variable $k$, est directement reliée à la création de turbulence : c'est cette variable qui permet de rendre compte des zones de l'écoulement où la turbulence est présente.

Voici ce que nous obtenons pour la zone proche du train :

 

On peut donc voir que, comme bien souvent, la turbulence trouve sa source dans le sillage d'un objet qui se déplace, ici le train. C'est en effet une zone de fort cisaillement, ce qui permet un fort développement de turbulence.

Aussi, en zoomant un peu plus sur la zone "inter-wagons", on peut observer aussi une importance source de turbulence.

Encore une fois, ce phénomène est plutôt cohérent dans la mesure où cette zone crée une vraie discontinuité dans la géométrie et empêche l'écoulement de s'organiser de manière ordonnée.

 

 

Simulation OpenFOAM avec le modèle k-oméga SST

 

Pour ce modèle $k-\omega SST$, nous avons également suivi les conseils de nos collègues. Par conséquent, nous avons encore imposé une intensité turbulente égale à 1% dans tout le domaine initial, ainsi qu'en tant que condition limite sur tous les bords du tunnel.

La simulation peut donc être lancée, pour une durée physique égale à 8 secondes et un train se déplaçant à une vitesse de 40 m/s.

Les autres paramètres numériques de la simulation sont récapitulés ci-dessous :

 

Maillage : 105525 cellules

Modèle : $k-\omega$

Intensité turbulente : 1%

$k=0.06 m²/s²$

$\omega = 0.35 m²/s²$

 Schéma temporel : Euler implicit

Schéma pour les termes de Laplacien/Gradient : Gauss ($2^{ème}$ ordre)

Schéma pour les termes de divergence : Upwind

$\Delta t=0.0001$s

$t_f=8s$

 

Avec ces deux animations de champ de vitesse et de pression, on peut voir des comportements assez similaires à ceux obtenus avec le modèle $k-\epsilon$. Il y a en effet développement d'une traînée ainsi qu'une apparition d'une surpression à l'avant du train. Par contre, quelques phénomènes assez étranges apparaissent sur le champ de vitesse au dessus du train qui, il faut le concéder, nous semblent très étonnants et difficiles à expliquer.

Comme pour le modèle $k-\epsilon$, on peut tracer par l'intermédiaire de Paraview le champ de l'intensité de l'énergie cinétique turbulente. Cela permet de se rendre de compte dans quelle zone de l'écoulement le modèle $k-\omega ~SST$ crée de la turbulence, et de comparer les valeurs par rapport au modèle $k-\epsilon$.

Voici ce que l'on obtient sous paraview :

​Comme précédemment, la turbulence est créée de manière importante dans le sillage du train ainsi que dans la zone séparant chaque wagon. Concernant l'amplitude de cette éergie cinétique turbulente, le modèle $k-\omega ~SST$ crée plus de turbulence dans le sillage, et moins entre les wagons, par rapport au modèle $k-\epsilon$. Cela est difficilement justifiable, étant donné que ce genre d'écart est tout à fait caractéristique d'un changement de modèle de turbulence dans un calcul CFD.

 

Comparaison avec une étude réalisée sous Ansys Fluent

Après avoir fait quelques remarques indépendamment sur chacune des deux simulations complètes sur openFOAM que nous avons réalisées, on peut essayer de comparer ces résultats entre eux. A cela, nos responsables industriels nous ont donné accès à des analyses déjà réalisées sur Ansys Fluent.

Les méthodes numériques de cette analyse étaient les mêmes que celles que nous avons utilisées avec OpenFOAM, mis à part la discrétisation des termes de divergence qui était effectué avec un schéma d'ordre 1 amont au lieu d'un ordre 2 de type Gauss.

Il s'agissait d'un modèle $k-\epsilon$ où on imposait aussi une intensité turbulente initiale, et aux bords, de 1%.

Aussi, et il s'agit d'une remarque plus importante, une différence a lieu concernant les conditions limites en pression. Dans nos simulations OpenFOAM, nous imposions une pression entrée et sortie du tunnel égale à 0 Pa. Sur le calcul Fluent, les conditions limites en entrée et sortie du domaine était du type "Pressure Inlet" et "Pressure Outlet" par conséquent, il se peut qu'il y ait des décalages dans les valeurs de la pression, il faudra donc analyser les résultats en termes de différence de pression entre deux points.

 

Le premier point auquel nous nous sommes intéressés est les profils de vitesse dans le sillage du train. Pour cela, avec OpenFOAM, nous avons défini des lignes verticales dans le sillage du train, avec un fichier classique sampleDict. Sur ces lignes, nous avons demandé à OpenFOAM d'écrire dans un fichier la vitesse horizontale. Ce résultat sera observé pour $t=3s$.
Le deuxième point sur lequel nous nous sommes penchés est le profil de pression le long du tunnel, au niveau du toi. Cela permet de se faire une idée assez précise des variations de pression dans le tunnel, et l'influence que détient le mouvement du train dans un espace confiné sur cette pression dans le tunnel. Le fait de tracer cette pression sur le toit du tunnel est un résultat très pertinent pour les concepteurs de tunnel qui souhaitent déterminer au mieux les surpressions et dépressions qu'engendre le mouvement d'un train dans un tunnel.

 

Voici schématiquement où sont tracés ces profils, en rappelant que nous avons relevé les résultats pour $t=3s$.

 


Concernant le profil vertical de vitesse horizontale, à 1m du train dans son sillage, voici ce que donnes les résultats (pour $t=3s$)

 

On peut donc remarquer que dans la partie basse du tunnel, les écoulements est bien entraîné par le train qui se déplace à 40 m/s. On peut même repérer certaines zones où la vitesse est encore plus importante, ce qui correspond, dans le cas où le train est fixe et un écoulement est imposé en sens inverse, à des recirculations.
Les modèles $k-\epsilon$ d'OpenFOAM et de Ansys Fluent donnent des résultats relativement similaires. Pour ces deux cas, on peut voir que l'écoulement se situant à une hauteur de plus de 5m, correspond à un écoulement turbulent avec un profil quasi-uniforme (mis à part au bord). Le modèle $k-\omega SST$, quant à lui, se démarque davantage avec un profil plus abrupte à la hauteur du train. Cela est relativement cohérent car il est bien connu que le modèle $k-\epsilon$ a souvent tendance à lisser les gradients de vitesse, ce qui étale au cours du temps la zone influencée par le déplacement du train.


Concernant le profil vertical de vitesse horizontale, à 5m du train dans son sillage, voici ce que donnes les résultats (pour $t=3s$)

On observe globalement les mêmes résultats que ceux obtenus pour le profil situé à 1m derrière le train. Seul le modèle $k-\omega SST$ présente des différences importantes avec l'autre profil. La zone influencée par le train est encore plus mince, ce qui confirme que les gradients de vitesse restent plus significatifs qu'avec le modèle $k-\epsilon$

 


Le dernier profil tracé est celui de la pression le long du tunnel, au niveau du toit.

Précision tout d'abord que pour les 3 cas, le train se situe entre $X=150m et X=220m$. On peut alors remarquer que, à chaque fois, une dépression importante a lieu au niveau du train. On peut aussi voir la présence des 4 wagons avec 4 paliers visibles sur le profil de pression.

Aussi, l'amplitude de cette dépression peut paraître plus importante avec le calcul Fluent. Cependant, si on prend en compte le fait que la pression dans la zone arrière du train est plus faible sur Fluent qu'avec les calculs OpenFOAM, cette première idée d'un saut de pression beaucoup plus grand devient discutable. Ce décalage de pression entre la simulation Fluent et les calculs OpenFOAM peut probablement venir de la différence dans les conditions aux limites de pression qui ne sont pas identiques, comme précisé plus haut.

Remarquons aussi la présence d'un pic de pression avec le calcul Fluent. Nous avons eu quelques difficultés à justifier la présence de ce dernier. Après quelques discussions avec notre responsable industriel, ce saut de pression pourrait s'expliquer par une discontinuité du maillage. En effet, à ce niveau, la taille caractéristique des mailles change radicalement et peut éventuellement expliquer cet étrange phénomène, probablement non physique.

Enfin, on retrouve encore le lissage des gradients du modèle $k-\epsilon$ puisque, là aussi, les pentes ont l'air atténuée par ce modèle. C'est notamment le cas à l'avant du train et surtout à l'arrière de celui-ci.

Conclusion et Perspectives

Synthèse


L'objectif de ce projet, qui était principalement la mise en place d'une méthode sous OpenFOAM pour pouvoir gérer le maillage mobile d'un train qui se déplace dans un tunnel peut être considéré comme atteint.

En effet, par l'élaboration du script qui automatise tout le processus, l'utilisateur n'a qu'à créer un dossier classique OpenFOAM avec tous les paramètres numériques de la simulation. Les simulations s'enchaînent alors toutes seules et l'utilisateur pourra post-traiter les résultats.

Le travail ainsi rendu à TECH-AM Ingénierie et à SETEC Bâtiment leur permettra ensuite de modifier s'ils le souhaitent la géométrie du train et du tunnel, les modèles et schémas numériques, dans de bonnes conditions. Les premiers résultats présentés sur ce site ont quand même permis d'être assez confiants dans les calculs effectués par le code. Les industriels pourront ainsi modifier les paramètres des modèles tout en étant assez optimiste sur la qualité des calculs.

Toutefois, l'utilisation des modèles nous éloigne forcément du cas réel, et le plus important reste d'être conscient de cet écart. Il faut pouvoir justifier les erreurs et en trouver les sources; Dans notre cas, on pourra toujours citer par exemple les interpolations entre chaque sous-simulation qui ont tendance à parfois lisser certains gradients importants.

 

Perspectives


Passer en 3D ...

Créer le domaine représentant le train et le tunnel en trois dimensions aurait tout à fait pu être une étape de ce BEI. OpenFOAM étant naturellement un solveur 3D des équations de Navier-Stokes, notre "seul" travail serait de créer la géométrie du domaine de calcul dans une troisème dimension. Cela concernerait à la fois le fichier blockMeshDict pour avoir plusieurs mailles dans la direction Z, et la géométrie du train décrite dans le fichier .obj et obtenue avec le logiciel FreeCAD.

Les temps de calcul auraient sans doute été beaucoup plus grands. De plus, les résultats n'auraient pas forcément été plus pertinents, comparés à l'effort supplémentaire nécessaire à fournir en ressources numériques et en temps de maillage (ressources humaines).

 

Entrée-sortie du tunnel ...

Un autre aspect qui aurait pu être étudié de façon numérique est l'entrée du train dans le tunnel ainsi que sa sortie. Cet aspect pourrait être en soi un travail qui nécessiterait un investissement conséquent, mais Mr Tekam ne désirait pas avoir un code résolvant ces problèmes d'entrée et de sortie du tunnel.

Cependant, nous avons eu quelques premières réflexions sur ce sujet et nous pouvons mettre en valeur le problème majeur. Vu la construction actuelle de notre maillage mobile, nous avons pu voir que les mailles, dans l'intégralité du domaine, sont mobiles. Or, si on souhaite modéliser aussi l'extérieur du tunnel, il est nécessaire d'avoir deux types de conditions limites. Un type de condition limite correspondrait à l'extérieur du tunnel, à l'air libre, et l'autre correspond à la paroi intérieure du tunnel. Ces conditions limites sont fixes dans l'espace alors que le maillage est mobile. C'est ici que le problème se pose puisqu'il n'est pas évident de gérer des conditions limites fixes avec un maillage mobile (même au niveau des parois).

Nous pouvons proposer deux pistes que l'on aurait pu poursuivre si cela faisait partie de nos objectifs.

La première proposition serait de trouver un moyen d'avoir un maillage fixe sur la partie haute du domaine, là où ne passe pas le train. On pourrait simplement envisager de créer un bloc fixe sur cette partie supérieure avec le fichier blockMeshDict. La partie inférieure serait en tout point similaire à ce que nous avons jusqu'à présent. Simplement, la condition limite pour modéliser soit la paroi du tunnel, soit l'air libre, devra être imposé cette fois sur le maillage fixe, ce qui se fait simplement. Le point compliqué est cependant de pouvoir gérer le lien entre les maillages fixe et mobile. Il faut pouvoir interpoler les valeurs à l'interface entre ces deux maillages,ce que OpenFOAM est capable de gérer si l'on en croit les résultats du tutoriel MixerVessel2D. Cependant, une très grande compréhension de cette interpolation est requise pour pouvoir l'appliquer à notre cas du train, rendant peut-être compliquée sa mise en place.

La deuxième solution serait d'utiliser un "package" que l'on peut trouver sur internet. En effet, un package du nom de "Swak4Foam" (couteau-suisse pour OpenFOAM) contient un outil appelé "Groovy Boundary Conditions". Celui-ci permet d'imposer des conditions limites en fonction de la position dans le domaine, alors que, classiquement, les conditions limites sont définies bloc par bloc. Ensuite, selon la position dans le domaine, on peut imposer une condition de type "Mur" ou de type "Air libre" selon qu'on se trouve dans ou en dehors du tunnel. Si l'on devait aller plus loin, nous serions partis sur cette piste puisqu'il assure une totale compatibilité avec notre solveur PimpleDyMFoam et l'outil de maillage mobile actuellement utilisé.

Remerciements

Remerciements


Nous tenions à remercier Mr. Sydney Tekam de Tech-Am Ingénierie ainsi que Mr. André-Marie Dogbo de SETEC Bâtiment pour nous avoir proposé ce sujet, et nous avoir guidés tout au long de celui-ci afin de le mener à bien. Nous avons en effet eu la chance d'obtenir ce sujet intéressant pour notre BEI, auquel nous avons eu le droit à une grande part de liberté et d'initiatives. 

Nous pensons également aux enseignants-chercheurs Dominique Legendre et Annaïg Pedrono qui nous ont régulièrement suivi pour répondre à nos questions et vérifier que nous suivions le bon chemin durant ce BEI.