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.