Lancer une simulation avec OpenFOAM


Lancer une simulation avec OpenFOAM


 

Comme on l'a expliqué précédemment, il y a trois grandes étapes pour le lancement d'une simulation sous OpenFOAM : le maillage, le paramétrage de la simulation, et le lancement du calcul. Cette partie va décrire ces trois étapes pour la simulation de notre goutte.

De plus, on décrira également comment lancer un calcul en parallèle, afin de réduire le temps de calcul.

 

Avant tout calcul, dans un terminal, il faut taper la ligne suivante :

source /mnt/hmf/OpenFOAM/OF22.sh

Ou, pour éviter de taper cette ligne à chaque fois, créer un alias dans le fichier bashrc. Pour cela il faut être dans son home et taper :  vi .bashrc (ou gedit ou geany), et dans ce fichier : alias OF='source /mnt/hmf/OpenFOAM/OF22.sh'. Il suffira ensuite de taper OF dans un terminal.

Ensuite, il est conseillé de travailler dans le /work afin d'économiser son espace disque (les résultats peuvent prendre de la place). Pour créer un dossier :

mkdir -p /work/${USER}/OpenFOAM/${USER}-{WM_PROJECT_VERSION}/run

Pour s'entrainer, ou même pour avoir un modèle, il est conseillé de copier les tutoriaux existants :

cp -r $WM_PROJECT_DIR/tutorials /work/${USER}/OpenFOAM/${USER}-${WM_PROJECT_VERSION}/run

Voir ici pour plus de détails.

 

 

 

Génération du maillage

Dans un premier temps, il faut mailler le domaine d'étude. OpenFOAM travaille avec le système de coordonnées cartésiennes en 3 dimensions et toutes les géométries sont générées en 3D. Il effectue les calculs en 3 dimensions par défaut, mais on peut le "forcer" à travailler en 2D en utilisant la condition limite empty pour la troisième dimension pour laquelle il n'y a aucune solution nécessaire.

Dans notre cas, le domaine est un rectangle de 2 cm de long par 0.4 de haut, avec une faible épaisseur (0.1 cm), il est représenté sur la figure suivante :

Figure : Domaine de calcul

Le fichier maillage se trouve dans le répertoire constant/polyMesh et se nomme blockMeshDict. On trouve plusieurs parties dans ce fichier :

  • vertices : liste des sommets avec leurs coordonnées (x,y,z)
  • edges : liste des arrêtes reliant deux sommets entre eux (par défaut : ligne droite, mais il y existe d'autres géométries : arc, simpleSpline...)
  • block : comme son nom l'indique, c'est la liste des blocs. L'ordre des points est important. Dans cette partie, se trouve :
    • le numéro des sommets définissant le bloc, dans l'ordre.
    • le nombre de cellules dans chaque direction (Nx, Ny, Nz)
    • le ratio de la taille des cellules (rapport entre la plus petite taille de maille et la plus grande) : il permet de raffiner le maillage dans la direction souhaitée (selon y dans notre cas)
  • boundary (ou patches) : ce sont les conditions limites du domaine
    • patch : type générique
    • wall : mur
    • symmetryPlane : plan de symétrie
    • empty : utile pour avoir des géométries 2D
    • wedge : pour les géométries axi-symétriques
    • cyclic : pour les conditions périodiques

Au final, une fois nos conditions limites et nombres de mailles choisis, notre fichier blockMeshDict est donné ci-dessous.

Figure : Fichier blockMeshDict

Ensuite, la génération du maillage se fait grâce à la commande blockMesh exécutée dans le dossier principal de travail. Afin de regarder à quoi il ressemble, on utilise paraFoam. Le résultat est donné sur la figure suivante.

Figure : Maillage du domaine

Lorsque le maillage est généré, il faut s'occuper des conditions initiales et limites, ainsi que des propriétés physiques.

Conditions Limites et Initiales

Conditions limites et initiales

Les conditions initiales et limites sont définies dans le dossier 0. Dans ce dossier, on trouve les fichiers :

  • U : conditions sur la vitesse
  • p : conditions sur la pression  (pas utilisée ici)
  • p_rgh : pression adimensionnalisée par la masse volumique
  • alpha1.org : 

De plus, dans chaque fichier OpenFOAM, les grandeurs sont données avec leurs unités : pour cela, on donne chaque exposant pour chaque unité du système international au début de chaque fichier. Les unités sont données dans cet ordre :

  1 2 3 4 5 6 7
Unité SI kg m s K mol A Cd
  masse longueur temps température quantité de matière intensité électrique intensité lumineuse

 

  • Vitesse : U

Dans notre cas, les murs sont immobiles, donc les conditions limites pour la vitesse seront juste des ​vitesses imposées à 0. Le fichier U est donné ci-dessous. Il y a bien le champ nul uniforme (initialement) est les murs immobiles.

  • Pression hydrostatique : p_rgh

​Comme pour la vitesse, les conditions limites sur la pression seront des conditions de flux nul à la paroi. Le fichier p_rgh est donné ci-dessous.

  • alpha1.org

Pour notre problème, la méthode Volume of Fluid est utilisée. Cette méthode utilise une fonction discrète qui représente la fraction volumique du liquide : alpha1. Alpha1 représente la distribution gaz/liquide dans le domaine (alpha1 = 0 pour le gaz et 1 pour le liquide). Cette distribution sera initialisée avec la commande funkySetFields qui va écraser le fichier alpha1, c'est pourquoi il n'y a que les conditions limites dans alpha1.org, qui sera alors une sauvegarde. 

Dans un premier temps, on impose un angle de contact statique (ici 90°), par la suite on modifiera cette condition pour imposer un angle de contact dynamique.

Initialisation d'une goutte grâce à l'outil funkySetFields

Dans cette partie, on va expliquer comment initialiser une goutte sphérique grâce à l'outil funkySetFields (issu de swak4Foam, téléchargeable ici). 

 

Pour cela, on modifie le fichier funkySetFieldsDict dans le répertoire system​ de la manière suivante :

Dans ce fichier, on définit dans un premier temps alpha1 = 0 partout, ce qui signifie qu'il y a du gaz dans tout le domaine (dans notre cas l'air). Puis on initialise une goutte centrée en 0, en définissant une condition sur un cercle de rayon 0.01 m. Pour chaque cellule, si elle est dans ce cercle alors alpha1 = 1, sinon alpha1 =0.

 

