<figure><img src="image_1.png" /></figure>

HAL
open science

# Développement d'un solveur tout-couplé parallèle 3D pour la simulation des écoulements diphasiques incompressibles

Mohamed El Ouafa, Stéphane Vincent, Vincent Le Chenadec, Benoît Trouette

► To cite this version:

Mohamed El Ouafa, Stéphane Vincent, Vincent Le Chenadec, Benoît Trouette. Développement d'un solveur tout-couplé parallèle 3D pour la simulation des écoulements diphasiques incompressibles. 25ème Congrès Français de Mécanique (CFM 2022), Aug 2022, Nantes, France. hal-04079918

HAL Id: hal-04079918

https://hal.science/hal-04079918

Submitted on 25 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Développement d'un solveur tout-couplé parallèle 3D pour la simulation des écoulements diphasiques incompressibles

Conference Paper · March 2023

CITATIONS

0

READS

63

4 authors:

<figure><img src="image_2.png" /></figure>

Elouafa Mohamed

Université Gustave Eiffel

7 PUBLICATIONS 10 CITATIONS

SEE PROFILE

<figure><img src="image_3.png" /></figure>

Vincent Le Chenadec

CentraleSupélec

30 PUBLICATIONS 374 CITATIONS

SEE PROFILE

<figure><img src="image_4.png" /></figure>

Stéphane Vincent

Université Gustave Eiffel

260 PUBLICATIONS 3,106 CITATIONS

SEE PROFILE

<figure><img src="image_5.png" /></figure>

Benoît Trouette

Université Gustave Eiffel

36 PUBLICATIONS 112 CITATIONS

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

<figure><img src="image_6.png" /></figure>

Block Preconditioning of Monolithic Solvers for Incompressible Two-Phase Flows at Large Contrast Ratios View project

<figure><img src="image_7.png" /></figure>

Multi-scale modeling and hybrid simulations View project
# Développement d'un solveur tout-couplé parallèle 3D pour la simulation des écoulements diphasiques incompressibles

M. El Ouafa, S. Vincent, V. Le Chenadec, B. Trouette

Univ Gustave Eiffel, Univ Paris Est Creteil, CNRS, UMR 8208, MSME, F-77454 Marne-la-Vallée, France

email : mohamed.elouafa@univ-eiffel.fr

## Résumé

Les équations de Navier-Stokes diphasiques jouent un rôle crucial dans plusieurs applications industrielles et environnementales. Cependant, leur résolution pose de nombreux problèmes et devient un enjeu important pour les mathématiciens ansi que les ingénieurs. Dans le contexte des écoulements diphasiques incompressibles, les difficultés sont multiples : d'une part, le couplage vitesse-pression induit par la contrainte d'incompressibilité, donne lieu à un problème de point selle pour lequel des solveurs dédiés doivent être développés. Et, d'autre part, les forts rapports de densité et de viscosité en présence de grandes déformations de l'interface conduisent à des systèmes linéaires très mal conditionnés. Nous décrivons ici un solveur tout-couplé parallèle 3D, didié à la simulation des écoulements diphasiques incompressibles. Les équations linéarisées de quantité de mouvement et de continuité découlant de la solution implicite des vitesses et de la pression dans les fluides sont résolues simultanément. Nous avons simulé un cas de l'instabilité d'un jet visqueux pour vérifier les modèles, notamment en présence de forts rapports de densité et de viscosité. Nous avons également réalisé un cas d'inversion de phase 3D, pour démontrer la capacité de notre nouveau solveur tout-couplé à résoudre des problèmes diphasiques sur un milliards de cellules.

## Abstract :

The Navier-Stokes equations play a crucial role in several industrial and environmental applications. However, their resolution causes many problems and remains a challenging task for mathematicians as well as engineers. In the context of incompressible two-phase flows, the difficulties are multiple : on the one hand, the velocity-pressure coupling induced by the incompressibility constraint, gives rise to a saddle point problem for which suitable solvers must be developed. And, on the other hand, the high density and viscosity ratios in the presence of large interface deformations lead to very ill-conditioned linear systems. We describe here a fully-coupled 3D parallel solver for the simulation of two-phase incompressible flows. The linearized momentum and continuity equations arising from the implicit solving of the fluid velocities and pressure are solved simultaneously. The instability of a 3D viscous jet is simulated, in order to validate the models, especially in the presence of strong density and viscosity ratios. We also performed a 3D phase inversion case, in order to exhibit the capacity of our new fully-coupled solver to deal with two-phase problems on a billion cells.

**Mots clefs :** Solveur tout-couplé, point-selle, 3D parallèle, écoulements diphasiques, fort ratios
# 1 Introduction

