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.
$$