Avant de lancer tout calcul, il faut copier le fichier alpha1.org dans le fichier alpha1 puisque l'initialisation avec funkySetFields va écraser le fichier alpha1. Il faut donc taper la commande (dans le répertoire principal):

cp 0/alpha1.org 0/alpha1

Ensuite, pour initialiser la goutte, il suffit de taper la commande ( toujours dans le répertoire principal) :

funkySetFields -time 0  ou funkySetFields -latestTime

​Cette fonction va donc modifier le fichier alpha1 et le remplir de 1 et de 0, une valeur pour chaque cellule du domaine de calcul, en se basant sur la condition posée dans le fichier funkySetFieldsDict. 

On visualise ensuite notre goutte initiale avec paraFoam.

On a bien une goutte centrée en 0, de rayon 0.01m.

Prise en compte de l'angle de contact

L'initialisation de la goutte avec funkySetFields écrase le fichier alpha1, comme on l'a dit précédemment. Mais dans certains cas, il peut être nécessaire d'imposer un angle de contact sur les wall du domaine. Or dans ce cas, l'option funkySetFields ne permet pas de garder l'angle de contact définit dans alpha1.org. Il faut alors modifier à la main le fichier alpha1 pour imposer l'angle de contact. Comme on lance beaucoup de simulation, on a choisi de faire un script qui prend toutes les commandes en compte sauf le lancement du calcul (interFoam).

  • Imposer un angle de contact statique : Un angle de contact statique doit être imposé lorsque la ligne de contact est immobile. Il faut, sur les frontières où l'on souhaite imposer cette contrainte, choisir :
    • type constantAlphaContactAngle : angle de contact statique,
    • theta0 : valeur de l'angle de contact en degré (ici 90°)
  • Imposer un angle de contact dynamique: Lorsque la ligne de contact bouge, il faut alors imposer un angle dynamique, ce qui est plus compliqué, il y a plus de champs à compléter. Cet angle dynamique a un angle d'avancement et de recul, et nécessite d'imposer la vitesse de la ligne de contact :
    • type dynamicAlphaContactAngle : angle de contact dynamique
    • theta0 : ici, c'est l'angle de contact "à l'équilibre", c'est -à-dire l'angle de contact final lorsque la goutte a fini de se propager
    • uTheta : vitesse de la ligne de contact
    • thetaA : angle d'avancement
    • thetaR : angle de recul (lors des phénomènes d'hysteresis, la goutte s'étale et se rétracte)

Pour  cet angle de contact dynamique, la formule sur laquelle se base OpenFOAM est $$ \theta = \theta_0  + (\theta_A - \theta_R) tanh (\frac{u_{wall}}{u_\theta}) $$

  • Imposer un angle de contact dynamique avec dépendance temporelle :
    • type dynamicAlphaContactAngle : angle de contact dynamique
    • thetae : angle de contact "à l'équilibre"
    • te : temps final
    • t0 : temps initial
    • theta0 : angle de contact initiale

La formule sur laquelle se base OpenFOAM est $$ \theta = \theta_0 + (t - t_0)* \frac{\theta_e - \theta_0}{t_e - t_0}$$.

Le choix de la condition sur l'angle de contact est délicat est demande une étude physique approfondie sur ce problème de mouillage de la goutte.

  • Script pour garder l'angle de contact dans alpha1 :Ce script simplifie les lancements des calculs :
    • il va nettoyer le dossier principal en supprimant tous les dossiers créés précédemment : foamCleanTutorials
    • puis générer le maillage : blockMesh
    • copier alpha1.org dans alpha : cp 0/alpha1.org 0/alpha1
    • initialiser la goutte : funkySetFields -latestTime
    • modifier le fichier alpha1 de manière à prendre en compte l'angle de contact

              