Les écoulements diphasiques incompressibles mettant en œuvre un liquide et un gaz sont omniprésents dans de nombreuses applications industrielles et académiques et, par conséquent, ils font l'objet de nombreuses études numériques. Les écoulements d'impact de gouttes sur des surfaces, la fabrication de dépots par jet de plasma thermique, le déferlement d'une vague sur une plage, l'atomisation des jets ou les sprays rencontrés lors du refroidissement de parois de bâtiments ou dans les moteurs sont, par exemple, des problèmes diphasiques et turbulents.

La modélisation et la simulation numérique de ce type d'écoulement reste un grand défi en raison, entre autres des aspects mutli-échelles (allant de la phase séparée à l'échelle des petites interfaces) et multi-physiques (phases séparées ou dispersées, frottements, effets capillaires, changement de phase, mouillabilité...). En effet, les forts rapports de densité et de viscosité dans les écoulements diphasiques engendrent souvent un cisaillement important et des discontinuités à l'interface fluide-fluide. Ces discontinuités compliquent particulièrement la discrétisation des équations de conservation, et conduisent à des instabilités numériques. En outre, le traitement du couplage vitesse-pression dans les écoulements monophasiques incompressibles a fait l'objet de nombreux travaux et est relativement bien maîtrisé. Au contraire, dans le cadre d'écoulements diphasiques avec de larges rapports de densité et de viscosité, des difficultés numériques apparaissent, liées à la présence de fortes déformations de l'interface et à des sauts de propriétés des fluides matériels qui y sont associés. Elles conduisent à des systèmes linéaires mal conditionnés, ce qui rend le problème beaucoup plus difficile à résoudre. Au niveau discret, le problème persiste sous la forme de recherche d'un point selle dont la résolution est couteuse. Pour faire face à ce problème, de gros progrès ont cependant été réalisés au cours des dernières décennies, et de nombreux méthodes ont été proposées, elles peuvent être classés ainsi en deux catégories. D'une part, les méthodes approchées, introduites par Chorin et Temam [6, 21] et leurs variantes [10, 4, 3, 11], approximent le système original en utilisant le splitting en temps, ce qui aboutit à deux équations découplées : l'une pour mettre à jour le champ de vitesse et l'autre le champ de pression. Cela soulève un certain nombre de problèmes, tels que l'apparition de couches numériques en raison de splitting en temps, ou la nécessité de conditions limites de pression qui n'existent pas dans le problème original. D'autre part, les méthodes exactes résolvent les deux champs (vietsse et pression) simultanément, préservant ainsi la consistance du système discrétisé avec les équations continues. Dans ces méthodes exactes, le système à point selle peut être résolu soit directement par un solveur tout-couplé [16, 1, 18, 17], soit par une méthode de Lagrangian augmenté (AL) [9, 23]. L'intérêt de ces méthodes réside dans le fait qu'elles permettent une résolution très précise et robuste, sans erreur de splitting (sauf celle liée à la linéarisation du terme d'inertie) des écoulements diphasiques incompressibles.

Nous développons ici une méthode exacte 3D parallèle, tout-couplée, pour simuler des écoulements supposés incompressibles et diphasiques. Nous décrivons d'abord le modèle, les schémas de discretisations et les solveurs seront présentées dans un seconde temps. Nous validons ensuite notre approche numérique de VOF lorsqu'elle est appliquée à la résolution de l'équation d'advection de la fraction volumique. Ensuite, Nous considérerons un cas 3D d'inversion de phase, pour montrer la capacité de la méthode à traiter des cas diphasiques sur plus d'un milliard de cellules. Enfin, un cas 3D d'un jet visqueux sera présenté à fin de mettre en valeur l'intérêt. du solveur tout-couplé à simuler des problémes raides à fort rapport de masse volumique et de viscosité.
# 2 Méthodes numériques

## 2.1 Approche diphasique de la résolution des équation de Navier-Stokes

Dans ce travail, l'écoulement diphasique est décrit par le modèle à un fluide de Katoaka [14], basé sur une approche eulérienne où les phases sont vues via des fractions de volume. Dans chacune des deux phases fluides, les équations à résoudre sont : l'équation de conservation de la masse et l'équation de conservation de la quantité de mouvement, avec un terme source additionnel dans l'équation de la quantité de mouvement (la force capillaire) et une équation d'advection supplémentaire pour la fonction de phase C. La fonction C, qui peut prendre des valeurs comprises entre 0 et 1 représente le taux de présence d'un des deux fluides dans un volume de contrôle. Les équations de modèle à un fluide dans le cas diphasique s'écritent ainsi :

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \overline{\overline{\mathbf{B}}}_{\text{penu}} (f(\mathbf{u}) - \mathbf{u}_{\infty}) = -\nabla p + \nabla \cdot \overline{\overline{\mathbf{T}}} + \rho \mathbf{g} + \mathbf{F}_s, \quad (1) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (2) $$

