# Simulation numérique d'écoulements liquide-gaz en approximation bas Mach

V. DARU<sup>a,b</sup>, M.-C. DULUC<sup>a,c</sup>, P. LE QUERE<sup>a</sup>, D. JURIC<sup>a</sup>

a. LIMSI-CNRS, UPR3251, BP133, Bât. 508, 91403 ORSAY CEDEX

b. Arts et Métiers ParisTech, Lab. SINUMEF, 151 bd de l'hôpital, 75013 PARIS

c. CNAM, 292 rue Saint Martin, 75141 PARIS Cedex 03

## Résumé :

On s'intéresse dans ce travail à la simulation numérique d'écoulements liquide-gaz. La phase liquide est considérée comme incompressible, tandis que la phase gazeuse est supposée compressible dans l'approximation d'un écoulement à bas Mach. On présente un modèle et une méthode numérique aptes à traiter ce genre de situation. La méthode utilise un suivi d'interface Lagrangien de type « front-tracking ». On effectue des simulations 2D sur une configuration de bulles d'air incluses dans un liquide confiné, et soumises initialement à une différence de pression thermodynamique. Un mouvement oscillatoire des bulles est alors engendré, au cours duquel les bulles subissent alternativement des compressions et des dilatations, associées à un échauffement ou un refroidissement. Les oscillations sont finalement amorties. On montre ici que la conductivité thermique du liquide a un effet significatif sur la fréquence et le temps d'amortissement des oscillations.

## Abstract :

This work is devoted to the numerical simulation of liquid-gas flows. The liquid phase is considered as incompressible, while the gas phase is treated as compressible in the low Mach number approach. We present a model and a numerical method aimed at the computation of such two-phase flows. The numerical model uses a lagrangian front-tracking method to deal with the interface. Numerical results in the 2-D case are presented. Two air bubbles are enclosed in a rigid cavity and surrounded with liquid water. As the initial pressure of the two bubbles is set to different values, an oscillatory motion is induced in which the bubbles undergo alternate compression and dilatation stages associated with heating and cooling processes. This oscillatory motion cannot be sustained and a damping is finally observed. It is shown in the present work that thermal conductivity of the liquid has a significant effect on both the frequency and the damping time scale of the oscillations.

## Mots clefs : écoulements diphasiques, interfaces, compressibilité, simulation numérique

## 1 Introduction

La simulation numérique d'écoulements diphasiques liquide-gaz à phases séparées pose de nombreuses difficultés, concernant par exemple la représentation des interfaces ou la prise en compte de la tension de surface. De nombreux travaux ont été menés sur ce sujet dans le cas incompressible isotherme. Cependant des applications spécifiques comme les écoulements en milieu confiné (comme c'est souvent le cas en microfluidique ou dans l'ingénierie nucléaire) nécessitent de prendre en compte la compressibilité de la phase gazeuse. En effet les variations de la pression thermodynamique dans la phase gazeuse jouent alors un rôle important. Ceci amène des difficultés supplémentaires car il faut prendre en compte simultanément des régimes d'écoulement compressible et incompressible.

Ce travail concerne la simulation numérique de tels écoulements, avec des propriétés physiques réalistes et des rapports de masse volumiques élevés entre phases. Dans beaucoup d'applications les vitesses d'écoulement dans le gaz restent faibles, et une approche de type bas Mach est bien adaptée. Le liquide, lui, peut être considéré comme incompressible, et on considérera donc un modèle dans lequel coexistent un écoulement purement incompressible et un écoulement bas Mach, séparés par une interface
dynamique à travers laquelle les propriétés physiques du fluide sont discontinues.

On présente en section 2 le modèle et la méthode numérique développés pour simuler ces écoulements. On considère ensuite en section 3 une application bidimensionnelle comportant deux bulles d'air incluses dans une enceinte fermée remplie d'eau. Les bulles étant initialement soumise à une différence de pression thermodynamique, un mouvement oscillatoire se développe avant d'être amorti. Les bulles subissent alternativement des compressions et des dilatations, associées à un échauffement ou un refroidissement. On mène une analyse permettant d'identifier les temps caractéristiques des phénomènes et de préciser l'influence de la conductivité thermique du liquide.

# 2 Modélisation et méthode numérique

## 2.1 Modélisation

