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.