$$ \frac{\partial C}{\partial t} + \nabla \cdot (\mathbf{u}C) = C\nabla \cdot \mathbf{u}, \quad (3) $$

$$ \rho(C) = \rho_1 C + (1 - C)\rho_2, \quad (4) $$

$$ \mu(C) = \mu_1 C + (1 - C)\mu_2 \quad (5) $$

Avec **u** le champ de vitesse, **p** le champ de pression, **t** le temps, $\rho$ et $\mu$ la masse volumique et la viscosité du fluide de mélange, données par une moyenne arithmétique, avec $\rho_1$, $\rho_2$, $\mu_1$ et $\mu_2$ les masses volumiques et les viscosités des fluides 1 et 2 respectivement. En outre, $\overline{\overline{\mathbf{T}}}$ représente le tenseur des taux de déformations pour un fluide Newtonien [5], tandis que $\overline{\overline{\mathbf{B}}}_{\text{penu}}$ est un tenseur d'ordre 2 pour spécifier les conditions aux limites [19], dont les composantes diagonales tendent vers l'infini le long de la frontière et sont identiquement nulles à l'intérieur du domaine de fluide. En outre, $\mathbf{F}_s$ est le terme capillaire agissant sur l'interface, modélisé dans cette étude par le modèle CSF de Brackbill [2].

## 2.2 Schémas de discrétisation et solveurs

Dans le présent article, les équations du modèle diphasique incompressible et instationnaire sont approximées via des volumes finis pénalisés (pour traiter différents types de conditions aux limites), sur un maillage décalé. Les variables discrètes de vitesse et de pression sont disposées de manière décalée, chacune ayant son propre volume de contrôle, comme le montre la Figure 1. Ce choix garantit la concordance des opérateurs différentiels tels que la divergence et le gradient, et il évite également les oscillations sur le champ de pression. Sachant que la solution de l'équation d'advection (3) est effectuée avec la méthode VOF conservative de Weymouth [12], qui conserve le volume à la précision machine. Ainsi, Le système semi-discret résultant de la discretisations des équations de conservation de la masse (Eq. 2) et de la quantité de mouvement est avancé dans le temps, en utilisant un schéma implicite De Gear du second ordre. En revanche, un schéma centré pour le gradient de pression, la divergence de vitesse, les termes convectifs et visqueux est utilisé. Pour éviter de mettre en œuvre la linéarisation de Newton et d'inverser le Jacobien complet, qui s'avère plus efficace mais plus coûteux pour les écoulements installationnaires, les termes convectifs sont linéarisés en utilisant l'extrapolation d'Adams-Bashforth comme
suit :

$$ \mathbf{u}^{n+1} \otimes \mathbf{u}^{n+1} \approx (2\mathbf{u}^n - \mathbf{u}^{n-1}) \otimes \mathbf{u}^{n+1}. \qquad (6) $$

<figure><img src="image_8.png" /><figcaption>FIGURE 1 – Grille décalée de vitesse et pression : représentation des nœuds supplémentaires de vitesse en dehors de la frontière physique.</figcaption></figure>

Ici, $n$ désigne un indice de temps associé au temps physique $n\Delta t$ et $\Delta t$ est supposé être le pas de temps constant à chaque itération résolue. Ainsi, après une linéarisation et une discrétisation appropriées, et, comme une résolution implicite tout-couplé (Fully Coupled) est visée, le système générique est un système linéaire non symétrique de type point selle. Les deux inconnues du problème, i.e. la vitesse $\mathbf{u}$ et la pression $p$, sont couplées par la contrainte d'incompressibilité $\nabla \cdot \mathbf{u} = 0$. Ce système de points selle peut être reformulé comme suit :

$$ \left( \begin{array}{c} \frac{1}{\Delta t} M_u^{(\rho)} + N_u^{(\rho)} + L_u^{(\mu)} \\ B_u \end{array} \quad \begin{array}{c} B_p^T \\ 0 \end{array} \right) \left( \begin{array}{c} \mathbf{u}^{n+1} \\ p^{n+1} \end{array} \right) = \left( \begin{array}{c} \mathbf{f} \\ 0 \end{array} \right) \qquad (7) $$

Où la matrice $B_u$ est la partie discrète de la divergence et $B_p^T$ représente l'opérateur discret du gradient de pression. En outre, $M_u^{(\rho)}$ est la matrice de masse de la vitesse, $L_u^{(\mu)}$ est le Laplacien discret de la vitesse, $N_u^{(\rho)}$ désigne la matrice associée au terme convectif et $\mathbf{f}$ représente le terme source. Le système réduit peut être réécrit sous la forme $Ax = b$ avec $A$ la matrice associée aux coefficients de discrétisation spatio-temporelle et $b$ le second membre associé à chaque composante du vecteur solution. La solution du problème diphasique est alors :