On considère donc un écoulement liquide-gaz, le liquide étant incompressible tandis que le gaz est compressible. Le nombre de Mach dans le gaz est supposé petit. Le modèle bas Mach proposé dans [1] est bien adapté pour représenter la phase gazeuse. Dans ce modèle, la pression est décomposée en une pression thermodynamique uniforme en espace, et une pression hydrodynamique qui correspond à la pression usuelle des écoulements incompressibles. On introduit une fonction caractéristique $H(\vec{x}, t)$ égale à 1 dans le gaz et 0 dans le liquide. $H$ est simplement advectée par l'écoulement, et obéit donc à l'équation :

$$ \frac{dH}{dt} = \frac{\partial H}{\partial t} + \vec{v} \cdot \nabla H = 0 \quad (1) $$

où $\vec{v}$ est la vitesse du fluide. On considère un modèle diphasique à un fluide, qui exprime la conservation de la masse, de la quantité de mouvement et de l'énergie. On l'écrit :

$$ \left\{ \begin{aligned} \nabla \cdot \vec{v} &= -\frac{1}{\rho} \frac{d\rho}{dt} \\ \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \nabla \vec{v} + \frac{1}{\rho} \nabla \pi &= \frac{1}{\rho} \nabla \cdot \vec{\tau} \\ \frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T &= \frac{1}{\rho c_p} \nabla \cdot (k \nabla T) + H \frac{\gamma-1}{\gamma} T \frac{1}{P} \frac{dP}{dt} \end{aligned} \right. \quad (2) $$

où $P(t)$ est la pression thermodynamique, $\pi$ la pression hydrodynamique, $T$ la température, et $\vec{\tau}$ est le tenseur des contraintes visqueuses. La gravité et la tension superficielle ne sont pas prises en compte dans cette formulation. Le gaz est supposé obéir à une loi d'état de gaz parfait, soit $P(t) = R\rho(\vec{x},t)T(\vec{x},t)$ avec

$R = c_p - c_v$ et $\gamma = \frac{c_p}{c_v}$, $c_p$ et $c_v$ étant les chaleurs spécifiques. Le liquide est considéré incompressible de masse volumique constante $\rho_l$. L'équation d'état généralisée valable pour le liquide et le gaz s'écrit alors :

$\rho = H \frac{P}{RT} + (1-H)\rho_l$. Les propriétés physiques (conductivité thermique $k$, viscosité $\mu$, chaleurs spécifiques $c_p$ et $c_v$) des deux fluides sont obtenues de manière similaire. Enfin la pression thermodynamique $P(t)$, définie seulement dans le gaz, est obtenue par intégration de l'équation de continuité sur le domaine gazeux considéré $\Omega(t)$, ce qui donne :

$$ \frac{1}{P} \frac{dP}{dt} = \frac{1}{\int_{\Omega(t)} H dx} \int_{\Omega(t)} \left( H \frac{1}{T} \frac{dT}{dt} - \nabla \cdot \vec{v} \right) dx \quad (3) $$

L'intégration temporelle de cette relation fournit directement la pression :
$$P(t + \delta t) = P(t) \cdot \exp\left(\int_{t}^{t+\delta t} \frac{1}{P} \frac{dP}{d\xi} d\xi\right) \quad (4)$$

Lorsque le domaine gazeux est non simplement connexe (cas de plusieurs bulles par exemple), une pression thermodynamique est calculée dans chaque composante connexe, puis une pression thermodynamique « étendue » est définie dans tout l'espace afin de reproduire les effets dynamiques de ces différentes pressions sur le liquide.

## 2.2 Méthode numérique

Ce modèle est couplé à une méthode de suivi Lagrangien de l'interface dite « front-tracking », développée dans les travaux de Tryggvason et al [2] et suivants (Shin et al [3]). L'interface est discrétisée et déplacée de façon Lagrangienne. La fonction de Heaviside $H$, définie sur le maillage Eulérien, est construite à partir de la définition Lagrangienne. Cette fonction est lissée sur 3 à 4 mailles pour des raisons de stabilité numérique. Ceci conduit à l'existence d'une zone de « mélange » entre les deux fluides dans laquelle le modèle n'a pas de signification physique.

La discrétisation spatiale est basée sur des différences finies centrées d'ordre 2 pour les termes de convection et de diffusion. Un maillage décalé est utilisé, la pression hydrodynamique et la température étant localisées au centre des mailles, tandis que les composantes de la vitesse sont localisées sur les faces. Un schéma explicite d'ordre un en temps est utilisé. Un algorithme classique de type prédiction-projection permet de calculer les vitesses. La résolution de l'équation de Poisson pour la pression utilise un solveur multigrille spécialement adapté pour le cas de coefficients discontinus [4].

## 3 Résultats et discussion

Le modèle a été validé dans le cas unidimensionnel par comparaison avec des simulations ALE (non présenté ici pour raisons de place). On considère dans la suite un problème bidimensionnel. Deux bulles d'air sont incluses dans une cavité fermée remplie d'eau. La cavité, dont les parois sont adiabatiques, a pour dimensions 50µm x 100µm. A l'état initial, la température dans la cavité est uniforme, égale à 293,15 K. La première bulle a un rayon de 12µm et est à la pression thermodynamique de 1,4 atm, tandis que la seconde a un rayon de $1,4^{1/2}$ x 12µm et est à la pression de 1 atm. Les bulles sont centrées respectivement en (25µm x 25µm) et (25µm x 75µm). Les propriétés physiques des fluides sont celles de l'air et de l'eau, soit : pour l'eau, $\rho_l=1000$ kg/m³, $\mu_l=10^{-3}$ Pa.s, $c_l=4180$ Jkg⁻¹K⁻¹ et $k_l=0,6$ Wm⁻¹K⁻¹; pour l'air $\mu_v=1,82$ 10⁻⁵ Pa.s. $k_v=0,0256$ Wm⁻¹K⁻¹, $R=287$ Jkg⁻¹K⁻¹ et $\gamma=1,4$. La tension superficielle est prise en compte (tension superficielle de l'eau : 0,07 Nm⁻¹). Le maillage utilisé comporte 48 x 96 mailles.