Proprietés physiques

Pour imposer les proprietés physiques du problème, il faut aller dans le dossier constant. Dans notre cas, on modifie 2 fichiers : transportProperties et g.

  • transportProperties

​​Dans ce fichier, on définit la viscosité cinématique et la masse volumique pour chaque phase. Les valeurs utilisées pour la phase liquide et la tension de surface ont été trouvées dans la publication de Chen (1988). 

De plus, dans ce fichier, on peut choisir d'imposer une loi à la viscosité, autre qu'une loi newtonienne. Toutefois dans notre cas, on utilise une viscosité constante. Il faut faire attention à bien prendre la viscosité cinématique et non la viscosité dynamique.

  • g

La gravité est imposée dans ce fichier.

Maintenant que toutes les constantes sont définies, il est temps de passer au paramétrage du calcul. 

 

Paramétrage du calcul

Pour toutes les simulations numériques, le pas de temps est important pour leur résolution. Celui-ci dépend du nombre de Courant :  $$ Co = v \frac{ \Delta t}{ \Delta x} $$. Ce nombre doit vérifier certains critères de stabilité. De plus, le pas de temps peut être fixé ou variable. Dans notre cas on choisit un pas de temps variable avec l'option adjustTimeStep, en définissant le nombre de courant maximal maxCo et le pas de temps maximal maxDeltaT .

Dans le fichier controlDict, on doit choisir :

  • startTime : le temps de départ de la simulation. Pour le premier calcul, le temps initial est 0. L'un des avantages de OpenFOAM est le fait de pouvoir reprendre un calcul là où il s'est arrêté : par exemple si le calcul s'arrête à 5s mais qu'il n'est pas terminé, il suffit de choisir startTime = 5 pour reprendre le calcul là où il s'est arrêté précédemment.
  • endTime : le temps de fin de simulation : c'est le temps auquel va s'arrêter la simulation.
  • writeControl : cette option contrôle l'instant auquel les données seront sauvegardées. Dans notre cas, comme dt varie, si OpenFOAM doit sauvegarder les résultat à un nombre donné de pas de temps, les temps de sauvegarde seront par conséquent arbitraire. Dans ce cas, il faut choisir l'option adjustableRunTime afin qu'il sauve les résultats à des temps "ronds".
  • adjustTimeStep : cette option permet d'ajuster le pas de temps en se basant sur un nombre de Courant maximal et un nombre de pas de temps maximal.

Le fichier final controlDict est donné ci-dessous.

 

Tous les paramètres sont définis, il est maintenant temps de lancer le calcul.

Lancement du calcul

  • Lancer le calcul sur 1 seul processeur

Il suffit de taper interFoam et le calcul est lancé. Pour ne pas bloquer un terminal à cause du calcul, et pour pouvoir récupérer plus facilement les infos affichées à l'écran, il est utile de les mettre dans un fichier log grâce à la commande :

interFoam > log &

Le calcul est alors en arrière-plan et on peut continuer à utiliser le terminal. Lorsque le calcul sera terminé, il suffira de regarder le fichier log pour garder une trace des résultats, pour avoir le temps de calcul ... Pour voir l'avancement du calcul :

watch ls -t  :  cette commande actualise le dossier toutes les 2 secondes, ce qui permet de voir les dossiers créés au fur et à mesure.

tail -f log  :  cette commande affiche en direct le fichier log, ce qui permet de voir les informations du calcul en cours.

 

  • Lancer le calcul sur plusieurs processeurs

Certains calculs peuvent être très longs, et cela peut être avantageux de les lancer sur plusieurs coeurs. Pour cela, il faut modifier le fichier decomposeParDict (dans le répertoire system) en choisissant :

   - numberofSubdomains : le nombre de sous-domaine du domaine initial, il correspond au nombres de coeurs disponibles pour le calcul.

  -  la méthode de décomposition, ici simpleCoeffs : le nombre de sous-domaines dans chaque direction est donné par le vecteur n, et comme dans notre cas on a une géométrie 2D, la composante z doit rester égale à 1. Comme les ordinateurs utilisés ont 4 coeurs, on décompose notre domaine en quatre sous-domaines, d'où notre n égal à (2,2,1).

Il faut ensuite lancer la décomposition grâce à :

decomposePar

La distribution entre les processeurs est alors montrée à l'écran (voir image). 4 dossiers sont alors créés (1 par processeur) dans le dossier principal. Les résultats seront alors dans chacun de ces dossiers.

      

Ensuite, le lancement du calcul se fait avec la commande :

mpirun -np 4 interFoam -parallel > log &

On peut voir l'avancement du calcul grâce à la commande tail -f log toujours. A la fin du calcul, les résultats seront dans chaque dossier pour chaque processeurs, il faut alors reconstruire pour avoir le domaine entier :

reconstructPar

Une fois le domaine reconstruit pour chaque instant, on a à nouveau les dossiers résultats pour chaque pas de temps.

 

Une fois le calcul lancé et terminé, il faut faire le post-processing.