$$ x = A^{-1} \times b \qquad (8) $$

Notre stratégie pour résoudre ce système linéaire creux est de faire appel au solveur BiCGStab(2) [13] sur la totalité du système, en choisissant la structure CSR (compressed storage raw) afin de ne stocker que les coefficients non nuls de la matrice, et ensuite de construire un préconditionneur efficace pour accélérer la convergence du solveur itératif. Cette stratégie peut être mise en place en utilisant le préconditionneur pression-convection-diffusion (PCD) [1] pour le bloc pression, et le préconditionneur triangulaire par
bloc pour le bloc de vitesse (nous invitions le lecteur à consulter [17, 18] pour plus de détails).

# 3 Résultats

## 3.1 Cas test pour l'équation d'advection

Afin de montrer les capacités de notre méthode VOF à décrire correctement l'interface, même lorsque celle-ci est fortement étirée et déformée, un cas test standard de transport d'interface est étudié. Ce cas test a été proposé par Enright et al. [7] à partir des travaux de Leveque [15]. Initiallement, une goutte sphérique de rayon $r = 0.15$ est placé aux coordonnées $(0.35, 0.35, 0.35)$ dans une boite periodique de coté $L = 1$. Dans cette configuration, l'interface est transportée avec le champ de vitesse suivant :

$$ 
\begin{align} 
u(x, y, z) &= 2 \sin^2(\pi x) \sin(2\pi y) \sin(2\pi z) \nonumber \\ 
v(x, y, z) &= -\sin(2\pi x) \sin^2(\pi y) \sin(2\pi z) \tag{9} 
\end{align} 
$$

$$ w(x, y, z) = -\sin(2\pi x) \sin(2\pi y) \sin^2(\pi z) $$

Ce champs de vitesse est modulé par une fonction de temps $f(t)$ définie par $f(t) = \cos(\pi t/T)$, avec $T = 3s$ la période. Ainsi, le champs de vitesse est réversible, en d'autres termes, la goutte revient à sa position initiale en $t = T$.

La figure 2 présente les résultats que nous avons obtenus pour une interface sphérique avec un maillage de 512³. Nous pouvons voir que nous arrivons à des interfaces qui sont similaires à celles proposées dans les travaux de Douglas Enright. Nous avons également réalisé une étude de convergence en calculant l'ordre de précision spatial. Pour ce faire, nous avons choisi deux types d'erreur : d'une part, l'erreur de conservation de masse intégrée sur le domaine de calcul. Elle est définie comme suit,

$$ \epsilon_{mass} = \frac{\int_V c(\mathbf{x}, t) \, dV - \int_V c(\mathbf{x}, 0) \, dV}{\int_V c(\mathbf{x}, 0) \, dV} \qquad (10) $$

Et d'autre part, l'erreur géometrique qui mesure la différence entre l'interface calculée et l'interface exacte :

$$ \epsilon_{geo} = \int_V [c(\mathbf{x}, T) - c_{exacte}(\mathbf{x}, T)] \, dV \qquad (11) $$

Nous dressons les résultats de l'étude de convergence dans le tableau 1. Nous constatons une conservation de la masse à la précision machine quelque soit le maillage. En outre, un ordre 2 est obtenu sur l'erreur géometrique ce qui valide notre méthode de transport.

<table><caption>TABLE 1 – Etude de convergence schéma VOF de Weymouth : Goutte sphérique dans un champs de vitesse réversible.</caption><thead><tr><th>Maillage</th><th>ε<sub>mass</sub></th><th>ordre</th><th>ε<sub>geo</sub></th><th>ordre</th></tr></thead><tbody><tr><td>32<sup>3</sup></td><td>2.05×10<sup>-12</sup></td><td>-</td><td>7.82×10<sup>-3</sup></td><td>-</td></tr><tr><td>64<sup>3</sup></td><td>4.69×10<sup>-13</sup></td><td>-</td><td>2.72×10<sup>-3</sup></td><td>1.52</td></tr><tr><td>128<sup>3</sup></td><td>5.72×10<sup>-14</sup></td><td>-</td><td>6.50×10<sup>-4</sup></td><td>2.07</td></tr><tr><td>256<sup>3</sup></td><td>2.18×10<sup>-14</sup></td><td>-</td><td>1.51×10<sup>-4</sup></td><td>2.11</td></tr></tbody></table>
<figure><img src="image_9.png" /><figcaption>FIGURE 2 – Maillage 512³ : Évolution de l'interface sphérique dans un champ de vitesse réversible à divergence nulle</figcaption></figure>

## 3.2 Cas d'application : L'inversion de phase eau-huile

