Simulation de goutte impactante sous OpenFoam

 Simulation de goutte impactante sous OpenFoam

                                                 


Equipe

Swann Thuillet (swann.thuillet@etu.enseeiht.fr)

Audrey Vaubourg (audrey.vaubourg@etu.enseeiht.fr)

 

Encadrant

LEGENDRE Dominique

PEDRONO Annaig

 


Partenaire

Jerome Helie (Continental)
 


 

 

Introduction


Introduction


Le Bureau d'Etudes Industrielles est un projet commun aux différentes options de troisième année hydraulique mécanique des fluides de l'ENSEEIHT. C’est un projet pluridisciplinaire qui permet de travailler sur des problématiques industrielles ou de recherche appliquée.

 


Contexte

Les normes de pollution automobiles tendent à réduire fortement les émissions de particules carbone pour les années 2020. En injection Directe Essence haute pression, une grande part des émissions vient des combustions riches dûes à la présence de liquide sur le piston. La compréhension et la modélisation des impacts spray–piston pendant l'injection est de la responsabilité mondiale du groupe de recherche avancée de Toulouse du groupe Continental Automobile, un des 4 grands fournisseurs de composants automobiles (injecteurs, pompe…).

                                        

Figure : Rayonnement des suies produites par un mouillage paroi dans un moteur     


Objectifs

Ce BEI en collaboration avec le groupe Continental s'intéresse àla simulation d'un impact goutte-paroi. Le groupe de recherche utilise OpenFOAM avec succès pour des études de pulvérisation. L’objectif de notre BEI est d'observer de façon exploratoire la capacité d'OpenFOAM à reproduire le comportement d’un impact goutte – paroi, afin de clarifier certains points tels que l'influence de la taille des mailles ou de l'angle de contact.

Les résultats seront comparés aux résultats expérimentaux des travaux de Chen (1988).

Figure : Profil de la goutte à différents instants (Chen, 1988)

Présentation de OpenFOAM


Présentation de OpenFOAM


 

OpenFOAM (Open Field Operation and Manipulation) est une bibliothèque C++ comprenant des outils pour la simulation numérique principalement axée sur la résolution des équations de la mécanique des fluides. Il possède une bibliothèque de solveurs qui permettent de décrire de nombreuses situations. Il existe par exemple des solveurs pour des écoulements compressibles, des écoulements multiphasiques, des écoulements turbulents. Il est gratuit et téléchargeable ici.

 

Une simulation sous OpenFOAM s'effectue en trois grandes étapes :

 

Un dossier type de calcul sous OpenFOAM se présente de la manière suivante. C'est dans ce répertoire que l'utilisateur travaille et fait les modifications nécessaires au lancement de sa simulation.

 

Figure : Structure d'un répertoire de calcul OpenFOAM 

Le dossier principal est le répertoire de travail, il faut se placer dans celui-ci pour lancer les calculs. Dans ce dossier on trouve :

 

Maintenant que nous avons présenté rapidement OpenFOAM, nous allons faire un tutoriel pour le lancement d'une simulation sous OpenFOAM.

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.

Vérification préliminaires

Avant de se lancer dans des calculs conséquents, il y a des vérifications importantes à faire, en particulier la convergence en maillage. On a en plus vérifié que l'angle de contact simulé avec OpenFOAM ne donne pas de résultats trop aberrants.

 

Convergence en maillage

La convergence en maillage est une partie importante pour un travail de simulation. En effet, un maillage raffiné donne des résultats précis, mais plus le nombre mailles est élevé, plus le temps de calcul sera long. Il est donc nécessaire de trouver un compromis entre un temps de calcul raisonnable et un maillage assez précis. Le principe est simple : on commence à lancer une simulation avec un maillage grossier (ici 200*50) et on relève les grandeurs importantes pour l'étude : le rayon et la hauteur de la goutte. Afin d'être sûrs que la goutte est immobile, les calculs lancés font 500 000 itérations. Ensuite, on augmente le nombre de mailles jusqu'à ce que les résultats convergent. Cela signifie qu'augmenter le nombre de mailles ne change rien à la taille du rayon ou de la goutte.