L'évolution en temps des pressions thermodynamiques dans les deux bulles est représentée sur la figure 1 (à gauche). La figure de droite montre l'évolution temporelle du rapport entre leurs surfaces. La différence de pression initiale induit un mouvement oscillatoire de compression et dilatation alternées des bulles. La figure 2 montre la fonction caractéristique $H$ à plusieurs instants au cours du processus. On constate que les bulles sont soumises à des déformations de grande amplitude et à une variation de surface. Le mouvement s'amortit au cours du temps jusqu'à l'établissement d'un état stationnaire uniforme en température et pression.

Les valeurs à l'état final stationnaire peuvent être comparées aux valeurs théoriques fournies par la thermodynamique. On obtient pour la pression et la température finales :

$$P_f = \frac{P_{10} R_{10}^2 + P_{20} R_{20}^2}{R_{10}^2 + R_{20}^2}, \text{ soit } P_f = \frac{7}{6} \approx 1,167 \text{ atm, } T_f = T_0 = 293,15 \text{ K et } S_2/S_1 = 1 \text{ (rapport des surfaces des bulles). La comparaison avec les valeurs obtenues par la simulation numérique donne un excellent accord pour la pression thermodynamique et la température (erreur inférieure à 1‰). Le temps nécessaire à l'établissement de l'équilibre thermique est néanmoins beaucoup plus long que la durée du processus oscillatoire à la fin duquel les pressions thermodynamiques sont égales dans les deux bulles. Le rapport final des surfaces des bulles (figure 1 droite) montre un accord satisfaisant avec la valeur théorique. L'erreur, de l'ordre de 2%, est très probablement associée à la méthode de suivi d'interface. Cette erreur décroit lorsque l'on raffine le maillage, et converge à l'ordre 1. Pour vérifier les propriétés de conservation de la masse de la simulation, on a tracé sur la figure 3 l'évolution temporelle de la surface gazeuse totale dans la cavité,}$$
normalisée par la surface initiale. Ce rapport doit rester constant et égal à 1 puisque le liquide est incompressible et la cavité rigide. On constate que l'erreur reste très faible tout au long de l'évolution, ce qui démontre les très bonnes propriétés de conservation de l'algorithme et de la méthode numérique.

Les plus grandes vitesses observées dans l'écoulement (de l'ordre de 1 m/s) le sont aux premiers instants, lorsque les différences de pression thermodynamiques sont les plus grandes. Ceci valide notre hypothèse d'écoulement bas Mach dans le gaz sur laquelle est basé le modèle.