Un écoulement diphasique, eau-huile par exemple, dans lequel une des phases est dispersée dans une autre, peut correspondre à une inversion de phase, c'est-à-dire, la phase dispersée devenant la phase continue. Ce phénomène instationnaire est ensuite accompagné de plusieurs processus interfaciaux à petite et grande échelle, tels que la déformation, la formation de ligaments, la rupture de l'interface et la coalescence (voir figure 3b).

Ici, nous menons une simulation numérique directe sur un cas d'inversion de phase étudié récement par Saeedipour et al [20], afin de montrer la capacité de notre solveur tout-couplé à faire des simulations massivement parallèle. Les quantités d'intérêt seront donc les grandeurs macrosopiques, telles que l'énergie potentielle, l'énergie cinétique, l'enstrophie, la dissipation visqueuse, l'aire interfaciale et le volume du fluide léger dans la partie supérieure du domaine. Dans ce qui suit, nous présentons la géometrie et les paramètres physiques du problème. Ensuite, nous discutons les résultats des simulations
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 3 – Problème de l'inversion de phase liquide/liquide dans une boîte fermée [8]</figcaption></figure>

numériques ainsi que leurs comparaisons avec les références.

### 3.2.1 Géometrie et paramètres physiques

Comme décrit sur la figure 3a, un volume cubique d'un fluide léger (identifié comme le fluide 2) est placé dans la partie inférieure d'une boîte remplie d'un fluide lourd (le fluide 1). Nous pourrons assimiler les fluides 1 et 2 comme étant respectivement de l'eau et de l'huile. Nous imposerons une viscosité égale pour les 2 fluides ($\mu_1 = \mu_2 = \mu$). La boîte est de dimension ($H, H, H$) et la goutte de fluide 1 ($H/2, H/2, H/2$). Nous considérerons un glissement libre et un mouillage partiel sur les murs extérieurs. Nous considérerons également que les angles de contact de 90° sont équivalents à la symétrie d'un miroir. L'accélération gravitationnelle est choisie comme étant $g = (0,0,-9.81)$. Les différents paramètres du cas d'étude sont donnés dans le tableau 2. Les nombres adimensionnels qui gouvernent la dynamique de ce problème sont donc le nombre de Reynolds (Re) et le nombre de Weber (We) définis de la manière suivante :

$$Re = \frac{\rho_1 H U_g}{\mu_1} \qquad (12)$$

$$We = \frac{\rho_1 H U_g^2}{\sigma} \qquad (13)$$

Dans cette étude, nous utiliserons la définition suivante pour la vitesse caractéristique :

$$U_g = \frac{\rho_2 - \rho_1}{\rho_1} \sqrt{\frac{gH}{2}} \qquad (14)$$

Ainsi, nous définissons l'énergie cinétique et l'énergie potentielle pour la phase n de la manière suivante :

$$Ek_n = \frac{1}{2} \int_{\Omega} C_n \rho_n \mathbf{u}^2 dV \qquad (15)$$

$$Ep_n = \int_{\Omega} C_n \rho_n gz dV \qquad (16)$$
avec $C_n$ la fonction couleur définie telle que $C_n = 1$ dans le fluide n et $C_n = 0$ ailleurs. Nous définissons également $C_n$ la fonction couleur de la phase n : $C_1 = C$ et $C_2 = 1 - C$. L'évolution en temps de l'intégrale de l'enstrophie dans les deux fluides sera également déterminée. Elle sera définie par :

$$Er_n = \int_{\Omega} C_n (\nabla \times \mathbf{u})^2 dV \qquad (17)$$

En ce qui concerne la dissipation visqueuse de la phase n, elle est définie comme suit :

$$Diss_n = \int_{\Omega} C_n 2\mu_n \bar{\bar{S}} : \bar{\bar{S}} dV \qquad (18)$$

Avec $\bar{\bar{S}} = 1/2 (\nabla \mathbf{u} + (\nabla \mathbf{u})^T)$. Quant à l'aire interfaciale, et contrairement aux études précédentes [20, 8], où elle est approximée comme étant $\sum_{int} S_{cell}$ où $S_{cell} = \Delta x \Delta y = \Delta y \Delta z = \Delta x \Delta z$, dans cette étude, elle est évaluée à partir de la normale à l'interface $\mathbf{n}$, en calculant dans un cube unitaire, la surface coupée par le plan $\mathbf{n} \cdot \mathbf{x} = \alpha$. Enfin, la grandeur $Vol$ représente le volume de fluide léger, elle permet de vérifier que le volume de fluide léger reste constant au cours du temps.