Afin de calculer le rayon et la hauteur de la goutte, grâce à Paraview, on exporte au format .csv les points correspondants aux coordonnées de la goutte (que l'on a obtenu en effectuant une slice de normale z et un contour pour une valeur de alpha de 0.5) :

Les résultats pour les différents maillages sont donnés dans le tableau suivant. Comme attendu, le temps de calcul augmente lorsque le nombre de mailles augmente. La faible différence de temps de calcul entre un maillage de 50 000 et 100 000 vient de la performance des ordinateurs utilisés : du fait qu'ils sont utilisés par d'autres, ils ne sont pas toujours à 100% de leur capacité et donc le temps de calcul peut varier. Pour les calculs les plus importants, on les a effectué en parallèle sur 4 processeurs afin de les terminer plus rapidement, mais le temps de calcul reste tout de même assez conséquent (85h pour 160 000 mailles).
 

Maillage Nb mailles Hauteur (cm) Rayon (cm) Temps calcul (s) Temps calcul (h)
200*50 10 000 0.0204 0.786 29 462 ~ 8
300*75 22 500 0.02 0.799 47 825 ~ 13
400*50 20 000 0.0203 0.789 42 232 ~ 11
400*100 40 000 0.0202 0.795 80 351 ~ 22
500*50 25 000 0.0203 0.788 53 622 ~ 15
500*100 50 000 0.0198 0.811 121 608 ~ 34
650*150 97 500 0.0198 0.812 135 593 (4 coeurs) ~ 38
800*200 160 000 0.0198 0.811 306 363 (4 coeurs) ~ 85

Ce tableau nous permet de constater une convergence en maillage à partir de 50 000 mailles, puisque les hauteurs sont les mêmes, et les rayons sont très peu différents.  Toutefois, on a constaté que plus le maillage est raffiné, plus les contours sont "bruités", en effet, ils ne sont plus très lisses au niveau de la hauteur, ce qui rend difficile l'interprétation et le calcul des valeurs utiles, on est obligé de faire une interpolation pour avoir le contour de la goutte. C'est pourquoi on a choisi de travailler par la suite avec un maillage 500*100 (soit 50 000 mailles), car c'est le plus lisse, et le plus rapide à converger (si on le fait tourner sur 4 coeurs, il sera à priori encore plus rapide à converger).

         

Figure : Convergence en maillage (zoom sur une zone peu bruitée)

Sur la figure précédente, le contour de la goutte est représenté, avec un zoom sur la zone la plus lisse (la zone la plus intéressante car la moins bruitée).

Figure : Convergence en maillage : hauteur

            

Figure : Convergence en maillage : rayon

Les figures de convergence en maillage pour le rayon et la hauteur permettent de mieux observer la convergence à partir de 50 000 mailles puisque la hauteur est constante et le rayon également. On constate également la corrélation rayon-hauteur, en effet lorsque le rayon augmente, la hauteur diminue, et réciproquement. Ceci est logique puisqu'il y a conservation du volume de la goutte lors du processus de propagation.

Validation de l'angle de contact

Lors du processus de propagation d'une goutte sur une surface solide, celle-ci forme un angle avec la paroi : $\theta$ , qui est l'angle de contact. Celui-ci rend compte de l'aptitude du liquide à s'étaler sur la surface, il représente le mouillage du fluide. Il est lié aux tensions interfaciales par l'équation d'Young.

                                       

Figure : Angle de contact (source)

Comme il a été décrit dans le chapitre précédent, OpenFOAM permet de modéliser plusieurs types d'angles de contact :
  • Angle de contact statique ( $\theta$ = cste)
  • Angle de contact dynamique avec dépendance temporelle
  • Angle de contact dynamique avec angle d'avancement et de recul

Les angles modélisés par OpenFOAM sont détaillés dans les parties suivantes.

Angle de contact statique

Dans un premier temps, on a effectué la validation de l'angle statique simulé par OpenFOAM. Pour cela, comme expliqué dans la partie lancement de simulation, on a entré des paramètres pour spécifier l'angle souhaité.

 

Dans un premier temps, on a fait une vérification simple avec une goutte initialisée avec la gravité et un angle de contact de 90°. Avec cet angle de contact, on est censé avoir le même résultat que pour une simulation sans angle de contact. Le résultat est donné sur la figure suivante.

                     

Figure : Goutte soumise à la gravité, avec et sans angle de contact

On remarque bien qu'il n'y a aucune différence entre les deux contours de goutte. Les points au centre de l'axe des abscisses correspondent à de l'air encapsulé par la goutte au départ de la simulation. Ainsi, l'angle de contact statique semble bien simulé par OpenFOAM.

 

On a testé une deuxième méthode de vérification en initialisant des gouttes sans gravité pour 3 angles de contact différents : 30, 60 et 90. Ces gouttes sont initialisées avec l'angle de contact souhaité et la simulation est lancée pendant 5s de temps physique afin de voir s'il y a des courants parasites ou non. Comme il n'y a pas de gravité, la goutte est censée garder sa forme initiale. Les trois figures suivantes montrent que c'est bien le cas puisque les courbes à 5s sont superposées exactement aux courbes initiales. On a tracé également les tangentes pour chaque angle de contact afin de montrer que l'angle simulé vaut la bonne valeur, et la forme de la goutte théorique (qui n'est qu'une cercle). Les figures montrent bien que les gouttes suivent la forme théorique et ne bougent pas au fil du temps, leur angle de contact est bien simulé et il n'y a pas de courant parasite.

                    

Figure : Goutte sans gravité, angle de contact 30°

       

Figure : Goutte sans gravité, angle de contact 60°

                  

Figure : Goutte sans gravité, angle de contact 90°

Ainsi, l'angle de contact statique simulé par OpenFOAM est bien l'angle statique réel comme le montrent ces quelques vérifications que l'on vient d'effectuer. Pour nos simulations, on va garder un angle statique dans un premier temps.

 

 

Angle de contact dynamique

Comme on l'a précisé dans la partie "lancement d'une simulation" concernant l'angle de contact, OpenFOAM permet de simuler plusieurs types d'angles dynamiques :

  • Angle dynamique avec angles de recul et d'avancement : cet angle nécessite la vitesse de la ligne de contact ($u_{\theta}$) ainsi que l'angle à l'équilibre et se base sur cette formule : 

$$ \theta = \theta_0  + (\theta_A - \theta_R) tanh (\frac{u_{wall}}{u_\theta}) $$

Figure : Angle dynamique avec angle d'avancement et de recul

  • Angle dynamique dépendant du temps : pour celui-ci, il faut connaître l'angle d'équilibre à atteindre mais aussi le temps nécessaire à atteindre cet angle. La relation entre l'angle de contact et le temps est linéaire et se base sur cette formule :

$$ \theta = \theta_0 + (t - t_0)* \frac{\theta_e - \theta_0}{t_e - t_0}$$

Figure : Angle dynamique avec dépendance linéaire temporelle

 

Pour les deux simulations précédentes, l'angle de contact à l'équilibre vaut 0.33° (d'après la publication de Chen), et le temps de simulation est de 1s (car c'est juste pour observer les deux comportements de ces angles), donc la goutte n'a pas fini de se propager entièrement. Les deux types d'angles dynamiques ont à peu près le même comportement au début de la propagation de la goutte (pendant 0.5 s environ), mais leur comportement diffère après ce temps. En effet, tandis que le premier angle a l'air "aplati" assez tôt lors de la propagation, pour le deuxième, il y a apparition de bourrelet à chaque extrémité. Ces deux bourrelets viennent des effets inertiels exercés sur la goutte, afin de satisfaire la conservation du volume, la goutte se creuse au centre. De plus, on observe que la goutte revient sur elle même, elle va se stabiliser plus tard avec le bon angle de contact.

Parmi ces deux options, la première nous parait la plus cohérente avec ce que l'on recherche même si l'idéal aurait été une loi puissance (d'après la publication de Chen). Pour avoir cette loi puissance il aurait fallu modifier les sources des angles de contact, en particulier l'angle de contact avec dépendance temporelle, afin de lui rentrer la loi souhaitée. Nous n'avons pas eu le temps de modifier les sources et avons préféré garder un angle de contact statique qui donnait tout de même des résultats assez cohérents.

Maintenant que les vérifications préliminaires sont effectuées, il est temps de passer au calcul de l'impact goutte paroi.

Calcul symétrique

Le calcul symétrique est souvent utilisé afin de réduire les temps de calcul. La première solution envisagée était la condition limite "symmetryPlane" avec interFOAM, mais l'essai ne fut pas concluant. Une solution a été trouvée a la fin du projet : la condition limite "wedge".

Présentation du problème

Dans notre cas, puisque le problème est symétrique, le domaine peut être divisé en deux afin de diminuer le temps de calcul. Pour cela, le domaine est toujours rectangulaire, mais la condition limite à gauche (sur l'axe) n'est plus un wall mais une symétrie. Le domaine ainsi que les conditions limites pour le problème symétrique sont donnés sur la figure suivante.

                         

Figure : Domaine symétrique

Le fichier funkySetFieldsDict n'a pas à être modifié car la goutte initiale sera au même endroit, mais il n'y aura que la moitié à simuler. La goutte initiale est donnée sur la figure suivante.

 

 

OpenFOAM propose deux solutions pour simuler le problème symétrique : symmetryPlane et wedge.

Problèmes rencontrés avec InterFOAM

Dans un premier temps, nous avons utilisé la condition symmetryPlane de OpenFOAM. Le fichier blockMeshDict est donné ci-dessous avec les conditions limites pour le problème symétrique.

         

Figure : Fichier blockMeshDict problème symétrique

Plusieurs simulations ont ensuite été lancées. Dans un premier temps, des simulations d'impact goutte paroi avec un temps physique de 10s ont été lancées, ce qui laissait à la goutte le temps de tomber et de se propager un peu à priori. Toutefois, plusieurs surprises sont apparues :

  • La goutte montait, malgré une gravité vers le bas et tous les paramètres identiques à la simulation normale
  • La goutte avait un comportement étrange (cf vidéo plus bas)

Ces comportement étranges peuvent avoir plusieurs origines :

- le caractère non carré des mailles peut poser des problèmes à la méthode VOF

- la condition physique imposé par SymmetryPlane

​- des courants parasites qui peuvent apparaître...

 

Afin de regarder si la condition limite symmetryPlane est en cause où non, nous avons initialisé des gouttes avec des angles de contact connus, sans gravité, (comme pour la vérification de l'angle de contact), mais avec une condition limite symmetryPlane à gauche, ce qui initialise seulement une partie de goutte. Le calcul est ensuite lancé pendant 5s de temps physique. Les résultats sont donnés sur les figure suivantes, pour des angles de 30°,60°,90°.

                  

Figure : Goutte sans gravité, angle de contact 30°

                   

Figure : Goutte sans gravité, angle de contact 60°

                   

Figure : Goutte sans gravité, angle de contact 90°

 

Ainsi, contrairement au problème non symétrique ou la goutte ne bougeait pas, ce n'est pas le cas ici. Pour 30 et 60° la différence est très faible et pourrait être due au fait que le contour est pris à 0.5. Mais pour 90° la différence est non négligeable, et on voit donc bien l'influence des courants parasites sur la goutte.

En effectuant des recherches sur le web, nous nous sommes rendu compte que la condition symmetryPlane utilisée avec interFoam posait de nombreux problèmes, et chaque problème était résolu au cas par cas, mais aucun cas ne correspondait au nôtre.

Une solution aurait été de modifier les sources du solveur, au niveau des conditions limites, mais dans un projet si court, nous avons choisi de nous intéresser au calcul 2,5D et 3D qui nous ont paru plus concrets et plus intéressants.

Alternative étudiée : wedge

L'autre méthode qui a été testée pour tirer avantage de la configuration symétrique du problème est d'utiliser la condition limite d'OpenFOAM nommée wedge.

Le domaine est alors défini comme une "part de gâteau" :

 

 

 

 

 

 

 

 

 

Pour cette configuration les points sont définis comme suit :

Le segment à gauche sur l'image est laissé avec le type empty, alors que les faces latérales sont définies avec le type wedge.

Pour cette configuration, nous avons imposé un angle de contact statique de 0,33 puis nous avons obtenu le résultat suivant :

L'utilisation de cette géométrie semble donc adaptée au problème de l'étalement de la goutte.

Il faut cependant noter que ce résultat qui correspond à 1s de temps physique a été obtenu après environ 23 heures de calcul sur 4 coeurs. La longueur du calcul est due au respect du CFL qui impose un pas de temps petit, en effet la cellule au niveau du coin en bas à gauche du domaine est de petite taille et c'est elle qui demande le pas de temps le plus contraignant.

Cette géométrie pourrait permettre de reconstituer en 3D le mouillage d'une goutte en évitant de créer un vrai domaine 3D. Pour cela, il faut s'assurer du comportement isotropique du phénomène, c'est-à-dire que la goutte s'étale de la même façon suivant les directions horizontales. Afin de comprendre mieux la condition de symétrie utilisée "wedge" il est possible de lire les sources du code OpenFoam dans le dossier /mnt/hmf/OpenFOAM/OpenFOAM-2.2.2/src/finiteVolume/fields/fvPatchFields/constraint/wedge

 

Calcul 2,5 D

Du fait que le calcul symétrique n'est pas possible, on va donc lancer des calculs 2,5D avec OpenFOAM (2D et une faible épaisseur suivant z). Mais les résultats fournis par les publications, en particulier les lois pour la hauteur et le rayon en fonction du temps ne sont plus adaptées. Afin de vérifier les résultats obtenus avec OpenFOAM et de les comparer à une loi théorique, on a décidé de reprendre l'étude mathématique du problème de la propagation de goutte pour un cylindre.

 

Lorsque un liquide non volatile est déposé sur une surface homogène, il se déforme jusqu'à atteindre un état d'équilibre. Dans notre cas, la goutte est supposée se propager de manière symétrique (et non axi symétrique comme dans le cas d'une goutte sphérique) sur la surface horizontale. Afin de satisfaire la condition de volume constant, alors que le rayon va augmenter, la hauteur va diminuer. La figure suivante montre un schéma d'une goutte se propageant sur une surface solide. Les grandeurs définies sur ce schéma seront utilisées par la suite dans les équations.

                           

Figure : Schéma d'une goutte de liquide se propageant sur une surface solide 

Formulation mathématique

Comme la goutte étudiée avec OpenFOAM est une goutte 2,5D, donc une goutte "cylindrique", les résultats fournis dans les différents papiers ne sont plus valables. L'analyse mathématique doit donc être refaite à nouveau. On utilise l'étude que Y.Gu et D.Li ont fait (Colloid Surfaces A) pour refaire notre étude.

Pour cette étude, la méthode OEB est employée (Overall Energy Balance), qui est un bilan d'énergie global du système, afin de modéliser le processus de propagation. En général, ce processus dépend des forces inertielles, gravitationnelles et visqueuses ainsi que des tensions interfaciales et de la mouillabilité du système. Dans un premier temps, on va effectuer un bilan d'énergie, qui va nous fournir une équation, puis on utilisera la relation de Laplace afin de trouver l'équation du profil fluide-liquide.

La première équation obtenue par le bilan d'énergie est :  $$ \frac{m}{2} \left[ (\frac{d x_m}{dt})^2 + (\frac{d y_m}{dt})^2 \right] + \gamma_{lf} \left[A_{lf} - A_{sl} cos\theta_e - 2\pi r_0 L\right] +  mg (y_m -r_0) + A* \int_0^t \frac{L}{\theta(t)} \left[\frac{dR}{dt}\right]^2 dt =0  $$

Son obtention est détaillée dans la sous-partie "Bilan d'énergie".

 

La deuxième équation obtenue par l'équation de Laplace est :

$$
    \left \{
    \begin{array}{ll}
    \frac{dy_1}{dx} = y_2 \\
    \frac{dy_2}{dx} = \left[ \frac{\Delta \rho g}{ \gamma_{lf} } (H(t) - y_1 ) + \frac{2}{R_0} - \frac{y_2}{x \sqrt{1 + y_2^2} } \right] * (1+y_2^2 )^{\frac{3}{2}}
    \end{array}
    \right.
$$

Son développement est décrit dans la sous-partie "Profil liquide-fluide".

 

Bilan d'énergie

Pour ce bilan d'énergie, quatre composantes sont prises en compte :

  • L'énergie cinétique : $ E_k (t) $
  • L'énergie potentielle due aux tensions interfaciales : $ E_p (t) $
  • L'énergie potentielle due à la force gravitationnelle : $ E_g (t) $
  • Le travail de dissipation visqueuse : $ L_f (t) $

La conservation de l'énergie implique que l'énergie totale du système ($ E_k + E_p + E_g $) doit décroître pendant la propagation à cause de la dissipation visqueuse ($L_f $). L'équation du bilan d'énergie global (OEB) est donc :

$$ \frac{d}{dt}  \left[ \underbrace{E_k(t)}_\text{énergie cinétique} + \overbrace{
\underbrace{E_p(t)}_\text{tension de surface}
\underbrace{E_g(t)}_\text{force gravitationnelle} }^\text{énergie potentielle} \right] 
+ \underbrace{\frac{d L_f (t)}{dt}}_\text{dissipation visqueuse} = 0 $$

Chaque terme va être calculé afin de déterminer l'équation finale.

  • L'énergie cinétique : $ E_k (t) $

Afin de calculer l'énergie cinétique, on considère que la goutte est composée de plusieurs segments, qui se propagent tous à la vitesse $ U_m (t)$ du centre de masse définie telle que :

$   U_{mx} (t) = \frac{dx_m}{dt} $ et $   U_{my} (t) = \frac{dy_m}{dt} $

Comme la masse volumique du fluide est constante, le centre de masse et le centre géométrique sont les mêmes, définies par :

$  x_m(t) =\frac{ \int_0^{R(t)} x y(x) dx }{ \int_0^{R(t) y(x) dx} } $ et $ y_m(t) =\frac{1}{2} \frac{ \int _0^{R(t)}  y^2(x) dx }{ \int_0^{R(t)}  y(x) dx } $, où $ y(x) $ représente le profil liquide-fluide, c'est-à-dire la la forme de la goutte.

Au final, l'énergie cinétique vaut :

$$ \boxed{ E_k (t) = \frac{m}{2}  \left[ (\frac{d x_m}{dt})^2 + (\frac{d y_m}{dt})^2 \right] } $$

A t= 0, $ \boxed{E_k(0) = 0} $ car la goutte est immobile initialement.

L'avantage d'utiliser le centre de masse est qu'il n'y a pas besoin de connaître le champ de vitesse à l'intérieur de la goutte, ou de supposer une distribution simple qui satisfait l'équation de continuité et les conditions limites; et c'est quasiment impossible de mesurer la vitesse à l'intérieur de la goutte. 

  • L'énergie potentielle due aux tensions interfaciales : $ E_p (t) $

L'énergie potentielle due aux tensions interfaciales  peut s'écrire : $ E_p (t) = E_{lf} (t) + E_{sf} (t) + E_{sl} (t)$ .

Avec $ E_i (t)  = {\gamma}_i*A_i (t) $, où $ E $ est l'énergie potentielle, $ \gamma $ la tension de surface et $ A $ la surface, et les indices $i = lf, sf, sl $ correspondant à l'interface Liquide-Fluide, Solide-Fluide et Solide-Liquide. En utilisant l'équation d'Young ( $ {\gamma}_{sf} = {\gamma}_{sl} + {\gamma}_{lf} cos {\theta}_e $ ) et en sachant que $ A_{sf} (t) + A_{sl} (t) = A_{total} $ (surface du solide) on a :

$$ \boxed{ E_p (t) = {\gamma}_{lf} \left[ A_{lf} (t) - A_sl(t) cos {\theta}_e \right] + {\gamma}_{sf} A_{total} } $$

Avec $ A_{lf} (t) = \int_{-R(t)}^{R(t)} L \frac{dy(x)}{dx} $ et $ A_{sl} = 2 R(t) L $ (où $L$ est l'épaisseur du domaine, suivant z).

A t=0, $ A_{LF} (0) = 2 \pi r_0 L $ et $ A_{sl} (0) =0 $, donc $ \boxed{ E_p (0) = 2 \pi r_0 L {\gamma}_{lf} + {\gamma}_{sf} A_{total}  }$

  • L'énergie potentielle due à la force gravitationnelle : $ E_g (t) $

On se base à nouveau sur l'approche du centre de masse, et donc l'énergie potentielle gravitationnelle de la goutte se propageant vaut donc :  $$  \boxed{ E_g (t) = m g y_m (t) } $$.

A t= 0 , $ \boxed{ E_g (0) = m g r_0} $.

  • Le travail de dissipation visqueuse : $ L_f (t) $

On utilise la formule fournie dans le papier (Y.Gu et D.Li) qui est issue d'une étude précédente, pour la force visqueuse : $ F_v (t) = \frac{3 {\mu}_l }{ {\theta}_d (t) } \frac{dR(t)}{dt} \ln{(\epsilon_\delta^{-1})} $, où $ \mu$ est la viscosité dynamique du fluide, $ \theta $ est l'angle de contact et $\epsilon $ une constante à modifier pour caler au mieux aux résultats numériques. 

La force visqueuse totale exercée par la ligne de contact sur la surface solide vaut $ 2 L F_v (t) $. Par conséquent, le travail total dû à la force visqueuse vaut : $ L_f(t) = \int_0^{R(t)} 2 L F_v (t) dR(t) $.

Finalement : $$ \boxed{ L_f (t) = A * \int_0^t \frac{ R(t) }{ {\theta}_d} \left[ \frac{dR(t)}{dt} \right]^2 dt } $$.

Avec $ A = 6 \pi {\mu}_l \ln{(\epsilon_\delta^{-1})} $, et $ \boxed{ L_f (0) = 0} $.

 

Finalement, l'équation d'équilibre d'énergie au temps t s'écrit $$ \Delta E_k (t) + \Delta E_p(t) + \Delta E_g (t) + \Delta L_f (t) = 0 $$.

En remplacent chaque terme par son expression calculée juste avant : $$ \Rightarrow \frac{m}{2} \left[ (\frac{d x_m}{dt})^2 + (\frac{d y_m}{dt})^2 \right] + \gamma_{lf} \left[A_{lf} - A_{sl} cos\theta_e - 2\pi r_0 L\right] +  mg (y_m -r_0) + A* \int_0^t \frac{L}{\theta(t)} \left[\frac{dR}{dt}\right]^2 dt =0  $$

Cette équation est une équation générale pour modéliser la propagation d'une goutte liquide cylindrique se propageant sur une surface solide. Elle permet de déterminer quantitativement comment le processus de propagation d'une petite goutte de liquide sur une surface solide dépend de la surface et des propriétés physiques du système liquide-fluide-solide, des conditions initiales et de la forme du profil à chaque instant.

Profil liquide-fluide : y(x)

La forme de la goutte est gouvernée par l'équation de Laplace : $$ \Delta P = \underbrace{P_L}_\text{Pression du liquide} - \underbrace{P_F}_\text{Pression du fluide} = \sigma_{LF}(\frac{1}{R_1} + \frac{1}{R_2}) $$ Avec $R_1$ et $R_2$ les rayons de courbure de l'interface.

Sur l'interface, la pression s'exprime grâce à $ P_i  = P_{i0} +\rho_i g \left[ H(t) - y(x) \right] $, où i = l ou f pour le liquide ou le fluide respectivement, $P_{i0}$ la pression de référence, $H(t)$ la hauteur du sommet de la goutte et $y(x)$ la forme de l'interface. La différence de pression peut être déterminée au point $x=0$ et $y=H(t)$ à chaque instant. On a alors  $ \Delta P_0 = P_{l0} - P_{f0} = \frac{ 2 \gamma_{lf}}{R_0 (t)} $, où $R_0(t)$ est le rayon de courbure au sommet de la goutte ( toujours positif dans notre cas).

En substituant les formules des pressions dans l'équation précédente, et en remplaçant les rayons de courbure par leur valeur, on obtient : $$ \frac{ \pm \frac{d^2 y}{dx^2} }{ \left[ 1 + (\frac{dy}{dx})^2 \right]^{\frac{3}{2} } } + \frac{ \pm \frac{dy}{dx} }{x \left[1 + (\frac{dy}{dx})^2 \right]^{\frac{1}{2}} } = \frac{2}{R_0 (t) } + \frac{ \Delta \rho g }{\gamma_{lf} } \left[ H(t) - y(x) \right] $$.

Il faut faire attention aux $ \pm $ pour pouvoir obtenir le profil, en effet, les signes vont changer selon que l'on est sur la partie inférieure ou supérieure de la goutte lorsque l'angle de contact est supérieur à 90° (soit au début de la propagation). Cette équation est une équation du deuxième degré et nécessite donc deux conditions limites :

$$
    \left \{
    \begin{array}{ll}
    y(x=0) = H(t) \\
    \frac{dy}{dx} (x=0) =0
    \end{array}
    \right.
$$

Cette équation doit être transformée en deux équations différentielles du premier ordre, en posant $y_1 = y$ et $ y_2 = \frac{dy}{dx} $.

L'équation devient alors:

$$
    \left \{
    \begin{array}{ll}
    \frac{dy_1}{dx} = y_2 \\
    \frac{dy_2}{dx} = \left[ \frac{\Delta \rho g}{ \gamma_{lf} } (H(t) - y_1 ) + \frac{2}{R_0} - \frac{y_2}{x \sqrt{1 + y_2^2} } \right] * (1+y_2^2 )^{\frac{3}{2}}
    \end{array}
    \right.
$$

Avec les conditions limites suivantes :

$$
    \left \{
    \begin{array}{ll}
    y_1(x=0) = H(t) \\
    y_2 (x=0) =0
    \end{array}
    \right.
$$.

 

Pour résoudre cette équation, on a besoin de $H(t)$ et de $R_0(t)$ . Une fois qu'un des 2 est connu, l'autre se déduit par la conservation du volume ( $ V_1 = \int_0^{H(t)} L x dy = cst $ ).

Une fois que le profil est connu, l'angle de contact est déterminé avec :

$$
    \frac{dy}{dx} (x=R(t)) = \left \{
    \begin{array}{ll}
   \tan[180 - \theta_d (t) ] \text{ si } \frac{dy}{dx} (x=R) > 0 \\
   - \tan[\theta_d (t) ] \text{ si } \frac{dy}{dx} (x=R) \leq 0
    \end{array}
    \right.
$$

Procédure numérique

La résolution des équations présentées précédemment doit permettre de connaître la forme de notre goutte en 2,5 D à tout moment. Cependant, ces équations étant directement liées il n'est pas possible de les résoudre analytiquement. C'est pourquoi une procédure numérique adaptée a été mise en place et présentée par Gu & Li dans leur article A model for a liquid drop spreading on a solid surface.

Les étapes à suivre sont les suivantes :

- Initialisation de la goutte et des grandeurs physiques (rayon, hauteur, angle de contact ...)

- Choix d'un $ \Delta R_0$ tel que $\Delta R_0=R_0(t+\Delta t)-R_0(t)$

- Choix d'un $H(t+\Delta t)$ associé tel que  $0<H(t+\Delta t)< 2r_0$

- Calcul du profil liquide-fluide à $t+\Delta t$ grâce à l'équation :

$$ \frac{ \pm \frac{d^2 y}{dx^2} }{ \left[ 1 + (\frac{dy}{dx})^2 \right]^{\frac{3}{2} } } + \frac{ \pm \frac{dy}{dx} }{x \left[1 + (\frac{dy}{dx})^2 \right]^{\frac{1}{2}} } = \frac{2}{R_0 (t) } + \frac{ \Delta \rho g }{\gamma_{lf} } \left[ H(t) - y(x) \right] $$

 

et les conditions limites :

$$
    \left \{
    \begin{array}{ll}
    y_1(x=0) = H(t) \\
    y_2 (x=0) =0
    \end{array}
    \right.
$$.

- En considérant la conservation du volume de la goutte   $V_l = \int_0^H 2Lxdy$   la hauteur H est ajustée.

- Il est alors possible de calculer $R(t+\Delta t)$, $A_{lf}(t+\Delta t)$, $x_m(t+\Delta t)$, $y_m(t+\Delta t)$ et $\theta_d(t+\Delta t)$

- Enfin, l'équation du bilan d'énergie permet de calculer le pas de temps $\Delta t$ qui correspond à cette variation  $ \Delta R_0$. Pour cela, il suffit d'incrémenter le pas de temps $\Delta t$ jusqu'à ce que l'équation bilan soit vérifiée. Le coefficient $\epsilon_\delta^{-1}$ est choisi au hasard.

- Le temps final où le profil de la goutte n'évolue plus peut être déterminé de deux façons : la valeur de l'angle de contact dynamique est à peu près égale à celle de l'angle de contact d'équilibre ou la vitesse du profil est inférieure à une vitesse minimale fixée.

- Enfin le coefficient $\epsilon_\delta^{-1}$ qui intervient dans le bilan énergétique est choisi de manière à ce que le résultat analytique déterminé corresponde le plus possible aux résultats de la simulation numérique.

 

Il existe plusieurs problématiques qui sont liées à cette procédure :

- Le choix de $ \Delta R_0$ doit être adapté au problème. Si la valeur prise est trop grande, le calcul ne pourra pas converger. A contrario, si cette valeur est trop faible le pas de temps sera lui aussi très faible et le calcul en sera d'autant plus long.

- L'équation qui nous permet d'obtenir le profil de la goutte est décrite avec les signes $\pm$. En réalité il s'agit donc de deux systèmes différentiels à résoudre.

- Pour ajuster la valeur de la hauteur H pour qu'elle vérifie la conservation du volume de fluide, il faut donc incrémenter la valeur de H jusqu'à ce que le volume soit à peu près égal à celui du pas de temps précédent. On a donc une boucle "while" ou "until" dans laquelle il faut résoudre à chaque itération l'équation qui donne le profil de la goutte.

- Pour déterminer le $\Delta t$ qui vérifie le bilan d'énergie, il faut aussi incrémenter une valeur de $\Delta t$ petite jusqu'à ce que cette équation soit vérifiée. Ainsi, on a ici aussi une boucle de type "while" ou "until" où il faut résoudre l'équation d'énergie. Enfin, la valeur de ce pas de temps change pour chaque itération temporelle.

Résolution numérique

Pour la résolution numérique, on va coder les équations obtenues précédemment avec Matlab afin de les résoudre et de déterminer R,H et $\theta$ à chaque instant en se basant sur la procédure numérique expliquée précédemment.

On code dans un premier temps l'initialisation de la goutte avant de coder la résolution à chaque instant.

Initialisation de la goutte

Pour l'initialisation de la goutte, on utilise l'équation trouvée grâce à l'équation de Laplace afin de trouver le profil $y(x)$ initial de la goutte. A cause des $\pm$ on sépare le calcul de la goutte en quatre parties :

                                       

Du fait de la symétrie de la goutte, on calcule uniquement la partie 1 et la partie 4, puis par symétrie, la forme globale de la goutte est déduite.

L'initialisation de la goutte se passe très bien. Le système d'équation du premier ordre est résolu avec la fonction ode45 de Matlab. La goutte initiale calculée par Matlab est donnée sur la figure suivante :

              

Figure : Goutte initiale, résolue mathématiquement

Ainsi, le système d'équation du premier ordre issu de l'équation de Laplace de la capillarité est bien résolu par Matlab puisqu'on retrouve bien une goutte initiale d'un rayon de 0.02 cm.

Boucle temporelle

Après avoir vérifié que le système d'équation du premier ordre pour la détermination du profil liquide-fluide est bien résolu grâce à l'initialisation de la goutte, on s'occupe désormais de la résolution de la première équation issue de l'équilibre énergétique.

Afin de résoudre la première équation, obtenue par le bilan d'énergie, on la discrétise en temps et en espace. On résout comme pour la goutte initiale, la partie inférieure et la partie supérieure. Le comportement de chaque partie se passe plutôt bien, le seul problème vient du raccordement entre les deux parties, qui ne se fait pas, comme on peut le constater sur la figure suivante.

                    

C'est donc ici que se situe la difficulté qui ne nous permet pas d'avoir une solution analytique du mouillage d'une goutte de forme "cylindrique" (en 2,5 D).

Le profil liquide-fluide n'ayant pas pu être déterminé, il n'est pas possible de mettre en place la suite de la procédure numérique.

 

Simulation

Le calcul est lancé avec les paramètres physiques donnés dans la publication de Chen, pour une goutte ayant un rayon initial 0.0337 cm, pour un temps physique de 120s afin d'obtenir une loi à partir des résultats de simulation, et pour voir si la différence avec les résultats de Chen sont très éloignés ou non. La propagation de la goutte est donnée sur l'animation suivante.

 

On a également tracé les évolutions du rayon, de la hauteur et de l'angle de contact afin d'obtenir des lois en fonction du temps. On suppose que ce sont des lois en puissance de type $ \alpha t^\beta $. La figure suivante montre les évolutions pour les trois grandeurs.

Figure : Évolution des différentes grandeurs 

Grâce à une fonction d'optimisation de Matlab, on détermine les coefficients des lois en puissance, et on constate que ces lois sont adaptées. Les coefficients a, c et e respectivement pour le rayon, la hauteur et l'angle de contact pourraient être modifiés pour mieux caler les courbes (avec le coefficient $\epsilon $), mais on constate que cela colle déjà plutôt bien. Le fit commence à 1 seconde car l'évolution avant n'est pas toujours régulière et ne varie pas toujours en loi puissance.

 

Nous avons aussi tracé les composantes énergétiques de l'équation d'équilibre énergétique en fonction du temps afin de voir leur évolution. Cette figure montre que les plus grandes variations d'énergie potentielle due aux tensions interfaciales et de dissipation visqueuse ont lieu au début de la propagation de la goutte. De plus, on constate que l'énergie cinétique et l'énergie potentielle gravitationnelle sont négligeables, même pendant la période initiale de la propagation. Ainsi, la diminution de l'énergie potentielle due aux tensions interfaciales compense uniquement la dissipation visqueuse afin que la goutte se propage.

            

Figure : Composantes énergétiques

Pour un même rayon de goutte initiale (0.0337 cm), on a comparé nos résultats pour la loi en puissance avec les résultats expérimentaux des travaux de Chen pour la hauteur et le rayon, ils sont résumés dans les tableaux suivants.

a b $\alpha$ $\beta$
0.09 0.195 0.087 0.107

Tableau : Simulation (a,b) vs Expérimental ($\alpha$, $\beta^$) pour le rayon

c d $\alpha$ $\beta$
0.029 -0.188 0.015 -0.231

Tableau : Simulation (c,d) vs Expérimental ($\alpha$ ,$\beta$) pour la hauteur

Ainsi, bien que la différence soit faible, les évolutions sont bien différentes, donc l'analyse mathématique était bien à refaire afin de trouver l'évolution théorique pour la hauteur et le rayon.

Calcul 3D

La possibilité d'un calcul en 3D a été envisagée alors que la configuration symétrique ne fonctionnait pas. L'objectif était d'obtenir un résultat qui serait comparable à ceux publiés par Chen. Dans un premier temps la géométrie créée est la suivante :

Ce maillage comporte 100 mailles dans chaque direction horizontale et 50 en hauteur.

La calcul pour cette configuration a demandé environ 15 heures pour 1s de temps physique. De plus, les résultats ne sont pas satisfaisants physiquement, le volume de la goutte ne se conservant pas. Cette erreur est sûrement due au maillage qui n'est pas assez raffiné. Pour faire un calcul avec un maillage qui correspondrait au maillage optimal obtenu pour le 2D (c'est-à-dire 500x500x100), l'estimation grossière du temps de calcul sur 4 coeurs d'un HP Z400 conduit à environ 750 heures.

Nous n'avons donc pas eu le temps dans ce BEI d'étudier le cas avec ce maillage mais il est possible d'envisager cette étude avec des machines plus puissantes et plus de temps.

 

Conclusion

Ce projet nous a permis de découvrir OpenFOAM, et plus particulièrement le solveur interFoam. Il est capable de bien reproduire le comportement d'une goutte se propageant sur une surface solide, en respectant l'angle de contact.

Ce travail consiste en un travail préliminaire qui nous a permis de déterminer un maillage optimal qui ne nécessite pas un temps de calcul trop élevé mais qui donne des résultats précis. Une solution pour un calcul symétrique (par conséquent plus rapide) a été proposée, bien qu'elle ait été trouvée vers la fin du projet et qu'elle n'a pas pu être approfondie. Enfin, le modèle pour l'angle de contact doit être modifié directement dans les sources afin de prendre en compte une loi en puissance.

Par ailleurs, le calcul 3D qui pourrait nous permettre de comparer nos résultats avec ceux des publications n'est pas réalisable dans le temps du projet, il faudrait plus de temps pour avoir des résultats.

Toutefois, des modèles ont pu être déterminés à partir des résultats des simulations. La résolution mathématique n'a pas abouti, mais la forme générale des évolutions du rayon et de la hauteur de la goutte ont été trouvées.

 

 

Nous souhaiterions remercier Annaig Pedrono pour son aide et ses conseils sur OpenFOAM tout au long du projet, ainsi que Jerome Helie pour avoir proposé ce sujet et nous avoir guidé durant celui-ci.

Références