<figure><img src="image_1.png" /><figcaption>FIG. 1 – Evolution temporelle de la pression thermodynamique dans les bulles (à gauche) et du rapport de surface des bulles (à droite). L'unité de pression est 1 atm.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>FIG. 2 – Représentation de la fonction H montrant les compressions et dilatations des deux bulles d'air (en rouge) dans le liquide incompressible (en bleu). Les temps correspondants sont, de gauche à droite :</figcaption></figure>

$$t = 0, t = 4\mu s = \tau_{osc} / 2, t = 8\mu s = \tau_{osc}, \text{ état stationnaire.}$$

<figure><img src="image_3.png" /><figcaption>FIG. 3 – Evolution temporelle de la surface gazeuse totale, normalisée par la surface initiale.</figcaption></figure>

Les oscillations observées durant le processus instationnaire sont caractérisées par leur fréquence et temps d'amortissement. Du point de vue énergétique, un cas adiabatique idéal correspondrait à une situation sans échange de chaleur entre les bulles et leur environnement, ie la phase liquide. Ceci se rencontrait avec un
liquide de conductivité thermique nulle. Comme le processus oscillatoire s'effectue sans perte thermique, l'énergie initialement disponible est principalement transférée d'une bulle à l'autre sous forme de travail de compression/dilatation. L'amortissement est uniquement dû aux effets de dissipation visqueuse dans le liquide ou aux parois. Si ces effets sont négligeables, le système subit un processus réversible et le temps caractéristique d'amortissement tend vers 0. Le problème réel concerne un liquide de conductivité thermique finie. Pendant le processus transitoire, une partie de l'énergie disponible est convertie en chaleur et transmise au liquide. En retour ceci réduit la quantité d'énergie mécanique pouvant être échangée entre les deux bulles d'air. Les oscillations sont graduellement amorties et le système subit un processus irréversible.

Des simulations ont été effectuées pour examiner trois situations, concernant le même système avec les mêmes propriétés physiques sauf pour la conductivité thermique du liquide. La première correspond au cas présenté sur la figure 1 et sera désigné dans la suite comme le cas « nominal ». Dans la seconde (resp. troisième) la conductivité thermique est prise égale à 2 Wm⁻¹K⁻¹ (0,02 Wm⁻¹K⁻¹). La comparaison est effectuée sur l'évolution de la pression dans la première bulle. Les courbes de la figure 4 montrent que la fréquence d'oscillation des bulles décroit lorsque la conductivité thermique du liquide croit. Les valeurs correspondantes sont regroupées dans le tableau 1. On note sur cette figure que la première oscillation est presque identique pour les trois cas. Un décalage temporel apparaît par la suite lorsque la chaleur diffuse au voisinage des interfaces. Les valeurs des fréquences peuvent être comparées avec un modèle analytique développé dans [5] sur un cas unidimensionnel analogue. Ce modèle prédit que la fréquence naturelle d'un tel système perturbé autour d'un état d'équilibre à pression P₀ est limitée par une valeur supérieure fₐ qui correspond à un processus adiabatique idéal :

$$f_a = \frac{1}{2\pi} \sqrt{\gamma \frac{P_0 L_1 + L_2}{m_l L_1 L_2}} \quad \text{où } m_l \text{ est la masse de liquide et } L_1 \text{ et } L_2 \text{ sont les longueurs des zones gazeuses.}$$

Pour effectuer une comparaison avec le cas bidimensionnel, on peut réaliser une analogie avec la situation 1D en considérant que les longueurs correspondantes sont les diamètres des bulles, la masse de liquide étant proportionnelle à la distance entre les bulles, soit ici 24µm. Nous prenons la pression P₀ comme celle atteinte à l'état stationnaire, soit P₀ = 7/6 atm. La fréquence adiabatique obtenue est fₐ = 116 kHz, qui est proche des valeurs du tableau 1. Une différence de l'ordre de 10% correspond aux valeurs non nulles de la conductivité thermique et aux effets multidimensionnels.

On discute maintenant l'effet de la conductivité thermique sur le temps d'amortissement. En considérant les courbes de la figure 4 nous choisissons de représenter l'amortissement de la pression par la fonction suivante :