Concernant les paramètres numériques, les simulations ont été réalisées sur une grille cartésienne allant de $128^3$ à $1024^3$ cellules, avec un résidu de $\epsilon = 10^{-5}$ du solveur BiCGStab(2). En ce qui concerne les dérivées temporelles, un pas de temps constant $\Delta t = 10^{-3} s$ a été choisi, et 10000 incréments temporels seront ainsi calculés, correspondant à 10 secondes de simulation de l'écoulement. Pour le cas $1024^3$, il faut noter que les simulations ont été réalisées sur 6696 processeurs du cluster HPC Occigen du CINES et ont pris environ 5 jours de temps de calcul.

<table><caption>TABLE 2 – Propriétés physiques des fluides et nombres adimensionnels</caption><thead><tr><th>ρ<sub>1</sub></th><th>ρ<sub>2</sub></th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>g</th><th>σ</th><th>H</th><th>Re</th><th>We</th></tr><tr><th>(kg.m<sup>-3</sup>)</th><th>(kg.m<sup>-3</sup>)</th><th>(Pa.s)</th><th>(Pa.s)</th><th>(m.s<sup>-2</sup>)</th><th>(kg.s<sup>-2</sup>)</th><th>(m)</th><th></th><th></th></tr></thead><tbody><tr><td>1000</td><td>900</td><td>1</td><td>1</td><td>9.81</td><td>0.45</td><td>1</td><td>221</td><td>121</td></tr></tbody></table>

### 3.2.2 Résultats

Dans nos simulations numériques présentées sur la figure 4, nous observons les mêmes phénomènes que ceux présentés dans les travaux de Saeedipour et al. [20]. Tout d'abord la montée du fluide léger vers le haut de la boite. Ensuite, le fluide léger recircule avec plus au moins de fragmentation et dispersion. Pour l'ensemble des grandeurs que nous avons calculées, nous remarquons que les convergences en maillage par rapport aux références sont atteintes à partir d'un nombre de cellules égal à $256^3$ (Figures 5, 6, 7, 8, 9, 10, 11, 12, 13, 14), à l'exception de l'aire interfaciale. Nous avons également remarqué que les courbes présentant des pics, comme l'énergie cinétique, l'énergie potentielle, l'enstro,phie et la dissipation, retrouvent ces pics aux mêmes instants que la référence. Enfin, nous constatons que le volume de fluide léger $Vol$ reste constant tout au long de la simulation, quelque soit le maillage.
<figure><img src="image_12.png" /></figure>

FIGURE 4 – Maillage 256³ : l'inversion de phase liquide/liquide au cours du temps. L'interface des deux fluides est représentée en jaune.

<figure><img src="image_13.png" /><figcaption>FIGURE 5 – $EK_1$ : l'énergie cinétique de l'eau</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 6 – $EK_2$ : l'énergie cinétique de l'huile</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 7 – $EP_1$ : l'énergie potentielle de l'eau</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 8 – $EP_2$ : l'énergie potentielle de l'huile</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 9 – $ER_1$ : l'enstrophie de l'eau</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 10 – $ER_2$ : l'enstrophie de l'huile</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 11 – $Dis_1$ : la dissipation de l'eau</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>FIGURE 12 – $Dis_2$ : la dissipation de l'huile</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 13 – L'aire interfaciale</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 14 – Vol : le volume de l'huile</figcaption></figure>

## 3.3 Caractérisation de l'instabilité 3D d'un jet visqueux à faible nombre de Reynolds

Le problème de l'instabilité d'un jet est un phénomène crucial, dans la mesure où il est impliqué dans de nombreuses applications industrielles, telles que l'industrie alimentaire ou les applications aéronautiques et aérospatiales, avec le remplissage des propulseurs de fusée avec du propergol. Les conditions d'écoulement qui y correspondent sont à faible nombre de Reynolds. Initialement étudié par Tomé et al. [22], qui ont permis entre autre de caractériser le développement des instabilités du jet en fonction du nombre de Reynolds Re et du rapport $H/D$ entre la hauteur de la cavité et le diamètre d'injection, Peu d'études numeriques sont consacrées à l'étude de l'injection d'un fluide visqueux en régime de Stokes dans une cavite. Dans ce contexte, l'intérêt de ce cas est double : d'une part, la possibilité de confronter
les solutions numériques aux solutions expérimentales. D'autre part, les forts gradients de viscosité et de densités près de l'interface engendrent des systèmes linéaires très mal conditionnés, ce qui représente un véritable défi au niveau algébrique.

### 3.3.1 Géometrie et paramètres physiques

Le probléme étudié est le suivant : une cavité cubique de coté $H$, remplie d'air est étudiée. Dans sa partie supérieure, on injecte un fluide très visqueux à un debit constant, tandis que la la force gravitationnelle agit avec une intensité constante $g$ dans le sens négatif de la direction $z$. Pour les bords du domaine, nous avons utilisé des conditions limites de glissement sur les parois solides, et des conditions limites de de sortie, dans la partie superieure de la cavité. Les caractéristiques de la simulation numérique sont données dans le tableau 3. En outre, Tomé et al [22] donnent deux paramètres déterminants pour l'ap-

<table><thead><tr><th>Test case</th><th>ρ<sub>1</sub></th><th>ρ<sub>2</sub></th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>Re</th><th>H/D</th><th>V<sub>0</sub></th><th>ρ<sub>2</sub>/ρ<sub>1</sub></th><th>μ<sub>2</sub>/μ<sub>1</sub></th></tr></thead><tbody><tr><td>1</td><td>1.1768</td><td>1800</td><td>1.85 × 10<sup>-3</sup></td><td>500</td><td>0.1728</td><td>12.5</td><td>0.6</td><td>1529</td><td>1.6 × 10<sup>5</sup></td></tr><tr><td>2</td><td>1.1768</td><td>1800</td><td>1.85 × 10<sup>-3</sup></td><td>300</td><td>0.6</td><td>10</td><td>1</td><td>1529</td><td>2.7 × 10<sup>5</sup></td></tr></tbody></table>

TABLE 3 – Les paramètres physiques décrivant le cas du jet visqueux. Ici, $ρ_2$ et $μ_2$ représentent respectivement la masse volumique et la viscosité du propergol, $ρ_1$ et $μ_1$ représentent respectivement la masse volumique et la viscosité du gaz, $Re = \frac{ρ_2 D V_0}{μ_2}$ est le nombre de Reynolds associé au propergol, $\frac{H}{D}$ est le rapport d'aspect du jet et $V_0$ est la vitesse d'entrée. - pour des raisons de simplicité, le propergol est supposé se comporter comme un fluide newtonien.

partition de l'instabilité : le rapport d'aspect du jet $H/D > 3\pi$ et $Re < 0,52$, où $Re$ est le nombre de Reynolds basé sur le diamètre d'entrée $D$ et la vitesse d'entrée $V_0$. Ainsi, le but de cette étude est de vérifier la capacité de la méthode entièrement couplée à traiter ce type de problème et de fournir une bonne prédiction de leur comportement physique associé. En ce qui concerne les simulations numériques, les matrices résultant de la discrétisation seront résolues avec le même solveur, le BiCGStab(2) avec un résidu de $\epsilon = 10^{-5}$.

### 3.3.2 Résultats

Les figures 15-16 présentent l'évolution du liquide au cours du temps pour deux nombres de Reynolds : $Re = 0.1728$ et $Re = 0.6$, en utilisant un maillage de $100^3$ rafiné au niveau de l'injecteur. De ces figures, nous pouvons constaté que nos simulations numériques avec le solveur tout-couplé procurent de très bons résultats quant à la physique du jet, qui sont en bon accord avec le critère théorique donné par Tomé et al [22]. En effet, le premier cas $Re = 0.1728$ correspond à une instabilté du jet, alors que dans le deuxième, $Re = 0.6$ l'oscillation du jet n'apparaît pas.

## 4 Conclusion

Ce travail nous a permis de mettre en place un solveur tout-couplé parallèle 3D pour la simulation des écoulements diphasiques incompressibles. Nous avons également développé une méthode permettant de suivre la dynamique de l'interface au cours du temps, tout en conservant le volume à la précision machine. À travers les deux cas test que nous avons passé, notre solveur tout-couplé est très robuste lorsqu'il fait face à des problémes raides de fort rapports de masse volumique et de viscosité (comme
<figure><img src="image_23.png" /><figcaption>FIGURE 15 – Maillage 100³, $Re = 0.6$ : Le jet 3D visqueux au cours du temps. L'interface des deux fluides est représentée en rose.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 16 – Maillage 100³, $Re = 0.1728$ : Le jet 3D visqueux au cours du temps. L'interface des deux fluides est représentée en vert.</figcaption></figure>

dans le cas du jet visqueux). Ainsi, sur un problème de l'inversion de phase 3D, le solveur tout-couplé s'est montré très efficace pour simuler des problèmes diphasique sur un millards de cellules. Les futures étapes à envisager sont l'extention du solveur tout-couplé à la simulation des écoulements diphasiques compressibles avec ou sans changement de phase.

# 5 Remerciements