$$P_1 = P_{10} \exp(-t / \tau) \cos(\omega t) \quad (5)$$

où ω = 2πf et τ désigne le temps d'amortissement. Les valeurs obtenues pour les trois cas sont listées dans le tableau 2. Elles confirment que l'évolution de ce temps caractéristique est monotone en fonction de la conductivité thermique.

Le champ de température dans la cavité au temps t = 12µs est représenté sur la figure 5 pour les trois cas. On constate que la conductivité thermique du liquide affecte fortement le champ dans les bulles. Bien que les différences de température entre le gaz et le liquide induisent un échange de chaleur au travers des interfaces, la température reste inchangée dans le liquide du fait de sa résistance thermique élevée. La continuité du flux thermique aux interfaces induit un champ de température différent dans les bulles suivant la valeur de kᵢ. Dans le cas où kᵢ est faible (kᵢ=0,02 Wm⁻¹K⁻¹), le gaz et le liquide ont des conductivités thermiques du même ordre, produisant un gradient de température presque nul dans les bulles, où la température est alors presque uniforme, inférieure ou supérieure à celle du liquide suivant qu'on se trouve en étape de compression ou de dilatation. Au contraire, quand kᵢ est grand (kᵢ=2 Wm⁻¹K⁻¹), le gradient de température devient significatif du côté gazeux des interfaces. Les bulles ayant quasiment les mêmes surfaces dans les trois cas, la loi d'état du gaz implique que la pression thermodynamique dans chaque bulle est proportionnelle à sa température moyenne. La différence de pression entre les deux bulles pilotant le processus oscillatoire, toute augmentation de la conductivité thermique du liquide accélère l'amortissement.
## References

[1] Chenoweth, D., Paolucci, S. 1986. Natural convection in an enclosed vertical air layer with large horizontal temperature differences. J. of Fluid Mechanics 169, 173-210.

[2] Tryggvason, G. et al. 2001. A front-tracking method for the computations of multiphase flow. J. of Computational Physics 169, 708-759.

[3] Shin, S., Juric, D. 2002. Modeling three-dimensional multiphase flows using a level contour reconstruction method for front-tracking without connectivity. J. of Computational Physics, 180, 427-470.

[4] Kwak, D.Y. and Lee, J.S., 2004. Multigrid algorithm for the cell-centered finite difference method II: discontinuous coefficient case. Num. methods for partial differential equations, 20 (5), 742-764.

[5] Duluc, M.-C., Le Maître O., Daru V., Le Quéré P., 2009. Numerical study of liquid inclusion oscillations inside a closed 1-d microchannel filled with gas. Microfluidics & Nanofluidics 6, 163-177.

<figure><img src="image_4.png" /><figcaption>FIG. 4 – Evolution temporelle de la pression thermodynamique dans la bulle 1 pour 3 valeurs de la conductivité thermique du liquide.</figcaption></figure>

<table><caption>TABLEAU 1 – Fréquence des oscillations pour 3 valeurs de la conductivité thermique du liquide</caption><tr><td>k<sub>l</sub> (Wm<sup>-1</sup>K<sup>-1</sup>)</td><td>0.02</td><td>0.6</td><td>2</td></tr><tr><td>f (kHz)</td><td>127.6</td><td>126</td><td>121.5</td></tr></table>

<table><caption>TABLEAU 2 – Temps caractéristique d'amortissement pour 3 valeurs de la conductivité thermique du liquide</caption><tr><td>k<sub>l</sub> (Wm<sup>-1</sup>K<sup>-1</sup>)</td><td>0.02</td><td>0.6</td><td>2</td></tr><tr><td>τ (s)</td><td>5.7 10<sup>-5</sup></td><td>2.1 10<sup>-5</sup></td><td>1.9 10<sup>-5</sup></td></tr></table>

<figure><img src="image_5.png" /><figcaption>FIG. 5 – Champ de température dans la cavité aux temps $t = 12 \mu s \approx 1,5 \tau_{osc}$. $k_l = 0,02 Wm^{-1}K^{-1}$ (gauche), $k_l = 0,6 Wm^{-1}K^{-1}$ (milieu), $k_l = 2 Wm^{-1}K^{-1}$ (droite). Les bulles sont matérialisées par le contour noir.</figcaption></figure>