Nous remercions le CINES (Centre Informatique National de l'Enseignement Supérieur) et le TGCC accordé par GENCI sous les numéros de projet A0092B06115 pour l'accès aux moyens de calcul. Nous remercions également le personnel technique et administratif de ces centres et agences de supercalculateurs pour leur aide.

# Références

[1] N. Bootland, A. Bentley, C. Kees, and A. Wathen. Preconditioners for Two-Phase Incompressible Navier-Stokes Flow. *SIAM Journal on Scientific Computing*, 41(4):B843-B869, jan 2019.

[2] J.U Brackbill, D.B Kothe, and C. Zemach. A continuum method for modeling surface tension. *Journal of Computational Physics*, 100(2):335-354, 1992.
[3] J-P. Caltagirone and J. Breil. Sur une méthode de projection vectorielle pour la résolution des équations de navier-stokes. *Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy*, 327(11):1179-1184, 1999.

[4] J-P. Caltagirone and S. Vincent. A kinematics scalar projection method (ksp) for incompressible flows with variable density. *Open Journal of Fluid Dynamics*, 05, 03 2015.

[5] J-P. Caltagirone and S. Vincent. Sur une méthode de pénalisation tensorielle pour la résolution des équations de navier-stokes. *Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics*, 329(8):607-613, 2001.

[6] A-J. Chorin. Numerical solution of the navier-stokes equations. *Mathematics of Computation*, 22(104):745-762, 1968.

[7] D. Enright, F. Losasso, and R. Fedkiw. A fast and accurate semi-lagrangian particle level set method. *Computers Structures*, 83(6):479-490, 2005. Frontier of Multi-Phase Flow Analysis and Fluid-Structure.

[8] J.-L. Estivaleses, W. Aniszewski, F. Auguste, Y. Ling, L. Osmar, J.-P. Caltagirone, L. Chirco, A. Pedrono, S. Popinet, A. Berlemont, J. Magnaudet, T. Ménard, S. Vincent, and S. Zaleski. A phase inversion benchmark for multiscale multiphase flows. *Journal of Computational Physics*, 450:110810, 2022.

[9] M. Fortin and R. Glowinski. Augmented lagrangian methods : Applications to the numerical solution of boundary-value problems. ZAMM - *Journal of Applied Mathematics and Mechanics* / *Zeitschrift für Angewandte Mathematik und Mechanik*, 65(12):622-622, 1985.

[10] K. Goda. A multistep physics technique with implicit difference schemes for calculating two- or three-dimensional cavity flows. *Journal of computational physics* 30, 1979.

[11] J.L. Guermond and Jie Shen P. Minev. 2006-guermond-p6011-an overview of projection methods for incompressible flows. *Comput. Methods Appl. Mech. Engrg.*, 195:6011-6045, 2006.

[12] G.Weymouth and D.Yue. Conservative volume-of-fluid method for free-surface simulations on cartesian-grids. *J. Comput. Physics*, 229:2853-2865, 04 2010.

[13] D-C. Sorensen J-J. Dongarra, L-S. Duff and Henk A. Vander Vorst. *Numerical Linear Algebra for High Performance Computers*. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1998.

[14] I. Kataoka. Local instant formulation of two-phase flow. *International Journal of Multiphase Flow*, 12(5):745-758, 1986.

[15] R. J. LeVeque. High-resolution conservative algorithms for advection in incompressible flow. SIAM J. Numer. Anal., 33:627-665, 1996.

[16] Neelesh A. Patankar N. Nangia, Boyce E. Griffith and A-P-S Bhalla. A robust incompressible navier-stokes solver for high density ratio multiphase flows. *Journal of Computational Physics*, 390:548-594, 2019.

[17] M. El Ouafa, S. Vincent, and V. Le Chenadec. Monolithic solvers for incompressible two-phase flows at large density and viscosity ratios. *Fluids*, 6(1), 2021.

[18] M. El Ouafa, S. Vincent, and V. Le Chenadec. Navier-stokes solvers for incompressible single- and two-phase flows. *Communications in Computational Physics*, 29(4):1213-1245, 2021.

[19] C-H. Bruneau P. Angot and P. Fabrie. A penalization method to take into account obstacles in viscous flows. *Numerische Mathematik*, 81:497-520, 02 1999.
[20] M. Saeedipour, S. Vincent, and J.-L. Estivales. Toward a fully resolved volume of fluid simulation of the phase inversion problem. *Acta Mechanica*, 232(7):2695-2714, July 2021.

[21] R. Témam. Sur l'approximation de la solution des équations de navier-stokes par la méthode des pas fractionnaires (i). *Archive for Rational Mechanics and Analysis*, 32(2):135-153, Jan 1969.

[22] M-F. Tomé, B. Duffy, and S. McKee. A numerical technique for solving unsteady non-newtonian free surface flows. *Journal of Non-Newtonian Fluid Mechanics*, 62(1):9-34, 1996.

[23] S. Vincent, A. Sarthou, J.P. Caltagirone, F. Sonilhac, P. Février, C. Mignot, and G. Pianet. Augmented lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. application to hydroplaning flows interacting with real tire tread patterns. *Journal of Computational Physics*, 230(4):956-983, 2011.