# Caractérisation numérique des courbures et des normales aux interfaces d'écoulements diphasiques

Désir-André Koffi Bi

► To cite this version:

Désir-André Koffi Bi. Caractérisation numérique des courbures et des normales aux interfaces d'écoulements diphasiques. Sciences de l'environnement. Université Gustave Eiffel, 2021. Français.

NNT: tel-04086090

HAL Id: tel-04086090

https://hal.science/tel-04086090

Submitted on 1 May 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
École doctorale n° 531 : Sciences, Ingénierie et Environnement

THÈSE

pour obtenir le grade de

Docteur de l'Université Gustave Eiffel

Spécialité “Mécanique”

présentée et soutenue publiquement par

Désir-André Koffi Bi

le 13 Décembre 2021

Caractérisation numérique des courbures et des normales aux interfaces d'écoulements diphasiques

Devant le jury composé de:

<table><tr><td>J.-L. Estivalezes</td><td>Pr associé</td><td>IMFT - INP Toulouse</td><td>Rapporteur</td></tr><tr><td>S. Zaleski</td><td>Pr</td><td>∂'Alembert - Sorbonne Université</td><td>Rapporteur</td></tr><tr><td>M. Massot</td><td>Pr</td><td>École Polytechnique</td><td>Examineur</td></tr><tr><td>T. Sayadi</td><td>CR CNRS</td><td>∂'Alembert - Sorbonne Université</td><td>Examinatrice</td></tr><tr><td>S. Vincent</td><td>Pr</td><td>MSME - UGE</td><td>Directeur</td></tr><tr><td>E. Chénier</td><td>MCF-HDR</td><td>MSME - UGE</td><td>Co-directeur</td></tr></table>
# Résumé

## Caractérisation numérique des courbures et des normales aux interfaces d'écoulements diphasiques.

Les écoulements diphasiques sont rencontrés dans des domaines variés tels que le transport, la santé, l'industrie. Ils font l'objet de nombreux travaux, car ils possèdent une richesse et une diversité tant au niveau des phénomènes physiques s'y déroulant, des paramètres physicochimiques impliqués et des échelles rencontrées.

Quelle que soit la stratégie de simulation des écoulements diphasiques, un élément clé dans leur modélisation est la prise en compte des forces de tensions de surface. Elles font intervenir les propriétés géométriques de l'interface, à savoir la courbure et la normale. La bonne évaluation de ces éléments est capitale, car hormis la discrétisation des termes de forces de tensions de surface, une mauvaise estimation de la courbure et de la normale induit la création de courants parasites à l'interface. Par ailleurs, la courbure et la normale sont utilisées dans d'autres applications telles que l'évaporation, le décompte de gouttes dans un spray ou l'établissement de PDF (Probability Density Function) de tailles de gouttes. Ces données fournissent des informations globales (aire interfaciale, distribution de taille de gouttes) mais ne renseignent pas sur les différentes formes de gouttes/bulles potentiellement rencontrées dans un écoulement diphasique comme le montre le diagramme de Clift. C'est dans ce cadre que s'inscrit ce travail qui se veut être une étape préliminaire dans une caractérisation numérique et automatique des courbures et des normales aux interfaces d'écoulements diphasiques. La finalité est de pouvoir constituer une base de données d'informations permettant de caractériser des interfaces et d'opérer leur identification automatique.

Un des premiers volets consiste en la mise en place d'une méthode pour décrire l'interface 3D et suivre son évolution topologique au cours du temps. La technique du front-tracking (FT) a été choisie, car elle offre la connaissance explicite de la position de l'interface, premier jalon pour en déduire sa courbure et sa normale. Tout un ensemble de procédures 3D de gestion du maillage et de conservation/contrôle de volume ont été entièrement développées pendant ce travail et implémentées dans un code maison FUGU. Des tests de convergence ont également été menés.

Par la suite, tout un ensemble de méthodes existantes permettant d'accéder à la courbure et la normale grâce au maillage de l'interface issu du FT a été considéré avec des intercomparaisons. Deux cas de figure sont présentés pour rendre compte de la précision et la robustesse des méthodes quant à l'évaluation de la courbure et la normale : un test statique regroupant une grande variété de courbures et un test dynamique original issu de ce travail.

Enfin, une quantité statistique permettant de caractériser les interfaces est introduite et analysée. Elle utilise les courbures évaluées par les méthodes les plus robustes de la partie précédente. Elle est appliquée à différentes interfaces relatives au diagramme de Clift et les résultats sont interprétés.

**Mots-clés** — suivi d'interfaces 3D; diphasique; tension de surface; courbure; caractérisation d'interface
# Abstract

## Numerical characterisation of curvatures and normal vectors at multiphase flows interfaces.

Multiphase flows are encountered in various fields such as transportation, healthcare and industry. They are the subject of numerous studies as they possess a richness and diversity of physical phenomena, involve many physicochemical parameters and scales.

Whatever the strategy employed to simulate multiphase flows, a key element in their modelling is taking into account surface tensions. They depend on the geometric properties of the interface, namely the mean curvature and the normal vector. The correct evaluation of these elements is crucial because apart from the discretisation of the surface tensions terms, a wrong estimation of the curvature and the normal vector induces the creation of parasitic currents at the interface. In addition, the curvature and the normal vector are used in other applications such as evaporation, counting drops in a spray or establishing the PDF (Probability Density Function) of drop sizes. These data provide global information (interfacial area, drop size distribution) but do not enlighten on the different droplet/bubble shapes potentially encountered in a multiphase flow as shown in Clift diagram. It is within this framework that this work is intended to be a preliminary step in a numerical and automatic characterisation of curvatures and normal vectors at multiphase flows interfaces. The aim is to be able to build a database allowing the characterisation of interfaces and their automatic identification.

One of the first components consists of setting up a method for describing the 3D interface and follow its topological evolution over time. The front-tracking (FT) technique was chosen as it offers explicit knowledge of the position of the interface which is the first step in deducing its curvature and normal vector. A whole set of 3D procedures for mesh management and volume control/conservation were fully developed during this work and implemented in the FUGU in-house code. Convergence tests were also carried out.

Subsequently, a set of existing methods for accessing the curvature and normal vector through the interface mesh from the FT was considered with intercomparison. Two cases are presented to account for the accuracy and robustness of the methods in the evaluation of the curvature and the normal vector : a static test with a large variety of curvatures and an original dynamic test developed in this work.

Finally, a statistical quantity to characterise the interfaces is introduced and analysed. It uses the curvatures evaluated by the most robust methods of the previous part. It is applied to different interfaces relative to Clift diagram and the results are interpreted.

**Keywords** — front-tracking 3D; multiphase; surface tension; curvature; interface characterisation
# Table des matières

## 1 Introduction 19

1.1 Écoulements diphasiques : diversité et physique multi-échelle 19

1.2 Modélisation et simulation numérique d'écoulements diphasiques 22

1.2.1 Modèle 2-fluide 22

1.2.2 Modèle 1-fluide 24

1.2.3 La méthode Volume of Fluid (VOF) 26

1.2.4 La méthode Moment-Of-Fluid (MOF) 27

1.2.5 La méthode Level-Set (LS) 27

1.2.6 La méthode Front-Tracking (FT) 28

1.2.7 Les méthodes mixtes 29

1.3 De l'utilité d'une bonne évaluation des données géométriques de l'interface 30

1.3.1 Formulations des forces de tensions de surface 30

1.3.2 Estimation de la courbure et la normale 32

1.3.3 Exemples de l'utilité des approches de traitement des forces de tension de surface et de l'évaluation de la courbure et la normale 34

1.4 Plan de la thèse 36

## 2 Gestion du maillage de l'interface 45

2.1 Notation 45

2.2 Transport des marqueurs de l'interface 48

2.3 Enrichissement 49

2.3.1 Un nouvel élément : $(n_e^{(n)} - 3) = 1$ 52

2.3.2 Deux nouveaux éléments : $(n_e^{(n)} - 3) = 2$ 53

2.3.3 Trois nouveaux éléments : $(n_e^{(n)} - 3) = 3$ 54

2.4 Désenrichissement 55

2.5 Amélioration de la qualité du maillage 57

2.5.1 Échange d'arête 57

2.5.2 Lissage 60

2.5.3 Suppression de mauvais triangles 61

2.6 Conservation du volume 61

2.7 Validation - Interface sphérique dans un champ tourbillonnaire 63

## 3 Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces 71

3.1 Introduction 72

3.2 Num. methods for eval. $\kappa$, $n$ and $f$ 74

3.2.1 Main mathematical notations 74

3.2.2 Discretization of the Laplace-Beltrami Operator (LBO-methods) 75

3.2.3 Integral Formulation of the surface tension (IF-methods) 76

3.2.4 Surface Reconstruction (SR-methods) 78

3.2.5 Discrete surfaces and solution accuracy 79
<table>
  <tbody>
    <tr>
      <td>3.3</td>
      <td>Analysis of a stationary surface</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.3.1</td>
      <td>Results</td>
      <td>81</td>
    </tr>
    <tr>
      <td>3.3.2</td>
      <td>Comparisons</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.3.3</td>
      <td>Summary</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Analysis of advected surfaces</td>
      <td>88</td>
    </tr>
    <tr>
      <td>3.4.1</td>
      <td>Problem presentation</td>
      <td>88</td>
    </tr>
    <tr>
      <td>3.4.2</td>
      <td>Ligament formation - Case Q &lt; 0</td>
      <td>90</td>
    </tr>
    <tr>
      <td>3.4.3</td>
      <td>Bubble shape - Case Q &gt; 0</td>
      <td>93</td>
    </tr>
    <tr>
      <td>3.4.4</td>
      <td>Summary</td>
      <td>96</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Conclusion</td>
      <td>96</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Identification d'objets grâce à leurs propriétés surfaciques</td>
      <td>103</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Précision des méth. de calc. de κ sur un nuage de gouttes de forme sphé.</td>
      <td>103</td>
    </tr>
    <tr>
      <td>4.1.1</td>
      <td>Mise en œuvre de l'étude</td>
      <td>104</td>
    </tr>
    <tr>
      <td>4.1.2</td>
      <td>Grandeurs enregistrées</td>
      <td>106</td>
    </tr>
    <tr>
      <td>4.1.3</td>
      <td>Résultats et analyse</td>
      <td>106</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Echelles intermédiaires - Caractérisation des ligaments</td>
      <td>112</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Distribution surfacique de courbure</td>
      <td>112</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Méthodologie de l'étude numérique</td>
      <td>113</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Résultats et analyses</td>
      <td>119</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Ellipsoïdes déformées par des fonctions sinuées et variqueuses</td>
      <td>126</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>Déformation sinueuse</td>
      <td>130</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>Déformation variqueuse</td>
      <td>134</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Conclusion générale et perspectives</td>
      <td>145</td>
    </tr>
    <tr>
      <td>A</td>
      <td>La méthode de Cardan</td>
      <td>149</td>
    </tr>
    <tr>
      <td>B</td>
      <td>Calcul de la courbure par la méthode LBO - Cas d'un sommet à 3 triangles</td>
      <td>151</td>
    </tr>
    <tr>
      <td>C</td>
      <td>Complément sur les comparaisons de méthodes</td>
      <td>157</td>
    </tr>
    <tr>
      <td>C.1</td>
      <td>Errors for 20%-disturbed meshes</td>
      <td>157</td>
    </tr>
    <tr>
      <td>D</td>
      <td>Complément sur les méthodes de Reconstruction de Surface</td>
      <td>161</td>
    </tr>
    <tr>
      <td>D.1</td>
      <td>Peculiar stencil for the calculation of the SR-N<sub>1</sub>-n<sub>m</sub>-0 method</td>
      <td>161</td>
    </tr>
    <tr>
      <td>E</td>
      <td>Effet de la pondération dans les méthodes SR</td>
      <td>165</td>
    </tr>
    <tr>
      <td>E.1</td>
      <td>Analyse de l'infl. de la pond. et des perturb. sur la distrib. surf.</td>
      <td>165</td>
    </tr>
    <tr>
      <td>F</td>
      <td>Calcul exact de la distribution surfacique df<sub>S</sub> pour des ellipsoïdes</td>
      <td>173</td>
    </tr>
  </tbody>
</table>
Table des figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Exemples d'écoulements diphasiques, mettant en jeu différentes échelles.</td>
      <td>20</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Coalescence de deux bulles [18].</td>
      <td>23</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Différents méthodes de reconstruction de l'interface à partir de la fonction couleur C<sub>k</sub> issue de modélisation VOF.</td>
      <td>26</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Isosurfaces de la fonction Level-set Φ. L'interface est située à Φ = 0. Le signe de Φ permet de connaître la phase où on se situe.</td>
      <td>28</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Bulle/goutte modélisée par la méthode Front-Tracking (FT). Des marqueurs re-couvrent la surface et forment des triangles dont la taille moyenne des arêtes est d. La grille eulérienne a une résolution d'espace h.</td>
      <td>29</td>
    </tr>
    <tr>
      <td>1.6</td>
      <td>Évolution du maximum d'intensité Ca<sub>max</sub> des courants parasites dans le domaine de simulation au cours du temps (Fig.3 de [1]).</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Le volume Ω est délimité par une surface ∂Ω. La normale extérieure n&#773;<sub>ext</sub> pointe vers le milieu extérieur à Ω.</td>
      <td>45</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Numérotation des marqueurs de l'élément e selon les valeurs de n<sub>e</sub>.</td>
      <td>47</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>L'élément triangulaire initial e, constitué des marqueurs m<sub>[e,i]</sub>, i ∈ [1, 3], a ses arêtes 1 et 2 trop grandes, ce qui entraîne l'ajout des marqueurs supplémentaires m<sub>[e,4]</sub> et m<sub>[e,5]</sub> pour les scinder. L'orientation ou ordre de la numérotation est donnée par la normale sortante n&#773;<sub>ext</sub>.</td>
      <td>48</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>La surface de la sphère est maillée par des triangles. Les sommets sont transportés et suivent les changements topologiques de l'interface.</td>
      <td>49</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Détermination de la vitesse des marqueurs par la méthode PERM. Elle est obtenue à partir des vitesses adjacentes issues du maillage eulérien.</td>
      <td>49</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Sphère au maximum de déformation (t = T/2) dans un champ incompressible de vitesse (Eq. (2.10)). Les deux sphères ont le même maillage initial. A gauche (a), nous n'appliquons pas d'enrichissement. A droite (b), l'enrichissement est activé.</td>
      <td>50</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>L'élément e est constitué de 4 segments (a). Deux nouveaux éléments e' et e'' sont formés et l'élément e est supprimé. Les connectivités sont mises à jour (b).</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>L'élément e est formé de 5 segments (a). Trois nouveaux éléments e', e'' et e''' sont créés et l'élément e est supprimé. Le choix de construire les trois éléments et de les disposer dans cet ordre est fixé. Les voisins sont également mis à jour (b).</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td>L'élément e possède 6 segments (a). Quatre nouveaux éléments e', e'', e''' et e'''' sont créés tandis que l'élément e est supprimé. Les voisins sont également mis à jour (b).</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>L'arête 1 (rouge) de l'élément e remplit les conditions pour être supprimer (a). Les marqueurs des voisins ne changent pas. Seules les coordonnées de m<sub>[e',2]</sub> et m<sub>[e'',3]</sub> sont changées. Ils ont pour coordonnées le milieu de l'arête 1 de l'ancien élément e. Dans cette configuration, l'élément e' se retrouve renversé, situation qu'il faut éviter (b).</td>
      <td>56</td>
    </tr>
  </tbody>
</table>
2.11 L'arête 1 de l'élément e est courte et remplit les conditions pour être supprimer (a). Après la procédure, e et e' sont effacés. Les voisins sont également mis à jour (b). 57

2.12 L'arête 1 (rouge) de l'élément e est échangée pour rendre les éléments e et e' les plus équilatéraux possibles. 58

2.13 Variation de volume engendrée par l'échange d'arête entre les éléments e et e'. 58

2.14 L'arête 1 (rouge) de l'élément e est échangée. Le triangle e' voit sa normale inchangée après l'opération. Ce qui n'est pas le cas du triangle e, dont la normale pointe après vers l'intérieur. 59

2.15 Le lissage met $m_{[e,1]}$ au barycentre des marqueurs de son premier voisinage. 60

2.16 Un élément en aiguille (a) et un élément plat (b) à retirer du maillage lagrangien. 61

2.17 Le volume se calcule en faisant la somme des volumes des tétraèdres construits à partir des triangles de l'interface et d'un point de référence. 62

2.18 Évolution d'une interface sphérique dans un champ tourbillonnaire, avec $\Delta t = 10^{-3}$ et $h = 1/80$. De gauche à droite, l'interface est représentée aux instants 0, $T/10$, $T/5$ et $T/2$, le maximum de déformation. 64

2.19 Différents instants de l'interface à 0, $T/15$, $2T/15$, $3T/15$, $4T/15$, $5T/15$, $6T/15$, $7T/15$, $9T/15$, $11T/15$, $13T/15$ et $T$, tirés de [30]. 64

2.20 Écart entre un marqueur du maillage grossier et l'interface du maillage de référence. 65

2.21 Le marqueur m se trouve dans la cellule $(i_0, j_0, k_0)$. Les premiers éléments $e \in \mathcal{E}_{ref}$ rencontrés sont localisés dans la première couronne (zone blanche). Cependant, la distance $r_m$ établie avec ces éléments est plus grande que celle construite avec les éléments de la deuxième couronne (zone verte), $r'_m$. 65

2.22 Le point d'intersection entre la droite partant de m, de vecteur directeur $n_e$, et la plan ($P_e$) peut se trouver dans trois différentes zones. S'il est dans la zone 1 (intérieur de e), alors $I = H_1$. S'il est dans la zone 2, alors $I = \text{projeté orthogonal de } H_2$ sur l'arête. Dans la zone 3, I est un sommet de e. $H_3$ est proche de $m_{[e,2]}$ donc $I = m_{[e,2]}$. 66

2.23 Erreurs spatiales en norme $L_2$ (triangle) et $L_\infty$ (carré) en fonction de la taille h des cellules de discrétisation du domaine. 67

3.1 Notations and stencil for the discrete LBO methods [31, 76] 75

3.2 Construction of the mixed area $A_{mixed}$ of the discrete LBO, according to MEYER et al. [76] and ESSADKI [31] 76

3.3 Stencil for the original [123] method of the Integral Formulation (IF-Trygg. method) 77

3.4 Stencil for the 'triangle version' of the Integral Formulation (IF-triang. method). 77

3.5 Top (left) and side (right) views of the discrete mesh projected onto the surface (3.23). 81

3.6 LBO methods - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$ 81

3.7 LBO methods - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$ 82

3.8 Integral Formulation - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$ 83

3.9 Integral Formulation - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$ 83

3.10 Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$, $p = 0\%$ 84

3.11 Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\mathbf{n}$, $p = 0\%$ 84

3.12 Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$, $p = 20\%$ 84
<table>
  <tbody>
    <tr>
      <td>3.13</td>
      <td>Surface Reconstruction - Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on n, p = 20%</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.14</td>
      <td>Synthesis - Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 0%</td>
      <td>85</td>
    </tr>
    <tr>
      <td>3.15</td>
      <td>Synthesis - Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 20%</td>
      <td>86</td>
    </tr>
    <tr>
      <td>3.16</td>
      <td>Variation of the relative errors as a function of the maximum amplitude p of the perturbations</td>
      <td>86</td>
    </tr>
    <tr>
      <td>3.17</td>
      <td>Relative errors on κ in function of the time in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms, with (a) p = 0% and (b) p = 20%.</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.18</td>
      <td>Choice of the reference value to compute the relative errors.</td>
      <td>90</td>
    </tr>
    <tr>
      <td>3.19</td>
      <td>Bottom and side views of the discrete mesh transported in a spherical velocity field, Q &lt; 0</td>
      <td>91</td>
    </tr>
    <tr>
      <td>3.20</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms for κ and n, p = 0% with exact transport - Q &lt; 0</td>
      <td>91</td>
    </tr>
    <tr>
      <td>3.21</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms for κ and n, p = 0% with Lagrangian transport - Q &lt; 0</td>
      <td>92</td>
    </tr>
    <tr>
      <td>3.22</td>
      <td>Relative errors on κ (a) and n (b) according to CPU time in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms, p = 0% - Q &lt; 0 and numerical transport.</td>
      <td>93</td>
    </tr>
    <tr>
      <td>3.23</td>
      <td>Top and side views of the discrete mesh transported in a spherical velocity field, Q &gt; 0</td>
      <td>94</td>
    </tr>
    <tr>
      <td>3.24</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 0% with exact transport - Q &gt; 0</td>
      <td>94</td>
    </tr>
    <tr>
      <td>3.25</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 0% with Lagrangian transport - Q &gt; 0</td>
      <td>95</td>
    </tr>
    <tr>
      <td>3.26</td>
      <td>Relative errors on κ (a) and n (b) according to CPU time in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms, p = 0% - Q &gt; 0.</td>
      <td>96</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Gouttes sphériques avec différentes résolutions de courbure d<sub>r</sub> = 1/2r, r = 1, 10.</td>
      <td>104</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Perturbation maximale du maillage de l'interface de 5% et 40%.</td>
      <td>104</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Cône dans lequel est déplacé le marqueur m<sub>0</sub>.</td>
      <td>105</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Influence de la résolution de la courbure : PDF de l'écart relatif Err(&lt; κ &gt;), pour différentes résolutions d<sub>r</sub> et une perturbation p = 5%.</td>
      <td>109</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Influence de la résolution de la courbure : PDF de l'écart relatif Err(&lt; κ &gt;), pour différentes résolutions d<sub>r</sub> et une perturbation p = 20%.</td>
      <td>110</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Influence de la perturbation p : PDF de l'écart relatif Err(&lt; κ &gt;), pour différentes perturbations p et une résolution de courbure d<sub>r</sub> = 1/20.</td>
      <td>111</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>Formation de la famille d'ellipsoïdes à partir d'une sphère de résolution d<sub>r<sub>0</sub></sub> = 1/16.</td>
      <td>115</td>
    </tr>
    <tr>
      <td>4.8</td>
      <td>Évolution de df<sub>G</sub> en fonction du nombre N<sub>c</sub> de classes de courbure pour un ellipsoïde de paramètres d<sub>r<sub>0</sub></sub> = 1/256, p = 0%, α = 140%.</td>
      <td>116</td>
    </tr>
    <tr>
      <td>4.9</td>
      <td>Représentation des N<sub>c</sub> = 20 classes de courbures sur un ellipsoïde d'étirement α = 20% et de résolution de courbure. Dans la figure de droite, un zoom sur une classe de courbure (petit rectangle noir) permet de voir le maillage de l'interface.</td>
      <td>116</td>
    </tr>
    <tr>
      <td>4.10</td>
      <td>Sur-enrichissement de la zone de forte courbure pour un ellipsoïde de paramètre d<sub>r<sub>0</sub></sub> = 1/16, p = 0%, α = 20%</td>
      <td>117</td>
    </tr>
    <tr>
      <td>4.11</td>
      <td>Distribution surfacique df<sub>S</sub> (trait continu) et distribution de décompte df<sub>⊥</sub> (trait discontinu) pour deux stratégies de maillage de l'interface sur un ellipsoïde de paramètres d<sub>r<sub>0</sub></sub> = 1/16, p = 0%, α = 20%.</td>
      <td>117</td>
    </tr>
    <tr>
      <td>4.12</td>
      <td>Distribution surfacique df<sub>S</sub> pour des ellipsoïdes de paramètres d<sub>r<sub>0</sub></sub> = 1/256, p = 0%, α ∈ {20%, 80%, 140%}.</td>
      <td>119</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>B.2</td>
      <td>Aire utilisée pour construire A<sub>mixed</sub> dans le cas d'un triangle obtus.</td>
      <td>153</td>
    </tr>
    <tr>
      <td>C.1</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 20% with exact transport - Q &lt; 0</td>
      <td>157</td>
    </tr>
    <tr>
      <td>C.2</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 20% with Lagrangian transport - Q &lt; 0</td>
      <td>158</td>
    </tr>
    <tr>
      <td>C.3</td>
      <td>Relative errors on κ (a) and n (b) according to CPU time in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms, p = 20% - Q &lt; 0 and numerical transport.</td>
      <td>158</td>
    </tr>
    <tr>
      <td>C.4</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 20% with exact transport - Q &gt; 0</td>
      <td>158</td>
    </tr>
    <tr>
      <td>C.5</td>
      <td>Relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ and n, p = 20% with Lagrangian transport - Q &gt; 0</td>
      <td>159</td>
    </tr>
    <tr>
      <td>C.6</td>
      <td>Relative errors on κ (a) and n (b) according to CPU time in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms, p = 20% - Q &gt; 0 and numerical transport.</td>
      <td>159</td>
    </tr>
    <tr>
      <td>D.1</td>
      <td>Subdivision of the discrete mesh in three zones.</td>
      <td>161</td>
    </tr>
    <tr>
      <td>D.2</td>
      <td>Decomposition per zone of the relative errors in L<sub>2</sub> (continuous line) and L<sub>∞</sub> (dash line) norms on κ, p = 0% with exact transport - Q &gt; 0</td>
      <td>162</td>
    </tr>
    <tr>
      <td>D.3</td>
      <td>Stencil of the marker where the curvature and normal vector are not well evaluated.</td>
      <td>162</td>
    </tr>
    <tr>
      <td>D.4</td>
      <td>Relative error on κ in function of the angle of rotation θ of the markers M<sub>i</sub> ∈ (N<sub>2</sub>\N<sub>1</sub>), for the variant "Glimm - N<sub>1-1</sub>". The vertical dash lines are the principal directions.</td>
      <td>163</td>
    </tr>
    <tr>
      <td>E.1</td>
      <td>Evolution de la densité surfacique de courbure pour une famille d'ellipsoïdes de divers étirements.</td>
      <td>166</td>
    </tr>
    <tr>
      <td>E.2</td>
      <td>Erreurs relatives en norme L<sub>2</sub> (trait continu) et L<sub>∞</sub> (tiret) sur la densité surfacique de courbure d'une famille d'ellipsoïdes, avec la méthode LBO.</td>
      <td>169</td>
    </tr>
    <tr>
      <td>E.3</td>
      <td>Influence de la pondération - Erreurs relatives en norme L<sub>2</sub> (trait continu) et L<sub>∞</sub> (tiret) sur la densité surfacique de courbure d'une famille d'ellipsoïdes, à perturbation p = 0%, avec la méthode SR-N<sub>1</sub>-n<sub>m</sub>-0.</td>
      <td>170</td>
    </tr>
    <tr>
      <td>E.4</td>
      <td>Influence de la perturbation - Erreurs relatives en norme L<sub>2</sub> (trait continu) et L<sub>∞</sub> (tiret) sur la densité surfacique de courbure d'une famille d'ellipsoïdes, à paramètre N = 1, avec la méthode SR-N<sub>1</sub>-n<sub>m</sub>-0.</td>
      <td>171</td>
    </tr>
  </tbody>
</table>
# Nomenclature

$\Delta\kappa$ Pas de discrétisation de la plage de courbure sur un objet (section 4.2.1, page 112)

$\delta\Omega$ volume de contrôle (section 1.2.2, page 25)

$\partial\Omega$ surface délimitant le volume $\Omega$ (section 2.1, page 46)

$\hat{\kappa}$ Classe de courbures comprise entre $\kappa - \Delta\kappa/2$ et $\kappa + \Delta\kappa/2$ (section 4.2.1, page 112)

$\kappa$ Courbure moyenne (section 4.2.1, page 112)

$\Omega$ volume à simuler (section 2.1, page 46)

$\mathcal{D}$ ensemble caractérisant la discrétisation d'une interface 3D (section 2.1, page 46)

$\mathcal{E}$ ensemble des $N_{\mathcal{E}}$ triangles discrétisant la surface (section 2.1, page 46)

$\mathcal{M}$ ensemble des $N_{\mathcal{M}}$ marqueurs du maillage lagrangien (section 2.1, page 46)

$\mathcal{N}_{[k,m]}$ $k$-ème voisinage du marqueur $m$ (section 2.1, page 47)

($\mathcal{P}_e$) plan contenant l'élément $e$

$\mathcal{T}$ ensemble des $N_{\mathcal{T}}$ $n_e$-uplets définissant le sens de parcours des marqueurs (section 2.1, page 46)

$\mathcal{V}_e$ volume du tétraèdre construit à partir de l'élément $e$ et d'un point de référence

$\mathcal{V}$ volume de la bulle

$\overline{dl}_{e,i}$ vecteur de déplacement du marqueur $i$ de l'élément $e$ pour corriger le volume

$d_{e,i}$ longueur de l'arête $i$ de l'élément $e$ (section 2.1, page 47)

$DF_1(\kappa)$ Distribution de l'occurrence de $\kappa$ (section 4.2.1, page 112)

$DF_G(\eta)$ Distribution associée à la fonction $G$ et la quantité $\eta$ (section 4.2.1, page 112)

$DF_S(\kappa)$ Distribution surfacique de courbure (section 4.2.1, page 112)

$e_{[i]}$ voisin de l'élément $e$ par l'arête $i$ (section 2.1, page 47)

$G$ barycentre des marqueurs de l'interface

$h$ taille des cellules de la grille eulérienne

$i_0$ indice du marqueur précédant immédiatement le marqueur $i$ (section 2.1, page 46)

$i_2$ indice du marqueur suivant immédiatement le marqueur $i$ (section 2.1, page 46)

$m_{[e,i]}$ marqueur $i$ de l'élément $e$ (section 2.1, page 46)

$\bar{n}_{e,i}$ normale unitaire au marqueur $i$ de l'élément $e$

$\bar{n}_{ext}$ normale extérieure à la surface $\partial\Omega$ (section 2.1, page 46)

$\bar{n}_e$ normale unitaire de l'élément $e$
$n_e$ nombre de segments ou de marqueurs constituant l'élément $e$ (section 2.1, page 46)

$Q_e$ qualité de l'élément $e$

$r_{[e,i]}$ rugosité de l'arête $i$ de l'élément $e$ (section 2.4, page 55)

$\bar{v}_{e,i}$ vitesse du marqueur $i$ de l'élément $e$ (section 2.1, page 46)

$\bar{x}_{e,i}$ coordonnées du marqueur $i$ de l'élément $e$ (section 2.1, page 46)

$\bar{x}_G$ coordonnées du barycentre $G$
# Chapitre 1

# Introduction

## Sommaire

<table><tr><td>1.1</td><td>Écoulements diphasiques : diversité et physique multi-échelle</td><td>19</td></tr><tr><td>1.2</td><td>Modélisation et simulation numérique d'écoulements diphasiques</td><td>22</td></tr><tr><td>1.2.1</td><td>Modèle 2-fluide</td><td>22</td></tr><tr><td>1.2.2</td><td>Modèle 1-fluide</td><td>24</td></tr><tr><td>1.2.3</td><td>La méthode Volume of Fluid (VOF)</td><td>26</td></tr><tr><td>1.2.4</td><td>La méthode Moment-Of-Fluid (MOF)</td><td>27</td></tr><tr><td>1.2.5</td><td>La méthode Level-Set (LS)</td><td>27</td></tr><tr><td>1.2.6</td><td>La méthode Front-Tracking (FT)</td><td>28</td></tr><tr><td>1.2.7</td><td>Les méthodes mixtes</td><td>29</td></tr><tr><td>1.3</td><td>De l'utilité d'une bonne évaluation des données géométriques de l'interface</td><td>30</td></tr><tr><td>1.3.1</td><td>Formulations des forces de tensions de surface</td><td>30</td></tr><tr><td>1.3.2</td><td>Estimation de la courbure et la normale</td><td>32</td></tr><tr><td>1.3.3</td><td>Exemples de l'utilité des approches de traitement des forces de tension de surface et de l'évaluation de la courbure et la normale</td><td>34</td></tr><tr><td>1.4</td><td>Plan de la thèse</td><td>36</td></tr></table>

## 1.1 Écoulements diphasiques : diversité et physique multi-échelle

Les écoulements diphasiques à phases séparées se réfèrent aux mouvements de substances de même composition chimique évoluant simultanément dans des phases différentes (eau liquide et vapeur d'eau) ou de nature hétérogène mais se présentant dans des phases similaires (eau liquide et huile) ou distinctes (gasoil et air). Cette description nous fait prendre conscience de la diversité des écoulements diphasiques et de leur présence dans de nombreux domaines, tant au niveau de notre vie quotidienne, de la nature que du monde industriel.

Les échanges entre les océans et l'atmosphère sont un des premiers exemples d'écoulements diphasiques. En effet, notre planète est recouverte à 70% d'eau et les océans constituent la plus grande expérience à ciel ouvert d'écoulements diphasiques. Les phénomènes d'évaporation, le déferlement de vagues et la formation d'un nuage de gouttelettes ou embruns (Fig. 1.1a), le transfert de CO₂ entre les océans et l'atmosphère sont autant de paramètres qui influent sur le climat, la météorologie planétaire et nécessitent encore de nombreuses études [27, 66, 126, 137].

Le déferlement des vagues de l'exemple précédent nous permet d'introduire un autre aspect des écoulements diphasiques qui est l'atomisation de jets. Elle consiste en la fragmentation du jet d'un liquide en des structures de plus en plus petites, allant des ligaments aux gouttes. On assiste à une cascade d'échelles. Le jet initial continu représente les grandes échelles. Il se scinde ensuite en ligaments plus ou moins longs : ces derniers forment l'échelle intermédiaire. Enfin, la
<figure><img src="image_2.png" /><figcaption>(a) Déferlement de vagues sur des rochers avec formation d'un nuage de fines gouttelettes (embruns). [https://besthqwallpapers.com/fr/download/original/52962]</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c) Format d'un film liquide avec apparition de digitations après l'impact d'une goutte dans une flaque de lait. [https://publicdomainpictures.net/pictures/50000/velka/dopad-kapky-z-mleka-1.jpg]</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Fragmentation d'un jet. Des ligaments de tailles variés se forment et des gouttes satellitaires apparaissent [70].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(d) Différents modes de génération de gouttelettes. Les dimensions des canaux sont de l'ordre de $10\mu m – 100\mu m$ [4].</figcaption></figure>

<caption>FIGURE 1.1 – Exemples d'écoulements diphasiques, mettant en jeu différentes échelles.</caption>

scission se poursuit et s'achève avec la formation de petites gouttes sphériques qui doivent leur forme à la prédominance des forces de tension de surface (Fig. 1.1b). La collection de gouttes formées lors de la fragmentation du jet s'appelle un spray. Partagée en deux étapes (atomisations primaire et secondaire), elle est le lieu de plusieurs phénomènes physiques (turbulence, instabilités de Kelvin-Helmholtz, Rayleigh-Taylor, Rayleigh-Plateau, vitesse de Taylor-Culick pour le percement d'une nappe liquide) qui influent son évolution. L'ensemble de ces phénomènes est piloté par des nombres adimensionnels comme le nombre de Reynolds *Re*, le nombre de Weber *We*, le nombre d'Ohnesorge *Oh* et également les rapports de quantité de mouvement entre le gaz et le liquide injecté. De nombreux travaux se sont penchés sur les différentes étapes de l'atomisation et les phénomènes physiques impliqués [8, 10, 28, 47, 49, 61, 62, 70, 102]. Nous pouvons citer ainsi FUSTER et al. [35], qui en combinant des approches VOF (Volume Of Fluid) et AMR (adaptive mesh refinement) ont pu identifier deux régimes lors de la retrait d'un film liquide après atomisation qui dépendent du nombre d'Ohnesorge. Outre la compréhension physique de l'atomisation et des paramètres la gouvernant, une manière d'étudier l'atomisation de jet est d'établir la distribution ou PDF (Probability Density Function) des tailles de gouttes constituant le spray, pour ne citer que celle-là. BABINSKY et al. [7] présentent trois approches pour établir la PDF d'un spray : la méthode empirique, la méthode du maximum d'entropie et celle de la fonction de probabilité discrète. Chacune des méthodes est adaptée pour décrire
les différentes étapes de l'atomisation. Dans le livre de LEFEBVRE et al. [63], différentes lois mathématiques de probabilité sont présentées pour modéliser la PDF d'un spray. Ces PDF sont également utiles dans d'autres domaines comme la combustion dans les moteurs, le milieu de la santé, l'agriculture, les industries pharmaceutique et agroalimentaire.

Dans le milieu de la santé, plus particulièrement de la transmission des maladies, connaître les sites d'émissions de gouttelettes générées par l'activité humaine de même que leurs tailles permet de mieux comprendre la propagation de virus aéroportés et par conséquent de mieux les combattre. Dans cette catégorie, nous pouvons citer les travaux de MORAWSKA et al. [78] sur la distribution de tailles de gouttes émises lors d'activités expiratoires humaines.

Dans l'industrie automobile, enrichir les modèles d'atomisation et établir les PDF de tailles de gouttes jouent un rôle primordial dans l'amélioration de l'efficacité des moteurs. En effet, un des points clé est l'optimisation de la surface d'échange entre le carburant et l'air. Avoir la plus grande aire interfaciale favorise la combustion du carburant pour récupérer le maximum d'énergie. Ainsi, dans les travaux de SOM et al. [112], la prise en compte de phénomènes de cavitation et de turbulences dans le jet sortant de l'injecteur lors de l'atomisation primaire permet d'améliorer les résultats des simulations et d'être en concordance avec les données expérimentales. CHEN et al. [16] montrent que différents carburants conduisent à des PDF de taille de gouttes distinctes, mais surtout que l'augmentation de la pression d'injection réduit la taille des gouttes. LUO et al. [68] ont investigué par DNS (Direct Numerical Simulation) la combustion d'un spray de n-heptane et ont montré toute la complexité du phénomène.

Cette cascade d'échelles se retrouve également dans les problèmes d'impacts de gouttes (Fig. 1.1c) ou de génération de gouttelettes dans des microcanaux (Fig. 1.1d). Dans le premier, il est intéressant de déterminer les paramètres influençant l'étalement sur les surfaces d'impact [15, 79, 80, 127, 129]. Dans le second, les phénomènes physiques prédominants ne sont plus les mêmes et le contrôle précis du volume des gouttelettes ainsi engendrées permet dans l'industrie pharmaceutique par exemple de créer de mini-réacteurs pour des synthèses chimiques [4, 111].

Jusqu'ici, nous avons présenté des exemples où la cascade d'échelle s'effectue des grandes échelles aux plus petites. Il existe des situations où cette cascade se passe dans le sens inverse. C'est le cas de la coalescence de gouttes après recouvrement d'une surface par un spray. Ce problème est d'importance dans les domaines agroalimentaires et pharmaceutiques avec les processus dits de séchage de sprays. Ils servent par exemple à la production de poudres initialement mélangées dans un liquide qui est atomisé puis séché par évaporation. L'état des parois sur lesquels le spray peut coalescer avec formation de dépôts, les distributions des tailles de gouttes, les aires interfaciales ont un rôle primordial sur la qualité finale du produit [55, 60, 86]. Nous retrouvons également cette cascade inverse dans les écoulements à bulles dans les conduites (échangeurs, centrales nucléaires). Selon le flux de chaleur échangé entre les parois et le liquide et selon le régime d'écoulement, différents types d'écoulements liquide/vapeur se forment [12]. Un des phénomènes crucial à éviter dans les centrales nucléaires est la crise d'ébullition, première étudiée expérimentalement par NUKIYAMA [82]. Lorsque le flux de chaleur transmis de la surface vers le fluide dépasse une valeur critique appelée Critical Heat Flux (CHF), les bulles générées coalescent et forment un film séparant la surface chauffée du liquide. La chaleur est moins bien évacuée et entraîne la surchauffe de la surface, pouvant engendrer sa détérioration voire sa fusion. Des études sont menées afin de mieux comprendre cette situation et éviter des accidents nucléaires [46, 88].

L'ensemble de ces exemples nous montrent la diversité et la complexité des écoulements diphasiques. Outre les enjeux scientifiques tels que la détermination et la compréhension des phénomènes les gouvernant, il existe également une dimension industrielle à leur étude pour l'obtention de paramètres macroscopiques permettant l'optimisation de systèmes (rétablissement de PDF, obtention d'aires interfaciales, de vitesse, détection de la forme des objets). Cette richesse et cette complexité se retrouvent également dans les approches développées pour leur
simulation et modélisation.

## 1.2 Modélisation et simulation numérique d'écoulements diphasiques

La modélisation et la simulation des écoulements diphasiques demeurent encore un challenge du aux différentes échelles qui coexistent simultanément. Diverses approches existent selon le point de vue choisi. A l'échelle moléculaire, ces sont des modèles dits de dynamique moléculaire qui permettent de passer de l'aspect discontinu de l'écoulement à des grandeurs macroscopiques par des analyses statistiques. A l'échelle macroscopique, le milieu est supposé continu. Les équations locales sont définies sur de petits volumes appelés particule-fluide. Notre travail porte sur ce dernier volet où nous considérons deux fluides non-miscibles. Les équations régissant leur évolution dans le cadre d'interfaces résolues sont abordées dans un premier temps. Puis différentes approches pour modéliser et suivre l'interface sont présentées.

### 1.2.1 Modèle 2-fluide

Les équations de Navier-Stokes sont les relations qui permettent de suivre l'évolution de grandeurs physiques (masse volumique, quantité de mouvement, température, énergie) pour les écoulements fluides. Dans un écoulement diphasique où deux phases sont mises en jeu, une manière naturelle de décrire leur évolution est la résolution dans chaque phase des équations de Navier-Stokes. Nous ne considérons que l'évolution isotherme d'écoulements diphasiques. Ainsi, pour chaque phase $k$, nous résolvons les équations de conservation de la masse et de la quantité de mouvement, respectivement présentées dans les relations (1.1) et (1.2) ; $\rho_k$ correspond à la masse volumique, $\mu_k$ la viscosité dynamique, $\mathbf{u}_k$ la vitesse et $\tau_k$ le tenseur des contraintes associés à la phase $k$.

$$ \frac{\partial}{\partial t}(\rho_k) + \nabla \cdot (\rho_k \mathbf{u}_k) = 0 \qquad (1.1) $$

$$ \frac{\partial}{\partial t}(\rho_k \mathbf{u}_k) + \nabla \cdot (\rho_k \mathbf{u}_k \otimes \mathbf{u}_k) = \rho_k \mathbf{g} + \nabla \cdot \tau_k \qquad (1.2) $$

Le tenseur des contraintes est donné par la relation (1.3).

$$ \tau_k = -p_k \mathbb{I} + \tau'_k \qquad (1.3) $$

Le tenseur $\mathbb{I}$ est la matrice identité, $p_k$ la pression et $\tau'_k$ est le tenseur des contraintes visqueuses. Pour un fluide newtonien, il s'exprime suivant la relation (1.4).

$$ \tau'_k = 2\mu_k S_k - \frac{2}{3}\mu_k (\nabla \cdot \mathbf{u}_k) \mathbb{I} \qquad (1.4) $$

$S_k$ est le tenseur des taux de déformations et est donné par la relation (1.5).

$$ S_k = \frac{1}{2} (\nabla \mathbf{u}_k + (\nabla \mathbf{u}_k)^t) \qquad (1.5) $$

Dans le cas d'écoulements incompressibles ($\nabla \cdot \mathbf{u}_k = 0$), le tenseur de contraintes visqueuses $\tau'_k$ est simplifié pour aboutir à la relation (1.6).
$$ \tau'_{k} = \mu_{k} (\nabla \mathbf{u}_{k} + (\nabla \mathbf{u}_{k})^{t}) \qquad (1.6) $$

Au niveau de l'interface entre les phases, un raccord est effectué à travers les relations de saut développées par DELHAYE [21]. Les équations (1.7) et (1.8) représentent respectivement les relations de saut de vitesses (ou continuité de la vitesse dans le cas $\dot{m}_i = 0$) et de saut de contraintes.

$$ [\rho_k(\mathbf{u}_k - \mathbf{V}_i)] \cdot \mathbf{n} = \dot{m}_i \qquad (1.7) $$

$$ [-\dot{m}_k \mathbf{u}_k + \tau_k \cdot \mathbf{n}_k] = \sigma \kappa \mathbf{n} + \nabla \sigma \qquad (1.8) $$

L'opérateur [•] dans les relations (1.7) et (1.8) est une différence effectuée dans le sens de la normale $\mathbf{n}$, avec $\mathbf{n}_k$ la normale à l'interface pointant à l'extérieur de la phase $k$. Ainsi si $\mathbf{n} = \mathbf{n}_1$, la différence est effectuée de la phase 1 vers la phase 2 (et inversement) ; $\sigma$ correspond au coefficient de tension de surface et $\kappa$ au double de la courbure moyenne, $\mathbf{V}_i$ est la vitesse de l'interface, $\dot{m}_k$ le flux de masse de la phase $k$ à l'interface, $\dot{m}_i$ la génération de masse à l'interface (1.9).

$$ \sum_{k=1}^{2} \dot{m}_k = \dot{m}_i \quad \text{avec} \quad \dot{m}_k = \rho_k(\mathbf{u}_k - \mathbf{V}_i) \cdot \mathbf{n}_k \qquad (1.9) $$

L'interface étant sans épaisseur et sans masse, le terme $\dot{m}_i$ est nul en absence de changement de phase.

<figure><img src="image_6.png" /><figcaption>FIGURE 1.2 – Coalescence de deux bulles [18].</figcaption></figure>

Cette approche de l'étude des écoulements diphasiques est beaucoup employée avec la méthode Arbitrary Lagrangian-Eulerian (ALE) [18, 20, 36, 37, 38, 64]. Un exemple est donné sur la figure 1.2 qui illustre la coalescence de deux bulles. La méthode ALE consiste à modéliser le milieu fluide par un maillage non-structuré qui se déforme et s'adapte au cours de l'évolution de l'interface. Elle a l'avantage de capter et s'adapter à des géométries complexes. Toutefois, le revers de la médaille se situe dans son coût numérique élevé. En effet, à chaque itération temporelle de la simulation, le maillage est de nouveau généré. Réaliser cette opération en trois dimensions d'espace pour des configurations complexes de rupture ou de coalescence d'interface n'est en général pas possible.
1.2.2 Modèle 1-fluide

L'idée est d'avoir un seul jeu d'équations qui rend compte de l'évolution des deux fluides. Le formalisme du modèle 1-fluide a été introduit par Kataoka [54] dans les années 80. Le lecteur est également orienté vers [69] pour une démarche détaillée de la démonstration. Dans un premier temps, le maillage non-structuré du domaine de simulation est abandonné pour laisser place à un maillage structuré sur lequel sont résolues les équations de Navier-Stokes. Pour mettre en place cet unique système d'équations décrivant l'évolution des deux fluides, Kataoka introduit une fonction caractéristique $\phi_k(\bar{x}, t)$ associée à la phase $k$. Elle prend la valeur 1 quand le point $\bar{x}$ où elle est évaluée se situe dans la phase $k$. Dans le cas contraire, elle est nulle (1.10).

$$
\phi_k(\bar{x}, t) = \begin{cases} 1, & \bar{x} \in \text{phase } k \\ 0, & \text{sinon} \end{cases} \tag{1.10}
$$

Par ailleurs, la fonction $\phi$ possède d'autres propriétés indiquées dans les relations (1.11), (1.12) et (1.13).

$$
\forall \bar{x}, \forall t, \quad \sum_{k=1}^{2} \phi_k(\bar{x}, t) = 1 \tag{1.11}
$$

$$
\frac{\partial}{\partial t}(\phi_k) + \mathbf{V}_i \cdot \nabla \phi_k = 0 \qquad (1.12)
$$

$$
\nabla \phi_k = -\mathbf{n}_k \delta_i \tag{1.13}
$$

Le vecteur **n**_k est la normale à l'interface pointant à l'extérieur de la phase k et δ_i est la fonction
de Dirac localisée au niveau de l'interface. En effet, la fonction φ_k est discontinue à l'interface,
passant de la valeur 1 à 0. Sa dérivée est donc un Dirac centré à l'interface entre les deux phases.
Par ailleurs la variation est dans la direction de la normale à l'interface. La première étape pour
aboutir à une formulation 1-fluide dans le cadre d'écoulements incompressibles (∇ · **u**_k = 0)
est la multiplication par φ_k des équations de conservation (1.1) (1.2) et leur sommation sur
l'ensemble des phases. Par la suite, l'utilisation des relations de saut (1.7) (1.8) avec l'hypothèse
ṁ_i = 0 et des propriétés de φ_k (1.11) (1.12) (1.13) permet une simplification des expressions et
l'aboutissement à une formulation des équations de conservation incluant les relations de saut
par la présence des forces de tension de surface (1.14) à (1.17).

$$
\nabla \cdot \left( \sum_{k=1}^{2} \phi_k \mathbf{u}_k \right) + \left( \sum_{k=1}^{2} \mathbf{u}_k \cdot \mathbf{n}_k \right) \delta_i = 0 \quad (1.14)
$$

$$
\frac{\partial}{\partial t} \left( \sum_{k=1}^{2} \phi_k \rho_k \right) + \nabla \cdot \left( \sum_{k=1}^{2} \phi_k \rho_k \mathbf{u}_k \right) = 0 \quad (1.15)
$$

$$
\frac{\partial}{\partial t} \left( \sum_{k=1}^{2} \phi_k \mu_k \right) + \nabla \cdot \left( \sum_{k=1}^{2} \phi_k \mu_k \mathbf{u}_k \right) = - \left( \sum_{k=1}^{2} \phi_k \mu_k (\mathbf{u}_k - \mathbf{V}_i) \cdot \mathbf{n}_k \right) \delta_i \quad (1.16)
$$

$$
\frac{\partial}{\partial t} \left( \sum_{k=1}^{2} \phi_k \rho_k \mathbf{u}_k \right) + \nabla \cdot \left( \sum_{k=1}^{2} \phi_k \rho_k \mathbf{u}_k \otimes \mathbf{u}_k \right) = \sum_{k=1}^{2} \phi_k \rho_k \mathbf{g} + \nabla \cdot \left( \sum_{k=1}^{2} \phi_k \tau_k \right) + \sigma \kappa \mathbf{n} \delta_i \quad (1.17)
$$
La dernière étape consiste à effectuer un filtrage spatial car les équations (1.15) (1.17) comportent des discontinuités introduites par la fonction caractéristique $\phi_k$. Ce filtrage permet de travailler avec des grandeurs continues. Ainsi, un filtre spatial $G$ est introduit et remplat la condition (1.18).

$$ \int_{\delta\Omega} G(\bar{x} - \bar{x}') d\bar{x}' = 1 \qquad (1.18) $$

Calculé sur un volume de contrôle $\delta\Omega$ (généralement la maille de la grille de discrétisation spatiale), il est appliqué à la fonction $\phi_k$ pour donner sa contrepartie filtrée $C_k$ (1.19), et également à toute autre grandeur $\Psi_k$ du problème (1.20), avec $\Psi_k \in \{\rho_k, \mu_k, \mathbf{u}_k, \tau_k\}$. La fonction filtrée $C_k$ correspond à la fraction volumique de la phase $k$ dans le volume de contrôle $\delta\Omega$.

$$ C_k(\bar{x}, t) = \int_{\delta\Omega} \phi_k(\bar{x}', t) G(\bar{x} - \bar{x}') d\bar{x}' \qquad (1.19) $$

$$ C_k < \Psi_k > (\bar{x}, t) = \int_{\delta\Omega} \Psi_k(\bar{x}', t) \phi_k(\bar{x}', t) G(\bar{x} - \bar{x}') d\bar{x}' \qquad (1.20) $$

On définit également la moyenne $\Psi \in \{\rho, \mu, \mathbf{u}, \tau\}$ des grandeurs filtrées (1.21).

$$ \Psi = \sum_{k=1}^{2} C_k < \Psi_k > \qquad (1.21) $$

A ce stade, des hypothèses doivent être émises, dont les deux premières sont fondamentales.

- L'absence de changement de phase.

- L'absence de glissement entre les phases. Cela se traduit par le fait que $< \mathbf{u}_1 > = < \mathbf{u}_2 >$.

- Négligence des corrélations des fluctuations de vitesse qui peuvent exister à une échelle plus petite que celle du filtrage spatial : $< \mathbf{u}_1 \otimes \mathbf{u}_1 > = < \mathbf{u}_2 \otimes \mathbf{u}_2 >$.

- En moyenne (au sens du filtrage), les gradients sont les mêmes en un point donné, quelle que soit la phase : $< S_1 > = < S_2 >$.

L'emploi de ces hypothèses avec l'application du filtre aux équations (1.16) et (1.17) nous donne finalement la formulation 1-fluide pour des écoulements diphasiques (1.22) (1.23).

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1.22) $$

$$ \frac{\partial}{\partial t}(\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = \rho \mathbf{g} + \nabla \cdot \tau + \sigma \kappa \mathbf{n} \delta_i \qquad (1.23) $$

A cela, il faut ajouter l'équation d'advection de la fonction $C_k$ (1.24).

$$ \frac{\partial}{\partial t}(C_k) + \mathbf{u} \cdot \nabla C_k = 0 \qquad (1.24) $$
Le terme $\sigma\kappa n\delta_i$ représente la contribution des tensions de surface dont les paramètres $\kappa$ et $n$ dépendent des propriétés géométriques de l'interface. Dans la formulation 1-fluide, ils sont donnés par la relation (1.25).

$$ \kappa = -\nabla \cdot n \quad \text{avec} \quad n = \nabla C_k / \|\nabla C_k\| \quad \text{et} \quad n\delta_i = \nabla C_k \qquad (1.25) $$

Quel que soit le modèle utilisé (1-fluide ou 2-fluides) pour la simulation des écoulements diphasiques, l'évaluation de la courbure $\kappa$ et de la normale à l'interface $n$ est importante pour la prise en compte des forces de tensions de surface. Elles dépendent de l'interface dont la modélisation regroupe diverses approches. On distingue ainsi les approches dites de capture d'interface qui accèdent à la position de l'interface de manière implicite. Elles font face aux approches dites de suivi d'interface qui ont la connaissance explicite de celle-ci. La section suivante est consacrée à une présentation non exhaustive de ces méthodes.

### 1.2.3 La méthode Volume of Fluid (VOF)

Elle fait partie des méthodes dites de capture d'interface et fut introduite par HIRT et al. [50]. La position de l'interface se déduit implicitement à partir d'une fonction couleur $C$. Dans une maille de la grille du domaine discrétisé, elle représente la fraction volumique de la phase à laquelle elle est associée. Ainsi, la couleur $C_k$ associée à la phase $k$ faut 1 à l'intérieur de la phase $k$, 0 à l'extérieur et une valeur intermédiaire au niveau de l'interface. L'évolution de l'interface est captée par l'équation d'advection de la fonction couleur $C_k$ (1.26). On retrouve une forme conservative lorsque l'écoulement est incompressible ($\nabla \cdot u = 0$).

$$ \frac{\partial}{\partial t}(C_k) + \nabla \cdot (uC_k) - C_k \nabla \cdot u = 0 \qquad (1.26) $$

La résolution numérique de cette équation pose toutefois des problèmes de diffusion de l'interface et de lissage des propriétés physiques à son voisinage à cause de la discontinuité de la fonction $C_k$. La résolution de ce problème s'effectue par une double approche.

<figure><img src="image_7.png" /><figcaption>FIGURE 1.3 – Différents méthodes de reconstruction de l'interface à partir de la fonction couleur $C_k$ issue de modélisation VOF.</figcaption></figure>

La première consiste à chaque fois à reconstruire l'interface avant son advection. Pour ce faire, différentes techniques ont été développées. Des exemples sont présentés sur la figure 1.3. La méthode SLIC (Simple Line Interface Calculation) introduite par [81] représente l'interface par des segments parallèles à la grille eulérienne. La position et l'orientation des segments dépend des couleurs $C_k$ dans un bloc 3 × 1 centré autour de la cellule d'intérêt (Fig. 1.3b). Une amélioration à la représentation de l'interface a été apportée par YOUNGS [135] avec la méthode PLIC
(Piecewise Linear Interface Calculation). Dans chaque cellule où se situe l'interface, celle-ci est toujours approximée par une droite dont l'orientation est ici déterminée par une approximation de la normale à l'interface. Des discontinuités subsistent cependant entre les interfaces reconstruites de cellules adjacentes (Fig. 1.3c). Dans le cadre de méthodes d'ordre élevé, nous pouvons citer les méthodes LVIRA (Least-squares Volume-of-fluid Interface Reconstruction Algorithm) de PUCKETT [96] et son amélioration ELVIRA (E pour Efficient) de PILLIOD [89] qui sont du second ordre. Dans la première, la pente du segment d'interface (ou la normale) est évaluée en minimisant une erreur ($L_1$, $L_2$ ou $L_\infty$ au choix de l'utilisateur) entre la fonction $C_k$ et une approximation $\tilde{C}_k$ de cette dernière. L'amélioration apportée par la méthode ELVIRA consiste à augmenter le choix possible de pentes (6 en 2D) en considérant en plus les différences finies décentrées amont et aval pour le calcul de gradient de $C_k$. La minimisation est ensuite opérée sur ces différents candidats et la pente retenue est celle qui présente le minimum d'erreurs. Le lecteur est orienté vers [90] pour une description plus détaillée des méthodes présentées auparavant. SCARDOVELLI et al. [103] pour sa part propose, en plus d'une approximation linéaire, une approximation quadratique de l'interface, permettant de la représenter par une succession de deux segments. Plus récemment, [5] propose une alternative dans la cadre de la méthode PLIC pour reconstruire l'interface, pouvant s'avérer jusqu'à 100 fois plus rapide. Nommée NPLIC (Neural Piecewise Linear Interface Calculation), elle utilise un réseau de neurones pour trouver l'orientation et la position de l'interface.

La seconde approche se focalise sur les schémas de discrétisation de l'équation d'advection (1.26). SCARDOVELLI et al. [103] propose ainsi un schéma mixte d'advection EI-LE (Eulerian Implicit - Lagrangian Explicite) qui combiné à leur méthode de reconstruction d'interface donnent de meilleurs résultats que des approches VOF-PLIC. Une version 3D est présentée dans [6]. D'autres approches par contre ne passent pas par une reconstruction de l'interface. Dans la méthode THINC (Tangent of Hyperbola for INterface Capturing) de XIAO et al. [130], le saut introduit par la fonction caractéristique $C_k$ est modélisé par une tangent hyperbolique, permettant le calcul de manière algébrique des flux dans la relation (1.26). Une amélioration en est proposée par [131].

Un des grands avantages de la méthode VOF est la conservation de volume lors des simulations. Elle permet également de gérer automatiquement les phénomènes de coalescence et de rupture. Toutefois, les évaluations des propriétés de l'interface, à savoir la courbure $\kappa$ et la normale $n$ demeurent de faible qualité ou complexes (exemple avec la méthode Smooth VOF de [48]).

### 1.2.4 La méthode Moment-Of-Fluid (MOF)

Introduite en 2005 par DYADECHKO et al. [29], la méthode MOF (Moment-Of-Fluid) se veut plus précise et d'ordre supérieure (ordre 2 de convergence) à la méthode VOF [3, 58]. En plus de la fonction couleur $C_k$ qui est transportée, son centroïde est également transporté. Par ailleurs, contrairement à la méthode VOF qui nécessite l'information des cellules adjacentes pour la reconstruction de l'interface, elle utilise uniquement l'information de la cellule eulérienne concernée. Un des avantages de cette méthode est également d'arriver à représenter des éléments sous-maille comme les ligaments [52]. Bien que l'approche MOF puisse traiter d'éléments multi-matériaux (supérieur à 2), elle peut devenir coûteuse numériquement. En ce sens récemment, il est proposé une alternative qui permet d'améliorer les performances numériques de la méthode [19].

### 1.2.5 La méthode Level-Set (LS)

L'idée derrière la méthode Level-Set (LS) est de pouvoir accéder facilement et avec précision aux propriétés de l'interface $\partial\Omega$, à savoir la courbure $\kappa$ et la normale $n$. Introduite par OSHER
<figure><img src="image_8.png" /><figcaption>FIGURE 1.4 – Isosurfaces de la fonction Level-set $\Phi$. L'interface est située à $\Phi = 0$. Le signe de $\Phi$ permet de connaître la phase où on se situe.</figcaption></figure>

et al. [85], elle fait partie des méthodes de capture d'interface car la position de l'interface est connue de manière implicite à travers une fonction $\Phi$. Cette dernière est une distance signée d'un point $\bar{x} \in \Omega$ à l'interface $\partial\Omega$. Elle prend la valeur 0 au niveau de l'interface et son signe positif ou négatif détermine la phase dans laquelle on se situe (Fig. 1.4). Ainsi, la fonction $\Phi(\bar{x}, t)$ peut être définir par la relation (1.27).

$$ \Phi(\bar{x}, t) = \begin{cases} = 0, & \bar{x} \in \partial\Omega \\ > 0, & \bar{x} \in \text{phase 1} \\ < 0, & \bar{x} \in \text{phase 2} \end{cases} \qquad (1.27) $$

L'évolution de l'interface s'effectue grâce à l'advection de la fonction $\Phi$ (1.28).

$$ \frac{\partial}{\partial t}(\Phi) + \mathbf{u} \cdot \nabla \Phi = 0 \qquad (1.28) $$

La résolution numérique de l'équation (1.28) pose certains problèmes car au cours de sa résolution, il n'y a aucune garantie que la fonction $\Phi$ demeure une fonction distance signée. Pour y remédier, SUSSMAN et al. [115] proposent une procédure de réinitialisation de la fonction $\Phi$ après chaque advection. Par la suite, OLSSON et al. [83] [84] apportent une amélioration à cette réinitialisation et en profite pour proposer une approche conservative CLS (Conservative Level-Set). En effet, contrairement à la méthode VOF, la méthode LS peine à conserver le volume. Concernant les schémas de discrétisation des équations d'advection et de réinitialisation de la fonction $\Phi$, nous guidons le lecteur vers les travaux de [39] pour plus de détails. Mise à part les difficultés pour conserver le volume, la méthode LS gère aisément les changements topologiques de l'interface et accède facilement à la courbure $\kappa$ et la normale $\mathbf{n}$ par les relations (1.29).

$$ \kappa = -\nabla \cdot \mathbf{n} \quad \text{avec} \quad \mathbf{n} = \nabla \Phi / \|\nabla \Phi\| \qquad (1.29) $$

### 1.2.6 La méthode Front-Tracking (FT)

L'approche Front-Tracking (FT), dans la catégorie des méthodes de suivi d'interface, offre la connaissance explicite de la position de l'interface. En effet, en plus de la grille eulérienne fixe pour la résolution des équations de conservation, elle dispose d'un maillage distinct dédié à l'interface. La figure 1.5 en donne une illustration. Les premiers travaux sur la méthode FT remontent à [40, 43, 53, 124]. L'évolution de l'interface se fait à travers le transport lagrangien des
<figure><img src="image_9.png" /><figcaption>FIGURE 1.5 – Bulle/goutte modélisée par la méthode Front-Tracking (FT). Des marqueurs recouvrent la surface et forment des triangles dont la taille moyenne des arêtes est d. La grille eulérienne a une résolution d'espace h.</figcaption></figure>

points ou marqueurs de l'interface connaissant la vitesse de cette dernière. Ainsi pour un mar-
queur $m \in \partial\Omega$ de l'interface et de coordonnées $\bar{x}_m$, son mouvement est régi par l'équation (1.30),
avec $\mathbf{V}_i$ la vitesse de l'interface.

$$
\frac{d}{dt}(\bar{x}_m) = \mathbf{V}_i \tag{1.30}
$$

Lors de leur mouvement, les marqueurs du maillage de l'interface peuvent s'éloigner ou se rap-
procher, conduisant à une mauvaise représentation de l'interface si rien n'est fait. Pour ce faire,
tout un ensemble de procédures de gestion du maillage de l'interface doivent être mis en place.
Des procédures pour les phénomènes de coalescence ou de rupture doivent être également prises
en compte car l'approche FT ne gère pas automatiquement ces processus. Deux philosophies
existent quant à la gestion des marqueurs de l'interface. D'un côté, il y a l'approche où les
marqueurs de l'interface sont connectés entre eux. Nous pouvons citer dans cette catégorie les
travaux de [26, 42, 122]. De l'autre côté, TORRES et al. [120] introduisent la méthode PSM
(Point-Set Method). Nous avons également la méthode LCRM (Level Contour Reconstruction
Method) développée par SHIN et al. [108] [107]. Une amélioration, la méthode LFRM (Local
Front Reconstruction Method), en est proposée par [110]. Bien que la méthode FT a l'avantage
de connaître explicitement la position de l'interface et d'en déduire ses propriétés géométriques
grâce à la géométrie différentielle, elle a l'inconvénient de ne pas conserver le volume. GLIMM et
al. [41] se sont penchés sur ce problème et proposent une approche conservative pour la méthode
FT.

1.2.7 Les méthodes mixtes

Les différentes approches précédemment présentées possèdent toutes leurs avantages et leurs
inconvénients. Les méthodes VOF, MOF et LS gèrent automatiquement la coalescence et la
rupture mais moyennent la connaissance de l'interface à l'intérieur d'une maille eulérienne, ce
qui n'est pas le cas de la méthode FT. La conservation de la masse est assurée pour les méthodes
VOF et MOF alors que les approches LS et FT ne conservent pas la masse. Ces dernières offrent
de meilleures résultats dans l'évaluation des propriétés géométriques de l'interface, contrairement
aux approches VOF et MOF. Au vu de ces éléments, différents auteurs choisissent de combiner
les méthodes entre elles afin de profiter des avantages de chacune d'entre elles. Ainsi, SUSSMAN
et al. [114] proposent le couplage entre LS et VOF avec la méthode CLSVOF (Coupled Level-Set/Volume-Of-Fluid) dans l'idée d'améliorer la conservation de la masse et de mieux évaluer les tensions de surfaces. MÉNARD et al. [75] [62] utilisent une approche LS/VOF/GFM (Ghost Fluid Method par [34]) pour investiguer l'atomisation primaire d'un jet liquide. SUN et al. [113] développent une version simplifiée du couplage VOF/LS, appelée VOSET, où le calcul de la fonction LS se fait au moyen d'un processus géométrique itératif. La méthode LCRM introduite TORRES et al. [120] est une combinaison des méthodes FT et LS.

## 1.3 De l'utilité d'une bonne évaluation des données géométriques de l'interface

Différentes méthodes pour modéliser l'interface et son évolution au cours du temps ont été présentées. Que ce soit par le transport d'une fraction volumique, d'une fonction distance ou de marqueurs positionnés à sur l'interface, chacune de ces approches possède ses avantages et ses inconvénients. Toutefois, elles doivent toutes faire le lien entre la représentation qu'elles ont de l'interface et ses propriétés géométriques, à savoir la courbure $\kappa$ et la normale $\mathbf{n}$. En effet, ces éléments sont les principaux facteurs définissant les forces de tension de surface qui différencient un écoulement monophasique d'un écoulement diphasique. Par conséquent, une bonne modélisation et évaluation des forces de tension de surface, et par extension de la courbure et de la normale, semble crucial dans la simulation d'écoulements diphasiques. De même qu'il existe une multiplicité dans la modélisation de l'interface, une grande variété de méthodes d'évaluation des courbures et normales de l'interface a été développée au cours des années. Après avoir présenté les formulations employées pour décrire les forces de tension de surface dans les équations de conservation, quelques approches pour évaluer la courbure et la normale sont exposées. Elles font par la suite place à des exemples illustrant l'importance du traitement à accorder à l'évaluation de la courbure et de la normale.

### 1.3.1 Formulations des forces de tensions de surface

Dans le modèle 2-fluide, les forces de tensions de surface apparaissent dans les relations de saut de contraintes (1.8) ou dans l'équation de conservation de la quantité de mouvement pour le modèle 1-fluide (1.23). Dans cette dernière, la force de tension de surface par unité de volume est donnée par le terme $\mathbf{f}_{\sigma} = \sigma \kappa \mathbf{n} \delta_{i}$. Intégrée sur un volume de contrôle $\delta\Omega$, la force de tension de surface $\mathbf{F}_{\sigma}$ s'exprime par la relation (1.31).

$$ \mathbf{F}_{\sigma} = \int_{\delta\Omega} \mathbf{f}_{\sigma} d\mathcal{V} \quad \text{avec} \quad \mathbf{f}_{\sigma} = \sigma \kappa \mathbf{n} \delta_{i} \qquad (1.31) $$

A partir de cette relation, deux formulations sont possibles.

- Une approche surfacique : la force de tension de surface agissant sur l'interface uniquement, l'approche naturelle est de l'exprimer à cet endroit. La distribution de Dirac $\delta_{i}$ joue ce rôle car elle est uniquement non-nulle sur l'interface et permet ainsi de transformer la relation volumique (1.31) en une intégrale de surface (1.32). La surface de l'interface est notée $S$.

$$ \mathbf{F}_{\sigma} = \int_{S} \sigma \kappa \mathbf{n} \, d\mathcal{S} \qquad (1.32) $$

A ce stade, la force de tension de surface $\mathbf{F}_{\sigma}$ est évaluée soit géométriquement en utilisant les relations de Frenet [122], soit par un calcul direct [117]. En effet, dans le premier cas,
les relations de Frenet (en 2D) lient entre elles la courbure $κ$, la normale $n$, la tangente $t$ et la variation $dS$ sur une ligne courbe (1.33).

$$dt = \kappa ndS \quad \text{et} \quad d\mathbf{n} = -\kappa t dS \qquad (1.33)$$

En considérant dans le cas 2D un contour $Γ$ délimité par les points A et B, la force de tension de surface est donné simplement par la relation (1.34).

$$\mathbf{F}_{\sigma} = \int_{\Gamma} \sigma dt = \sigma_{B}t_{B} - \sigma_{A}t_{A} \qquad (1.34)$$

Ainsi la force de tension de surface s'obtient simplement par le calcul des tangentes à la surface. Dans le cas 3D, le terme $κn$ est plutôt remplacé par $(n \times \nabla) \times n$ (démonstration dans [123]) et permet d'aboutir grâce au théorème de Stokes à la relation (1.35).

$$\mathbf{F}_{\sigma} = \int_{S} \sigma \mathbf{t} \times \mathbf{n} \, dS \qquad (1.35)$$

Comme les équations de conservation sont résolues sur une grille fixe eulérienne, certaines approches transfèrent les informations de la surface vers la grille eulérienne et/ou vice-versa. TRYGGVASON et al. [122] utilisent ainsi des fonctions de PESKIN [87] pour interpoler les forces de tension de surface sur la grille eulérienne.

**Une approche volumique** : Un des éléments principaux dans l'approche volumique est le fait que la fonction Dirac $\delta_i$ est la dérivée au sens des distributions de la fonction Heaviside $H$. La force de tension de surface (1.31) se réécrit donc pour donner la relation (1.36).

$$\mathbf{F}_{\sigma} = \int_{\delta\Omega} \mathbf{f}_{\sigma} \, dV \quad \text{avec} \quad \mathbf{f}_{\sigma} = \sigma\kappa\nabla H \qquad (1.36)$$

L'étape suivante consiste à trouver une approximation $H_\epsilon$ de la fonction Heaviside. Le paramètre $\epsilon$ correspond à l'épaisseur donnée à l'interface. Aussi, lorsque $\epsilon$ tend vers 0, la fonction Heaviside est retrouvée (1.37)

$$\lim_{\epsilon \to 0} H_{\epsilon} = H \qquad (1.37)$$

Le tableau 1.1 donne quelques exemples d'approximations de la fonction Heaviside [93]. Ainsi dans une modélisation VOF de l'interface, la fonction couleur $C$ est utilisée comme une approximation de la fonction Heaviside. Cela s'apparente à une approche Continuum Surface Force (CSF) de BRACKBILL et al. [10] où l'épaisseur de l'interface est de l'ordre de la maille de la grille eulérienne. Dans une approche Level-Set, l'approximation $H_\epsilon$ est construite à partir de la fonction distance signée $\Phi$.

Ces différentes approches sont plus au moins adaptées à la modélisation choisie pour décrire l'interface. Ainsi pour une modélisation VOF ou LS, l'approche volumique est plus naturelle car la construction de l'approximation de la fonction Heaviside est aisée grâce à l'utilisation de la fonction couleur $C$ (VOF) ou la fonction distance signée $\Phi$ (LS). La modélisation FT permettant l'accès direct à la position de l'interface est plus encline à une approche surfacique. Toutefois cela
<table><caption>TABLE 1.1 – Exemples d'approximations de la fonction Heaviside [93].</caption><thead><tr><th>Méthode</th><th>Fonction H<sub>ε</sub></th></tr></thead><tbody><tr><td>CSF (VOF)</td><td>C</td></tr><tr><td>smooth (Level-Set/Front-Tracking)</td><td>$$ \begin{cases} 0 & , \Phi < -\epsilon \\ 1 & , \Phi > +\epsilon \\ (1 + \Phi/\epsilon + \sin(\pi\Phi/\epsilon)/\pi)/2 & , \sin\epsilon \end{cases} $$</td></tr><tr><td>GFM (Level-Set/VOF)</td><td>$$ \begin{cases} 1 & \text{si } \Phi > 1 \text{ (Level-Set) ou } C > 0.5 \text{ (VOF)} \\ 0 & \text{sinon} \end{cases} $$</td></tr></tbody></table>

n'empêche pas de combiner les approches surfacique et volumique entre elles. Ainsi SHIN et al. [107] proposent une approche hybride dans laquelle la force volumique $f_σ$ de tension de surface est exprimée dans une approche volumique avec une courbure calculée à partir d'une formulation surfacique. Plus de détails sur les calculs pour cette formulation hybride sont données dans [106, 109].

Par ailleurs, une attention particulière doit être apportée lors de la discrétisation des termes de tensions de surface et de gradient de pression. En effet, dans le cas d'une bulle/goutte statique ($\mathbf{u} = 0$) ou en translation ($\mathbf{u} = cste$) et en l'absence de gravité, les termes de tensions de surface équilibrent le gradient de pression (1.38).

$$ -\nabla p + \sigma \kappa \mathbf{n} \delta_i = -\nabla p + \sigma \kappa \nabla H = 0 \qquad (1.38) $$

Par conséquent, afin d'assurer cet équilibre, les schémas de discrétisation de la pression et des forces de tension de surface (ou de la fonction Heaviside) doivent être consistants et la courbure doit rester constante sur cette sphère. La rupture de cette équilibre conduit à l'apparition des courants parasites.

Quelle que soit l'approche employée pour formuler les forces de tensions de surface, il n'en demeure pas moins que la courbure $\kappa$ et la normale $\mathbf{n}$ à l'interface doivent être évaluées.

### 1.3.2 Estimation de la courbure et la normale

La courbure $\kappa$ et la normale $\mathbf{n}$ sont des propriétés géométriques de l'interface et interviennent dans le calcul des forces de tensions de surface. On peut ainsi distinguer les méthodes utilisant l'information implicite de l'interface pour en déduire ses propriétés et celles qui se basent sur des éléments de géométrie différentielles car elles ont accès à la position explicite de la surface.

Une des méthodes les plus classiques d'évaluation de la courbure et la normale à l'interface est la méthode Continuum Surface Force (CSF). Elle fut introduite par BRACKBILL et al. [10] dans le cadre de la méthode VOF. La fonction couleur $C$ (ou une version lissée $\tilde{C}$) est utilisée pour évaluer la courbure et la normale (1.39)

$$ \kappa = -\nabla \cdot \mathbf{n} \quad \text{avec} \quad \mathbf{n} = \nabla C / \|\nabla C\| \qquad (1.39) $$

Cependant avec cette formulation, les erreurs sur la courbure $\kappa$ tendent à augmenter avec la résolution spatiale. Bien que initialement développée dans le cadre de la méthode VOF, l'approche CSF peut se décliner à toute autre méthode de modélisation de l'interface, du moment où une fonction caractéristique peut être construite.

Une alternative est la méthode Continuous Surface Stress (CSS) de SCARDOVELLI et al. [104]. Elle permet de calculer les forces de tension de surface sans passer par l'évaluation de la courbure. Elle introduit le tenseur $T$ défini par la relation (1.40).
$$T = -\sigma (\mathbb{I} - \mathbf{n} \otimes \mathbf{n}) \delta_i \qquad (1.40)$$

avec $\mathbb{I}$ la matrice identité. La normale $\mathbf{n}$ est évaluée de la même manière que dans la relation (1.39) et à partir d'une fonction couleur lissée $\tilde{C}$. La force de tension de surface $\mathbf{f}_\sigma$ est calculée grâce à la relation (1.41).

$$\mathbf{f}_{\sigma} = \sigma \kappa \mathbf{n} \delta_{i} = -\nabla \cdot T \qquad (1.41)$$

On aboutit ainsi à une formulation conservative des forces de tension de surface.

Toujours dans le cadre de la modélisation VOF, la courbure peut être évaluée avec précision en construisant une fonction hauteur $y = f_y(x)$. La courbure $\kappa$ est calculée facilement grâce à l'équation (1.42).

$$\kappa = \frac{f_y''}{(1 + (f_y')^2)^{3/2}} \qquad (1.42)$$

La fonction hauteur $y = f(x)$ se construit en sommant dans une direction les valeurs de la fonction couleur $C$. Si ce choix de construction ne convient pas, l'orientation peut être changée pour construire cette fois la fonction hauteur dans l'autre direction $x = f_x(y)$. Afin de gagner en précision sur l'évaluation de $\kappa$, il suffit dans l'équation (1.42) de choisir des discrétisations d'ordre élevé pour le calcul des dérivées de la fonction hauteur $f_y$.

En ce qui concerne la modélisation FT, la connaissance explicite de la position de l'interface permet l'usage d'éléments de géométrie différentielle pour estimer avec précision les courbures et les normales de l'interface. Ainsi, MEYER et al. [76] utilise l'opérateur de Laplace-Beltrami pour évaluer la courbure et la normale sur un maillage de triangles. Cet opérateur qui est un laplacien surfacique, est appliqué à un point $\bar{x} \in \partial\Omega$ de l'interface puis intégré sur une région autour de ce point. La direction du vecteur résultant fournit celle de la normale et son intensité est liée à la courbure. La connaissance explicite de la position de l'interface permet également d'en construire une équation. L'approche locale est la plus commune [26]. L'approximation locale de la surface est établie généralement en minimisant un résidu rms dont la forme est donnée par la relation (1.43).

$$rms = \sum_k \omega_k (\bar{x}_k \cdot \mathbf{M}\bar{x}_k + \mathbf{n} \cdot \bar{x}_k + b)^2 \quad \text{et} \quad \nabla(rms) = 0 \qquad (1.43)$$

avec $\bar{x}_k$ les points utilisés pour effectuer les moindre-carrés, $\mathbf{n}$ une approximation de la normale, $\omega_i$ des poids judicieusement sélectionnés pour rendre compte de l'influence des points dans le calcul de l'équation locale de la surface et $\mathbf{M}$ la matrice symétrique de la forme quadratique de la surface. Si la surface a pour équation générique $Ax^2 + By^2 + Cz^2 + 2Dyz + 2Ezx + 2Fxy + 2Gx + 2Hy + 2Iz + J$, la matrice $\mathbf{M}$ de la forme quadratique (uniquement les termes $x^my^n$ de degré $m+n=2$) est donnée par la relation (1.44).

$$\mathbf{M} = \begin{pmatrix} A & F & E \\ F & B & D \\ E & D & C \end{pmatrix} \qquad (1.44)$$

Un nombre suffisant de points est choisi pour rendre le système défini en (1.43) inversible. Les inconnues sont les coefficients de l'équation de la surface. La courbure $\kappa$ peut être évaluée en
utilisant le gradient et la hessienne de l'équation de la surface [45]. Il est à noter que reconstruire localement une interface n'est pas seulement restreint à la méthode FT. Ainsi, RENARDY et al. [98] développent la méthode PROST (Parabolic Reconstruction Of Surface Tension) dans laquelle la minimisation est opérée entre les fractions volumiques C issues de l'advection de l'interface et celles obtenues en se servant de l'équation locale de l'interface.

Cette présentation montre les différentes approches développées pour traiter des forces de tensions de surface et de l'évaluation de leurs paramètres que sont la courbure et la normale dans la modélisation et la simulation des écoulements. Afin de montrer la validité et l'utilité de ces méthodes, des exemples sont analysés dans la section suivante.

### 1.3.3 Exemples de l'utilité des approches de traitement des forces de tension de surface et de l'évaluation de la courbure et la normale

Différentes techniques sont employées pour décrire les forces de tensions de surface et calculer les courbures et les normales. Toutefois, pour vérifier leur validité et fiabilité à la réalité, elles doivent être testées sur des cas simples. Un de ces cas tests est celui de la bulle/goutte statique (ou en translation) qui met en lumière l'efficacité des schémas de traitement des forces de tensions de surface et des méthodes de calcul de la courbure et la normale. Les propriétés géométriques de l'interface ne servent pas seulement à estimer les forces de tensions de surface, mais également à la description de phénomène d'atomisation. Une approche utilisant ces informations sera présentée.

#### Cas test de la bulle/goutte statique

Cette expérience permet de vérifier si les approches de calculs de la courbure et la normale, de même que les techniques de discrétisation des forces de tensions de surface et des gradients de pression respectent la loi de Laplace. Cette dernière se résume en l'équilibre entre le gradient de pression et les forces de tensions de surface pour une bulle/goutte statique (1.38). La proposition reste également vraie pour une bulle/goutte en translation. Si l'opérateur rotationnel est appliqué à l'équation (1.38), la condition se transforme et donne la relation (1.45).

$$ (\nabla \kappa) \times (\nabla H) = 0 \qquad (1.45) $$

Ainsi pour une bulle/goutte, cette relation est respectée si la courbure reste constante. En effet, la variation de la fonction Heaviside $H$ est dans le sens de la normale à l'interface. Quant à celle de la courbure $\kappa$, elle se fait tangentiellement à la surface. Par conséquent, le seul moyen de vérifier la condition (1.45) est d'avoir $\nabla \kappa = 0$ donc $\kappa = \text{cste}$. Deux éléments sont donc à vérifier pour la loi de Laplace : la consistance des schémas de discrétisation de forces de tensions de surface et du gradient de pression d'une part, et d'autre part la précision de l'évaluation de la courbure. C'est la violation de ces conditions qui conduit à l'apparition de courants parasites. Plusieurs auteurs se sont penchés sur ce test [1, 32, 92, 117]. Nous présentons la configuration employée par [1] mais l'esprit reste le même dans les autres approches.

Le test (2D) consiste en un cercle de rayon $D = 0.8$ placé dans un domaine carrée. La symétrie du problème fait qu'on peut se restreindre qu'à la simulation d'un quart du cercle, placé dans un coin (ici en bas à gauche) d'un domaine carré de côté 1. Les conditions limites sont de type symétrie sur les côtés du carré touchés par le cercle et de non glissement sur les deux autres. Les deux phases ont les mêmes densité $\rho$ et viscosité $\mu$. Les paramètres du problème sont le nombre capillaire $Ca_{max}$ (1.46) permettant de mesurer l'intensité des courants parasites et le nombre de Laplace $La$ (1.47).

$$ Ca_{max} = \frac{\| \mathbf{u}_{max} \| \mu}{\sigma} \qquad (1.46) $$
$$La = \frac{\sigma \rho D}{\mu^2} \qquad (1.47)$$

A cela, il faut ajouter deux temps caractéristiques $T_σ$ (1.48) et $T_μ$ (1.49).

$$T_{\sigma} = \sqrt{\frac{\rho D^3}{\sigma}} \qquad (1.48)$$

$$T_{\mu} = \frac{\rho D^2}{\mu} \qquad (1.49)$$

Ils correspondent respectivement au temps caractéristique des ondes capillaires et à celui de la dissipation visqueuse. Ces deux derniers paramètres sont d'importance car il faut que la simulation soit effectuée pendant un temps suffisamment long ($t > T_σ$ et $t > T_μ$) pour arriver à un état d'équilibre [92, 93]. Sans quoi, la simulation est encore dans un régime transitoire où les courants parasites sont dus plutôt à l'évolution d'une surface initialement perturbée. Ces perturbations sont amorties sur des temps de l'ordre de $T_μ$.

<figure><img src="image_10.png" /><figcaption>FIGURE 1.6 – Évolution du maximum d'intensité $Ca_{max}$ des courants parasites dans le domaine de simulation au cours du temps (Fig.3 de [1]).</figcaption></figure>

Dans les travaux de ABADIE et al. [1], différentes méthodes VOF et LS sont testées. Dans un premier temps, l'équation (1.38) est résolue avec une courbure fixée et constante dans le but de vérifier la consistance des schémas de discrétisation. Le nombre capillaire atteint alors la valeur $Ca_{max} = 5.08 \times 10^{-18}$, dénotant le respect de la loi de Laplace pour les schémas de discrétisation employés. Par conséquent, dans le cas d'une courbure évaluée numériquement, les courants parasites observés seront dus à la méthode utilisée pour calculer la courbure et à l'existence d'un gradient de courbure à l'interface de la bulle. Différentes techniques d'évaluation de la courbure sont ensuite comparées et les résultats montrant l'évolution de $Ca_{max}$ en fonction du temps adimensionnel $t^*_σ = t/T_σ$ sont répertoriés sur la figure 1.6. Dans le cadre des méthodes VOF, les techniques de calcul de courbure utilisant une approche de type CSF n'arrivent pas à réduire les courants parasites. Par contre, l'emploi d'une fonction hauteur permet de les réduire grandement, avec l'approche VOF-PLIC-HFCSF offrant l'erreur la plus faible dans cette catégorie (Fig. 1.6). L'approche CSF combinée avec une méthode LS (LS-CCSF) offre de mesures résultats que sa contrepartie VOF. Toutefois, elle ne converge pas, de même que les techniques LS-HFCSF et
LS-SSF. Cependant, ce sont ces deux dernières techniques qui arrivent à réduire les courants parasites à l'ordre de l'erreur machine lorsque l'étape de réinitialisation de la fonction distance n'est pas opérée. L'auteur attribue ce gain de performance dans cette configuration au fait que l'étape de réinitialisation de la fonction distance change la position de l'interface, altérant ainsi le calcul de la courbure. Ces conclusions sur l'impact d'une approche CSF quant aux courants parasites se retrouvent dans d'autres travaux avec une modélisation FT [32, 117] ou VOF [92]. L'idée principale à retenir étant de découpler le calcul de la courbure du gradient de la fonction couleur C, et que ce calcul soit le plus localisé possible sur la surface.

Caractérisation de l'évolution d'une atomisation à partir des données géométriques de l'interface

Les propriétés géométriques de l'interface sont également utiles dans d'autres domaines tels que l'atomisation de jets. Elle met en jeu une multitude d'échelles d'interface rencontrées dans différentes régions du jet atomisé. A la base du jet, les deux phases sont clairement distinctes : on parle alors de régime à phases séparées. Loin de la base du jet où toute une population de gouttes évolue constitue le régime à phase dispersée. La transition entre les deux est le lieu où les deux échelles coexistent, avec la présence de ligaments et gouttes de tailles moyennes. Tandis que le régime à phases séparées peut être résolu avec une approche 2-fluides ou 1-fluide, cela peut être prohibitif dans le régime à phase dispersée à cause de la grande quantité de gouttes/bulles présentes. Une des approches pour simuler ce régime est l'approche probabiliste qui se base sur la Number Density Function (NDF). Elle permet de suivre l'évolution de la population de gouttes grâce à l'équation de Williams-Boltzmann et d'avoir accès à des grandeurs comme la vitesse, la taille, la température pour ne citer que celles-ci. Toutefois, NDF n'est définie que dans le régime à phase dispersée et n'apporte pas d'informations sur le régime à phases séparées. En effet, elle dénombre le nombre de gouttes remplissant certains critères choisis par l'utilisateur, gouttes qui ne se sont pas encore formées dans le régime à phases séparées. Massot et al. [13, 77] introduisent alors la Surface Density Function (SDF) pour d'écrire l'atomisation tout au long de son déroulement. Comme son nom l'indique, elle rend compte de l'aire présentant une caractéristique donnée (courbure moyenne $κ$, courbure de Gauss $κ_G$, vitesse de l'interface). Contrairement à NDF, elle est définie tout le long du processus d'atomisation et elle permet ainsi de suivre l'atomisation du régime à phases séparées à celui à phase dispersée. Par ailleurs, un lien est établi entre SDF et NDF grâce au théorème de Gauss-Bonnet qui lie l'intégrale de la courbure de Gauss sur la surface d'un objet à sa caractéristique d'Euler (elle vaut 2 dans le cas d'une sphère). Comme la SDF utilise les données géométriques de l'interface, l'attention quant à l'évaluation de ces données est un élément important. A ce titre, sur un maillage des triangles de l'interface, c'est la méthode de [76] dans sa version barycentrique qui est employée.

1.4 Plan de la thèse

Nous avons vu au travers de cette partie la richesse, la diversité et la complexité des écoulements diphasiques à phases séparées, éléments qui se reflètent également dans les approches pour les modéliser. Quel que soit la méthode utilisée, les forces de tensions de surfaces demeurent un élément clé dans la simulation d'écoulements diphasiques. Elles font intervenir des propriétés géométriques de l'interface que sont la courbure $κ$ et la normale $n$ à l'interface. Leur bonne évaluation est primordiale tant pour éviter la création de courants parasites, que dans l'établissement de grandeurs macroscopiques (PDF, aire interfaciale) utiles dans des domaines tels que la combustion, l'agroalimentaire ou le milieu pharmaceutique. Ces grandeurs macroscopiques ne fournissent cependant pas suffisamment d'information sur les types de gouttes/bulles rencontrés dans les écoulements diphasiques. C'est dans ce cadre que s'inscrit ce travail dont la finalité est l'établissement d'une base de données permettant la caractérisation et l'identification automatique des courbures et des normales aux interfaces d'écoulements diphasiques. Pour ce faire,
nous avons agencés notre travail comme suit.

* Notre élément principal de travail étant l'interface, l'objectif est d'avoir dans le maximum d'information sur elle. Notre choix s'est porté donc sur le développement d'une méthode Front-Tracking 3D qui à l'avantage de connaître explicitement la position de l'interface. Ainsi le chapitre 2 est consacré à la présentation de cette méthode. Afin de suivre l'évolution topologique et complexe de l'interface, différentes procédures de gestion tant au niveau de la richesse et de la qualité du maillage de l'interface sont présentées. Une procédure de contrôle du volume est également détaillée. Elles ont été entièrement développées et implémentées dans le code de simulation maison FUGU de l'équipe TCM. L'ensemble des ces procédures est ensuite analysé dans un test impliquant la déformation d'une goutte/bulle dans un champ de vitesse tourbillonnaire.

* Le chapitre 3 est méthodes permettant l'évaluation des propriétés de l'interface dans le cadre de l'approche FT. Trois familles de méthodes sont présentées et comparées entre elles pour déterminer leur précision mais également leur robustesse quant à la qualité du maillage de l'interface. Deux tests constituent le creuset de leur comparaison : un test statique employant une surface qui présente une large variétés de type de courbure, et un test original dynamique qui permet de faire le lien avec des phénomènes réels tels que l'impact d'une goutte sur une surface ou la rupture d'une goutte. Ce chapitre fait l'objet d'un article en cours de révision dans *Journal of Computational Physics*.

* Le chapitre 4 est consacré à l'introduction et à l'étude d'une grandeur utile pour la caractérisation et l'identification de type de gouttes/bulles comme rencontrés dans le diagramme de Clift. Cette grandeur est analysée puis testée sur des structures simples, à commencer par des sphères et des ellipsoïdes, avant de s'orienter vers des structures complexes de type wobbling. Ce type est modélisé dans ce travail par des perturbations sinueuses et variqueuses dont une étude en fonction des paramètres de la perturbation (amplitude, longueur d'onde, déphasage) est effectuée.

Les travaux entrepris dans ce travail ont permis la publication de quatre articles et la mise en
place d'une méthode 3D de suivi d'interface dans le code maison FUGU développé dans l'équipe
TCM du laboratoire MSME.

— D.-A. KOFFI BI et al. "Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces". In: submitted in *J. Comput. Phys.* (2021)

— D.-A. KOFFI BI et al. “A Review of Geometrical Interface Properties for 3D Front-Tracking Methods”. In : *Turbulence and Interactions*. Sous la dir. de M. DEVILLE et al. Cham : Springer International Publishing, 2021, p. 144-149

— M. TAVARES et al. “A Front-Tracking Method for Multiphase Flows with a Sharp Interface Representation”. In : *Turbulence and Interactions*. Sous la dir. de M. DEVILLE et al. Cham : Springer International Publishing, 2021, p. 189-195

— M. TAVARES et al. “A Two-Dimensional Second Order Conservative Front-Tracking Method with an Original Marker Advection Approach Based on Jump Relations”. In : *Communications in Computational Physics* 27.5 (2020), p. 1550-1589

— Implémentation d'une méthode 3D de Front-Tracking dans le logiciel-maison FUGU,
développé par l'équipe Transferts de Chaleur et de Matière (TCM) du laboratoire Mo-délisation et Simulation Multi Échelle (MSME UMR 8208)
# Bibliographie

## Introduction

[1] T. ABADIE, J. AUBIN et D. LEGENDRE. "On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks". In: *J. Comput. Phys.* 297 (2015), p. 611-636 (cf. p. 34, 35).

[3] H. T. AHN et M. SHASHKOV. "Multi-material interface reconstruction on generalized polyhedral meshes". In: *J. Comput. Phys.* 226 (2007), p. 2096-2132 (cf. p. 27).

[4] S. L. ANNA. "Droplets and Bubbles in Microfluidic Devices". In: *Annu. Rev. Fluid Mech.* 48 (2016), p. 285-309 (cf. p. 20, 21).

[5] M. ATAEI, M. BUSSMANN, V. SHAAYEGAN, F. COSTA, S. HAN et C. B. PARK. "NPLIC : A machine learning approach to piecewise linear interface construction". In : *Computers and Fluids* 223 (2021), p. 104950 (cf. p. 27).

[6] E. AULISA, S. M. and R. SCARDOVELLI et S. ZALESKI. "Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry". In : *J. Comput. Phys.* 225 (2007), p. 2301-2319 (cf. p. 27).

[7] E. BABINSKY et P. SOJKA. "Modeling drop size distributions". In: *Prog. Energy Combust. Sci.* 28 (2002), p. 303-329 (cf. p. 20).

[8] S. BALACHANDAR et J. K. EATON. "Turbulent Dispersed Multiphase Flow". In: *Annu. Rev. Fluid Mech.* 42 (2010), p. 111-133 (cf. p. 20).

[10] J. BRACKBILL, D. KOTHE et C. ZEMACH. "A continuum method for modeling surface tension". In: *J. Comput. Phys.* 100.2 (1992), p. 335-354 (cf. p. 31, 32, 73).

[11] N. BREMOND et E. VILLERMAUX. "Atomization by jet impact". In: *J. Fluid Mech.* 549 (2006), p. 273-306 (cf. p. 20).

[12] C. E. BRENNEN. *Fundamentals of Multiphase Flows*. Cambridge University Press, 2005 (cf. p. 21).

[13] R. CANU, S. PUGGELLI, M. ESSADKI, B. DURET, T. MÉNARD, M. MASSOT, J. REVEILLON et F. DEMOULIN. "Where does the droplet size distribution come from ?" In: *Int. J. Multiphase Flow* 107 (2018), p. 230-245 (cf. p. 36, 72, 73).

[15] S. CHANDRA et C. T. AVEDISIAN. "On the Collision of a Droplet with a Solid Surface". In: *Proceedings : Mathematical and Physical Sciences* 432.1884 (jan. 1991), p. 13-41 (cf. p. 21).

[16] P.-C. CHEN, W.-C. WANG, W. L. ROBERTS et T. FANG. "Spray and atomization of diesel fuel and its alternatives from a single-hole injector using a common rail fuel injection system". In: *Fuel* 103 (2013), p. 850-861 (cf. p. 21).

[18] G. COMPÈRE, E. MARCHANDISE et J.-F. REMACLE. "Transient adaptivity applied to two-phase incompressible flows". In: *J. Comput. Phys.* 227 (2008), p. 1923-1942 (cf. p. 23).
[19] M. CUTFORTH, P. T. BARTON et N. NIKIFORAKIS. "An efficient moment-of-fluid interface tracking method". In : *Computers and Fluids* 224 (2021), p. 104964 (cf. p. 27).

[20] M. DAI et D. P. SCHMIDT. "Adaptive tetrahedral meshing in free-surface flow". In : *J. Comput. Phys.* 208 (2005), p. 228-252 (cf. p. 23).

[21] J. DELHAYE. "Jump conditions and entropy sources in two-phase systems. Local instant formulation". In : *Int. J. Multiphase Flow* 1.3 (1974), p. 395-409 (cf. p. 23, 72).

[26] J. DU, B. FIX, J. GLIMM, X. JIA, X. LI, Y. LI et L. WU. "A simple package for front tracking". In : *J. Comput. Phys.* 213 (2006), p. 613-628 (cf. p. 29, 33, 73, 78, 79, 83, 85, 90, 92, 107, 161, 163).

[27] L. DUCHEMIN, S. POPINET, C. JOSSERAND et S. ZALESKI. "Jet formation in bubbles bursting at a free surface". In : *Phys. Fluids* 14.9 (sept. 2002), p. 3000 (cf. p. 19).

[28] C. DUMOUCHEL. "On the experimental investigation on primary atomization of liquid streams". In : *Exp. Fluids* 45 (2008), p. 371-422 (cf. p. 20).

[29] V. DYADECHKO et M. SHASHKOV. "Moment-of-fluid interface reconstruction". In : *Tech. Rep. LA-UR-05-7571* (2005) (cf. p. 27).

[32] M. FEBRES et D. LEGENDRE. "Enhancement of a 2D front-tracking algorithm with a non-uniform distribution of Lagrangian markers". In : *J. Comput. Phys.* 358 (2018), p. 173-200 (cf. p. 34, 36).

[34] R. FEDKIW, T. ASLAM, B. MERRIMAN et S. OSHER. "A non-oscillatory Eulerian approach to interfaces in multi-material flows (the Ghost Fluid method)". In : *J. Comput. Phys.* 152 (1999). Sous la dir. de J. of COMPUTATIONAL PHYSICS, p. 457-492 (cf. p. 30, 73).

[35] D. FUSTER, G. AGBAGLAH, C. JOSSERAND, S. POPINET et S. ZALESKI. "Numerical simulation of droplets, bubbles and waves : state of the art". In : *Fluid Dyn. Res.* 41 (2009), p. 065001 (cf. p. 20).

[36] D. E. FYFE, E. S. ORAN et M. J. FRITTS. "Surface Tension and Viscosity with Lagrangian Hydrodynamics on a Triangular Mesh". In : *J. Comput. Phys.* 76 (1988), p. 349-384 (cf. p. 23).

[37] S. GANESAN et L. TOBISKA. "An accurate finite element scheme with moving meshes for computing 3D-axisymmetric interface flows". In : *Int. J. Numer. Methods Fluids* 57 (2008), p. 119-138 (cf. p. 23).

[38] S. GANESAN et L. TOBISKA. "Arbitrary Lagrangian-Eulerian finite-element method for computation of two-phase flows with soluble surfactants". In : *J. Comput. Phys.* 231 (2012), p. 3685-3702 (cf. p. 23).

[39] F. GIBOU, R. FEDKIW et S. OSHER. "A review of level-set methods and some recent applications". In : *J. Comput. Phys.* 353 (2018), p. 82-109 (cf. p. 28).

[40] J. GLIMM, M. GRAHAM, J. GROVE, X. LI, T. SMITH, D. TAN, F. TANGERMAN et Q. ZHANG. "Front tracking in two and three dimensions". In : *Computers & Mathematics with Applications* 35.7 (1998). Advanced Computing on Intel Architectures, p. 1-11 (cf. p. 28).

[41] J. GLIMM, X. LI, Y. LIU, Z. XU et N. ZHAO. "Conservative Front Tracking with Improved Accuracy". In : *SIAM J. Numer. Anal.* 41.5 (2003), p. 1926-1947 (cf. p. 29).

[42] J. GLIMM, J. W. GROVE, X. L. LI et D. C. TAN. "Robust Computational Algorithms for Dynamic Interface Tracking in Three Dimensions". In : *SIAM J. Sci. Comput.* 21.6 (2000), p. 2240-2256 (cf. p. 29).

[43] J. GLIMM, J. W. GROVE, X. L. LI, K.-M. SHYUE, Y. ZENG et Q. ZHANG. "Three-dimensional front tracking". In : *SIAM J. Sci. Comput.* 19.3 (1998), p. 703-727 (cf. p. 28).
[45] R. GOLDMAN. "Curvature formulas for implicit curves and surfaces". In: *Comput. Aided Geom. Des.* 22 (2005), p. 632-658 (cf. p. 34, 78, 79).

[46] S. GONG, W. MA et H. GU. "An experimental investigation on bubble dynamics and boiling crisis in liquid films". In : *Int. J. Heat Mass Transfer* 79 (2014), p. 694-703 (cf. p. 21).

[47] M. GOROKHOVSKI et M. HERRMANN. "Modeling Primary Atomization". In : *Annual Review of Fluid Mechanics* 40 (2008), p. 343-366 (cf. p. 20).

[48] R. GUILLAUMENT, S. VINCENT et J.-P. CALTAGIRONE. "An original algorithm for VOF based method to handle wetting effect in multiphase flow simulation". In : *Mech. Res. Commun.* 63 (2015), p. 26-32 (cf. p. 27).

[49] M. HERRMANN. "Detailed Numerical Simulations of the Primary Atomization of a Turbulent Liquid Jet in Crossflow". In : *J. Eng. Gas Turbines Power* 132 (juin 2010), p. 061506 (cf. p. 20).

[50] C. W. HIRT et B. D. NICHOLS. "Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries". In : *J. Comput. Phys.* 39 (1981), p. 201-225 (cf. p. 26).

[52] M. JEMISON, M. SUSSMAN et M. SHASHKOV. "Filament capturing with the Multimaterial Moment-of-Fluid method". In : *J. Comput. Phys.* 28 (2015), p. 149-172 (cf. p. 27).

[53] D. JURIC et G. TRYGGVASON. "Computations of boiling flows". In : *Int. J. Multiphase Flow* 24.3 (1998), p. 387-410 (cf. p. 28).

[54] I. KATAOKA. "Local instant formulation of two-phase flow". In : *Int. J. Multiphase Flow* 12.5 (1986), p. 745-758 (cf. p. 24).

[55] S. KESHANI, W. R. W. DAUD, M. NOUROUZI, F. NAMVAR et M. GHASEMI. "Spray drying : An overview on wall deposition, process and modeling". In : *J. Food Eng.* 146 (2015), p. 152-162 (cf. p. 21).

[56] D.-A. KOFFI BI, M. TAVARES, E. CHENIER et S. VINCENT. "Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces". In : *submitted in J. Comput. Phys.* (2021) (cf. p. 37).

[57] D.-A. KOFFI BI, M. TAVARES, É. CHÉNIER et S. VINCENT. "A Review of Geometrical Interface Properties for 3D Front-Tracking Methods". In : *Turbulence and Interactions*. Sous la dir. de M. DEVILLE, C. CALVIN, V. COUAILLIER, M. DE LA LLAVE PLATA, J.-L. ESTIVALÈZES, T. H. LÊ et S. VINCENT. Cham : Springer International Publishing, 2021, p. 144-149 (cf. p. 37).

[58] M. KUCHARIK, R. V. GARIMELLA, S. P. SCHOFIELD et M. J. SHASHKOV. "A comparative study of interface reconstruction methods for multi-material ALE simulations". In : *J. Comput. Phys.* 229 (2010), p. 2432-2452 (cf. p. 27).

[60] T. LANGRISH et D. FLETCHER. "Spray drying of food ingredients and applications of CFD in spray drying". In : *Chem. Eng. Process.* 40 (2001), p. 345-354 (cf. p. 21).

[61] J. C. LASHERAS, E. VILLERMAUX et E. J. HOPFINGER. "Break-up and atomization of a round water jet by a high-speed annular air jet". In : *J. Fluid Mech.* 357 (1998), p. 351-379 (cf. p. 20).

[62] R. LEBAS, T. MENARD, P. BEAU, A. BERLEMONT et F. DEMOULIN. "Numerical simulation of primary break-up and atomization : DNS and modelling study". In : *Int. J. Multiphase Flow* 35 (2009), p. 247-260 (cf. p. 20, 30).

[63] A. H. LEFEBVRE et V. G. McDONELL. *Atomization and Sprays*. Second. Combustion : An International Series. CRC Press, 2017 (cf. p. 21).

[64] D. LEGENDRE et J. MAGNAUDET. "The lift force on a spherical bubble in a viscous linear shear flow". In : *J. Fluid Mech.* 368 (1998), p. 81-126 (cf. p. 23).
[66] H. L'HUISSIER et E. VILLERMAUX. "Bursting bubble aerosols". In : *J. Fluid Mech.* 696 (2012), p. 5-44 (cf. p. 19).

[68] K. LUO, H. PITSCH, M. PAI et O. DESJARDINS. "Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor". In : *Proc. Combust. Inst.* 33 (2011), p. 2143-2152 (cf. p. 21).

[69] J. MAGNAUDET. *Méthodes numériques pour les problèmes d'inter interfaces*. Lalonde les Maures : Ecole de Printemps, mai 1997 (cf. p. 24).

[70] P. MARMOTTANT et E. VILLERMAUX. "On spray formation". In : *J. Fluid Mech.* 498 (2004), p. 73-111 (cf. p. 20).

[75] T. MÉNARD, S. TANGUY et A. BERLEMONT. "Coupling level set/VOF/ghost fluid methods : Validation and application to 3D simulation of the primary break-up of a liquid jet". In : *Int. J. Multiphase Flow* 33 (2007), p. 510-524 (cf. p. 30).

[76] M. MEYER, M. DESBRUN, P. SCHRODER et A. H. BARR. "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds". In : *Visualization and Mathematics III*. Sous la dir. de H.-C. HEGE et K. POLTHIER. Springer, 2003, p. 35-57 (cf. p. 33, 36, 75, 76, 81, 85, 90, 92, 107, 152, 154).

[77] E. MOHAMED, D. FLORENCE, de CHAISEMARTIN STÉPHANE, L. ADAM, M. THIBAULT et M. MARC. "Statistical modeling of the gas-liquid interface using geometrical variables : Toward a unified description of the disperse and separated phase flows". In : *Int. J. Multiphase Flow* 120 (2019), p. 103084 (cf. p. 36, 73).

[78] L. MORAWSKA, G. JOHNSON, Z. RISTOVSKI, M. HARGREAVES, K. MENGERSEN, S. C. B. C. CHAO, Y. LI et D. KATOSHEVSKI. "Size distribution and sites of origin of droplets expelled from the human respiratory tract during expiratory activities". In : *Aerosol Science* 40.3 (2009), p. 256-269 (cf. p. 21).

[79] C. MUNDO, M. SOMMERFELD et C. TROPEA. "Droplet-wall collisions : Experimental studies of the deformation and breakup process". In : *Int. J. Multiphase Flow* 21.2 (1995), p. 151-173 (cf. p. 21).

[80] N. NIKOLOPOULOS, A. THEODORAKAKOS et G. BERGELES. "A numerical investigation of the evaporation process of a liquid droplet impinging onto a hot substrate". In : *Int. J. Heat Mass Transfer* 50 (2007), p. 303-319 (cf. p. 21).

[81] W. F. NOH et P. WOODWARD. "SLIC (Simple line interface calculation)". In : *Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics June 28 – July 2, 1976 Twente University, Enschede*. Sous la dir. d'A. I. van de VOOREN et P. J. ZANDBERGEN. T. 59. Lecture Notes in Physics. Berlin, Heidelberg : Springer, Berlin, Heidelberg, juil. 1976 (cf. p. 26).

[82] S. NUKIYAMA. "The maximum and minimum values of the heat Q transmitted from metal to boiling water under atmospheric pressure." In : *Int. J. Heat Mass Transfer* 9 (1966), p. 1419-1433 (cf. p. 21).

[83] E. OLSSON et G. KREISS. "A conservative level set method for two phase flow". In : *J. Comput. Phys.* 210 (2005), p. 225-246 (cf. p. 28).

[84] E. OLSSON, G. KREISS et S. ZAHEDI. "A conservative level set method for two phase flow II". In : *J. Comput. Phys.* 225 (2007), p. 785-807 (cf. p. 28).

[85] S. OSHER et J. A. SETHIAN. "Front Propagating with curvature dependent speed : Algorithms based on Hamilton-Jacobi formulations". In : *J. Comput. Phys.* 79.1 (nov. 1988), p. 12-49 (cf. p. 27).

[86] A. PAUDEL, Z. A. WORKU, J. MEEUS, S. GUNS et G. V. den MOOTER. "Manufacturing of solid dispersions of poorly water soluble drugs by spray drying : Formulation and process considerations". In : *Int. J. Pharm.* (2012), p. 12732 (cf. p. 21).
[87] C. S. PESKIN. "Numerical analysis of blood flow in the heart". In: *J. Comput. Phys.* 25.3 (1977), p. 220-252 (cf. p. 31).

[88] M. PEZO et V. STEVANOVIC. "Numerical prediction of critical heat flux in pool boiling with the two-fluid model". In: *Int. J. Heat Mass Transfer* 54 (2011), p. 3296-3303 (cf. p. 21).

[89] J. E. PILLIOD. *An Analysis of Piecewise Linear Interface Reconstruction Algorithms for Volume-of-fluid Methods*. Davis: University of California, sept. 1992 (cf. p. 27).

[90] J. E. PILLIOD et E. G. PUCKETT. "Second-order accurate volume-of-fluid algorithms for tracking material interfaces". In: *J. Comput. Phys.* 199 (2004), p. 465-502 (cf. p. 27).

[92] S. POPINET. "An accurate adaptive solver for surface-tension-driven interfacial flows". In: *J. Comput. Phys.* 228.16 (2009), p. 5838-5866 (cf. p. 34-36).

[93] S. POPINET. "Numerical Models of Surface Tension". In: *Annu. Rev. Fluid Mech.* 50.1 (2018), p. 49-75 (cf. p. 31, 32, 35).

[96] E. G. PUCKETT. "A volume-of-fluid interface tracking algorithm with applications to computing shock wave refraction". In: *Proceedings of the Fourth International Symposium on Computational Fluid Dynamics*. Sous la dir. de H. DWYER. Davis, CA, 1991, p. 933-938 (cf. p. 27).

[98] Y. RENARDY et M. RENARDY. "PROST : A Parabolic Reconstruction of Surface Tension for the Volume-of-Fluid Method". In: *J. Comput. Phys.* 183 (2002), p. 400-421 (cf. p. 34).

[102] S. SAZHIN. *Droplets and Sprays*. Springer, 2014 (cf. p. 20).

[103] R. SCARDOVELLI et S. ZALESKI. "Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection". In: *Int. J. Numer. Methods Fluids* 41 (2003), p. 251-274 (cf. p. 27).

[104] SCARDOVELLI et ZALESKI. "Direct Numerical Simulation of free-surface and interfacial flow". In: *Annu. Rev. Fluid Mech.* 31 (jan. 1999), p. 567-603 (cf. p. 32).

[106] S. SHIN. "Computation of the curvature field in numerical simulation of multiphase flow". In: *J. Comput. Phys.* 222.2 (2007), p. 872-878 (cf. p. 32).

[107] S. SHIN, S. I. ABDEL-KHALIK, V. DARU et D. JURIC. "Accurate representation of surface tension using the level contour reconstruction method". In: *J. Comput. Phys.* 203 (2005), p. 493-516 (cf. p. 29, 32).

[108] S. SHIN et D. JURIC. "Modeling Three-Dimensional Multiphase Flow Using a Level Contour Reconstruction Method for Front Tracking without Connectivity". In: *J. Comput. Phys.* 180.2 (2002), p. 427-470 (cf. p. 29, 73).

[109] S. SHIN et D. JURIC. "A hybrid interface method for three-dimensional multiphase flows based on front tracking and level set techniques". In: *Int. J. Numer. Methods Fluids* 60.7 (2009), p. 753-778 (cf. p. 32).

[110] S. SHIN, I. YOON et D. JURIC. "The Local Front Reconstruction Method for direct simulation of two- and three-dimensional multiphase flows". In: *J. Comput. Phys.* 230.17 (2011), p. 6605-6646 (cf. p. 29).

[111] L. SHUI, J. C. EIJKEL et A. van den BERG. "Multiphase flow in microfluidic systems – Control and applications of droplets and interfaces". In: *Adv. Colloid Interface Sci.* 133 (2007), p. 35-49 (cf. p. 21).

[112] S. SOM et S. AGGARWAL. "Effects of primary breakup modeling on spray and combustion characteristics of compression ignition engines". In: *Combust. Flame* 157 (2010), p. 1179-1193 (cf. p. 21).
[113] D. SUN et W. TAO. "A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows". In: *Int. J. Heat Mass Transfer* 53 (2010), p. 645-655 (cf. p. 30).

[114] M. SUSSMAN et E. G. PUCKETT. "A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows". In: *J. Comput. Phys.* 162 (2000), p. 301-337 (cf. p. 29).

[115] M. SUSSMAN, P. SMEREKA et S. OSHER. "A level set approach for computing solutions to incompressible two-phase flow". In: *J. Comput. Phys.* 114 (1994), p. 146-159 (cf. p. 28).

[117] M. TAVARES. "Simulation et modélisation multi-échelle d'écoulements diphasiques". Thèse de doct. Université Paris-Est, juin 2019 (cf. p. 30, 34, 36).

[118] M. TAVARES, D.-A. KOFFI BI, E. CHENIER et S. VINCENT. "A Two-Dimensional Second Order Conservative Front-Tracking Method with an Original Marker Advection Approach Based on Jump Relations". In: *Communications in Computational Physics* 27.5 (2020), p. 1550-1589 (cf. p. 37, 45).

[119] M. TAVARES, D.-A. KOFFI BI, E. CHÉNIER et S. VINCENT. "A Front-Tracking Method for Multiphase Flows with a Sharp Interface Representation". In: *Turbulence and Interactions*. Sous la dir. de M. DEVILLE, C. CALVIN, V. COUAILLIER, M. DE LA LLAVE PLATA, J.-L. ESTIVALÈZES, T. H. LÊ et S. VINCENT. Cham: Springer International Publishing, 2021, p. 189-195 (cf. p. 37).

[120] D. J. TORRES et J. U. BRACKBILL. "The Point-Set Method: Front-Tracking without Connectivity". In: *J. Comput. Phys.* 165 (2000), p. 620-644 (cf. p. 29, 30, 47).

[122] G. TRYGGVASON, B. BUNNER, A. ESMAEELI, D. JURIC, N. AL-RAWAHI, W. TAUBER, J. HAN, S. NAS et Y.-J. JAN. "A Front-Tracking Method for the Computations of Multiphase Flow". In: *J. Comput. Phys.* 169.2 (2001), p. 708-759 (cf. p. 29-31, 76, 77).

[123] G. TRYGGVASON, R. SCARDOVELLI et S. ZALESKI. *Direct Numerical Simulations of Gas-Liquid Multiphase Flows*. 1re éd. Cambridge University Press, 2011 (cf. p. 31, 77, 82, 85, 90, 93, 95, 107).

[124] S. O. UNVERDI et G. TRYGGVASON. "A front-tracking method for viscous, incompressible, multi-fluid flows". In: *J. Comput. Phys.* 100.1 (1992), p. 25-37 (cf. p. 28).

[126] E. VILLERMAUX. "Fragmentation". In: *Annu. Rev. Fluid Mech.* 39 (2007), p. 419-446 (cf. p. 19).

[127] E. VILLERMAUX et B. BOSSA. "Drop fragmentation on impact". In: *J. Fluid Mech.* 668 (2011), p. 412-435 (cf. p. 21).

[129] S. WILDEMAN, C. W. VISSER, C. SUN et D. LOHSE. "On the spreading of impacting drops". In: *J. Fluid Mech.* 805 (2016), p. 636-655 (cf. p. 21).

[130] F. XIAO, Y. HONMA et T. KONO. "A simple algebraic interface capturing scheme using hyperbolic tangent function". In: *Int. J. Numer. Methods Fluids* 48 (2005), p. 1023-1040 (cf. p. 27).

[131] F. XIAO, S. LI et C. CHEN. "Revisit to the THINC scheme: A simple algebraic VOF algorithm". In: *J. Comput. Phys.* 230 (2011), p. 7086-7092 (cf. p. 27).

[135] D. L. YOUNGS. "Time-dependent multi-material flow with large fluid distortion." In: *Numerical methods for Fluid Dynamics* 24(2) (1982), p. 273-285 (cf. p. 26).

[137] D. ZHAO, Y. TOBA, Y. SUZUKI et S. KOMORI. "Effect of wind waves on air-sea gas exchange: proposal of an overall CO2 transfer velocity formula as a function of breaking-wave parameter". In: *Tellus* (2003), p. 478-487 (cf. p. 19).
# Chapitre 2

# Gestion du maillage de l'interface

## Sommaire

<table><tr><td>2.1</td><td>Notation</td><td>45</td></tr><tr><td>2.2</td><td>Transport des marqueurs de l'interface</td><td>48</td></tr><tr><td>2.3</td><td>Enrichissement</td><td>49</td></tr><tr><td>2.3.1</td><td>Un nouvel élément : ($n_e^{(n)}$ - 3) = 1</td><td>52</td></tr><tr><td>2.3.2</td><td>Deux nouveaux éléments : ($n_e^{(n)}$ - 3) = 2</td><td>53</td></tr><tr><td>2.3.3</td><td>Trois nouveaux éléments : ($n_e^{(n)}$ - 3) = 3</td><td>54</td></tr><tr><td>2.4</td><td>Désenrichissement</td><td>55</td></tr><tr><td>2.5</td><td>Amélioration de la qualité du maillage</td><td>57</td></tr><tr><td>2.5.1</td><td>Échange d'arête</td><td>57</td></tr><tr><td>2.5.2</td><td>Lissage</td><td>60</td></tr><tr><td>2.5.3</td><td>Suppression de mauvais triangles</td><td>61</td></tr><tr><td>2.6</td><td>Conservation du volume</td><td>61</td></tr><tr><td>2.7</td><td>Validation - Interface sphérique dans un champ tourbillonnaire</td><td>63</td></tr></table>

## 2.1 Notation

Dans la méthode front-tracking (suivi d'interface), plusieurs éléments sont à prendre en compte pour effectuer un suivi efficace de l'évolution complexe de l'interface. Des marqueurs de l'interface aux éléments triangulaires discrétisant la surface, en passant par la désignation de leurs voisins, ce sont autant de grandeurs que nous définissons ici de manière formelle pour une identification plus simple par la suite. Nous nous appuyons sur la nomenclature utilisée par TAVARES et al. [118], que nous étendons au front-tracking 3D.

<figure><img src="image_11.png" /><figcaption>FIGURE 2.1 – Le volume $\Omega$ est délimité par une surface $\partial\Omega$. La normale extérieure $\bar{n}_{ext}$ pointe vers le milieu extérieur à $\Omega$.</figcaption></figure>
Ainsi, soit un volume $\Omega$, délimité par la surface $\partial\Omega$. Cette surface est l'interface séparant l'intérieur de $\Omega$ du milieu extérieur. Nous définissons par $\bar{n}_{ext}$ la normale extérieure à $\partial\Omega$ (Fig. 2.1). Nous notons $\mathcal{D} = (\mathcal{E}, \mathcal{M}, \mathcal{T})$ l'ensemble caractérisant la discrétisation de l'interface $\partial\Omega$. Il est constitué de :

— $\mathcal{E}$ l'ensemble des $N_{\mathcal{E}}$ triangles discrétisant la surface. Les **éléments** de $\mathcal{E}$ sont notés $e$ et sont indistinctement appelés triangles ou éléments.

1. Le terme **arête** est employé pour définir géométriquement les éléments $e \in \mathcal{E}$. Ainsi, les éléments $e$ étant des triangles, ils ont 3 arêtes.

2. Le terme **segment** relève de la discrétisation des éléments $e$. En effet, les arêtes de $e$ peuvent être scindés en deux segments par l'ajout d'un point au milieu de l'arête dans l'optique de garder une discrétisation fine de l'interface.

3. Le **nombre de segments** constituant alors le triangle $e$ est noté $n_e$, avec $n_e \in [3, 6]$.

4. Le maillage est dit **non-conforme** quand il existe un élément $e \in \mathcal{E}$ tel que $n_e > 3$.

— $\mathcal{M}$ l'ensemble des $N_{\mathcal{M}}$ marqueurs $m \in \mathcal{M}$ du maillage de l'interface.

1. Les **marqueurs** constituent les extrémités des segments composant les triangles $e \in \mathcal{E}$. Par ailleurs, la quantité $n_e$ désigne aussi le nombre de marqueurs de l'élément $e$.

2. Le **i-ème marqueur** $m$ de l'élément $e$ est identifié par $m_{[e,i]}$, avec $i \in [1, n_e]$.

3. Les **coordonnées** associées à $m_{[e,i]}$ sont désignées par $\bar{x}_{e,i}$ et sa **vitesse** par $\bar{v}_{e,i}$.

4. Les **sommets du triangle** $e$ sont désignés par $m_{[e,i]}$, avec $i \le 3$. Les marqueurs ajoutés au **milieu des arêtes** du triangle $e$ sont indiqués par $m_{[e,i]}$, avec $i > 3$.

— $\mathcal{T}$ l'ensemble des $N_{\mathcal{T}}$ $n_e$-uplets tels que $\forall e \in \mathcal{E}$, $(m_{[e,1]}, \dots, m_{[e,n_e]}) \in \mathcal{T}$. Les $n_e$-uplets étant ordonnés, **l'orientation** ou sens de parcours local des marqueurs est défini de telle sorte que :

$$ \forall e \in \mathcal{E}, \forall i \in [2, n_e - 1], \left( (\bar{x}_{e,i+1} - \bar{x}_{e,i}) \wedge (\bar{x}_{e,i-1} - \bar{x}_{e,i}) \right) \cdot \bar{n}_{ext} > 0 \quad (2.1) $$

Ainsi, cet ordre se traduit par le fait que toute normale à l'élément $e$, construite à partir de trois marqueurs d'indices consécutifs, hormis les premiers et derniers marqueurs, possède le même sens que la normale extérieure. Cela nous permet d'avoir une **numéro-tation locale** des marqueurs. Cette numérotation est temporaire et adaptée au nombre de marqueurs de l'élément $e$. Nous nous ramènerons aux numérotations présentées sur la figure 2.2, selon les valeurs de $n_e$.

L'ensemble de discrétisation $\mathcal{D}$ ainsi défini, nous introduisons d'autres notions qui nous sont utiles dans la gestion du maillage de l'interface.

— **marqueurs adjacents** : la numérotation que nous avons établie nous permet de trouver, à partir d'un marqueur $i$ de l'élément $e$, le marqueur qui le précède ou qui le suit immédiatement dans le triangle. Ainsi, le **marqueur précédant** immédiatement $m_{[e,i]}$ dans l'élément $e$ a pour indice $i_0$. Le **marqueur suivant** est repéré par l'indice $i_2$. Nous définissons l'opérateur $\text{mod}(a, b)$ comme le reste de la division euclidienne de $a$ par $b$. En d'autres termes :

$$ \text{mod}(a, b) = r, \text{ avec } a = bq + r, \quad (q, r) \in \mathbb{N}^2, \ r < b $$

Ainsi, les indices $i_0$ et $i_2$ sont obtenus de la manière suivante :

$$ i_0 = \begin{cases} i + 3 & , \text{ si } (i + 3) \le n_e \\ \text{mod}(i + 1, 3) + 1 & , \text{ sinon.} \end{cases} \quad (2.2) $$

$$ i_2 = \begin{cases} \text{mod}(i, 3) + 1 & , \text{ si } n_e = 3 \\ \text{mod}(i, 4) + 1 + 3 \times \text{E}[(n_e - i - 1)/3] \times \text{E}[(6 - i)/3] & , \text{ sinon.} \end{cases} \quad (2.3) $$

où $\text{E}[\bullet]$ désigne la partie entière de $\bullet$.
<figure><img src="image_12.png" /><figcaption>FIGURE 2.2 – Numérotation des marqueurs de l'élément e selon les valeurs de $n_e$.</figcaption></figure>

- **identification des segments** : l'orientation ou sens de parcours imposé par $\bar{n}_{ext}$ induit une orientation aux segments constituant les triangles. Nous pouvons donc identifier un marqueur initial et un marqueur final pour chaque segment de $e$. Par conséquent, nous appelons le segment $i$ de l'élément $e$, le segment formé du marqueur final $m_{[e,i]}$ et du marqueur initial $m_{[e,i_0]}$, avec $i_0$ défini par l'équation (2.2). Par ailleurs, la **longueur du segment** $i$ de $e$ est notée $d_{e,i}$. Elle est définie par l'équation (2.4) :

$$d_{e,i} = \|\bar{x}_{e,i} - \bar{x}_{e,i_0}\| \qquad (2.4)$$

- **connectivité du maillage** : elle correspond à la connaissance des triangles et de leurs voisins, à savoir les éléments ayant un segment en commun. Il existe des méthodes front-tracking sans connectivité [120]. Nous avons choisi de garder une connectivité, car elle facilite l'accès à des propriétés géométriques de l'interface telles que la normale et la courbure, utiles pour le calcul des tensions de surfaces. Le voisin de $e \in \mathcal{E}$ par son segment $i$ est noté $e_{[i]} \in \mathcal{E}$. Ainsi, dire que les éléments $e, e' \in \mathcal{E}$ sont voisins se traduit par l'équation (2.5) :

$$\exists i \in [1, n_e], \exists k \in [1, n_{e'}] \mid e_{[i]} = e' \text{ et } e'_{[k]} = e \qquad (2.5)$$

L'ensemble de ces notations est synthétisé sur la figure 2.3. Elle illustre une situation où l'élément initial $e$ possède deux arêtes trop grandes. Elles sont donc coupées en ajoutant deux marqueurs supplémentaires, faisant ainsi de $e$ un triangle composé de 5 segments ($n_e = 5$). Cette configuration non conforme est transitoire mais nécessaire à la gestion et l'évolution du maillage de l'interface.

- **voisinage d'un marqueur** : c'est l'ensemble des couronnes successives d'éléments autour du sommet considéré. Ainsi, le $k$-ème voisinage du marqueur $m \in \mathcal{M}$ est noté
<figure><img src="image_13.png" /><figcaption>FIGURE 2.3 – L’élément triangulaire initial $e$, constitué des marqueurs $m_{[e,i]}$, $i \in [1,3]$, a ses arêtes 1 et 2 trop grandes, ce qui entraîne l’ajout des marqueurs supplémentaires $m_{[e,4]}$ et $m_{[e,5]}$ pour les scinder. L’orientation ou ordre de la numérotation est donnée par la normale sortante $\bar{n}_{ext}$.</figcaption></figure>

$\mathcal{N}_{[k,m]} \subset \mathcal{E}$. De manière précise, nous avons :

$$ 
\begin{aligned}
\mathcal{N}_{[1,m]} &= \{e \in \mathcal{E} \mid \exists! j \in [1, n_e], m_{[e,j]} = m\} \\
\mathcal{N}_{[k,m]} &= \{e \in \mathcal{E} \setminus \mathcal{N}_{[k-1,m]} \mid \exists e' \in \mathcal{N}_{[k-1,m]}, e \cap e' \neq \emptyset\} \cup \mathcal{N}_{[k-1,m]}, \quad k \geq 2
\end{aligned} 
\quad (2.6) 
$$

Nous voyons, d'après l'équation (2.6), que le premier voisinage $\mathcal{N}_{[1,m]}$ du marqueur $m$ est l'ensemble des triangles ayant $m$ en commun. Son second voisinage $\mathcal{N}_{[2,m]}$ comprend, en plus des triangles du premier voisinage, les éléments qui ont un segment ou un sommet en commun avec les triangles du premier voisinage. Et ainsi de suite. La **connectivité du marqueur** $m$ correspond au cardinal de son premier voisinage $\mathcal{N}_{[1,m]}$. — **évolution temporelle** : l'ensemble de discrétisation $\mathcal{D} = (\mathcal{E}, \mathcal{M}, \mathcal{T})$ évolue au cours du temps. Afin de prendre en compte le caractère instationnaire, cet ensemble et ses composantes sont munis d'un exposant $n$ pour définir son état à l'instant $t^{(n)}$. Ainsi, l'ensemble caractérisant la discrétisation de l'interface à l'instant $t^{(n)}$ est noté $\mathcal{D}^{(n)} = (\mathcal{E}^{(n)}, \mathcal{M}^{(n)}, \mathcal{T}^{(n)})$.

## 2.2 Transport des marqueurs de l'interface

Dans la méthode front-tracking, l'interface est suivie de manière explicite. Pour ce faire, en plus du maillage fixe utilisé pour la résolution des équations de Navier-Stokes, l'interface bénéficie d'un maillage qui lui est propre. Ce maillage est constitué de triangles et ce sont les sommets de ces triangles qui sont transportés pour rendre compte de l'évolution de l'interface. La figure 2.4 nous montre un exemple :

Le transport des sommets des triangles, ou encore marqueurs, se fait de manière lagrangienne. Ainsi, pour tout marqueur $i$ de $e$, $m_{[e,i]}$, l'évolution de sa position au cours du temps est définie par :

$$ \frac{d}{dt}(\bar{x}_{e,i}) = \bar{v}_{e,i} \quad (2.7) $$
<figure><img src="image_14.png" /><figcaption>FIGURE 2.4 – La surface de la sphère est maillée par des triangles. Les sommets sont transportés et suivent les changements topologiques de l'interface</figcaption></figure>

Afin d'avoir la position des marqueurs, nous utilisons un schéma d'interpolation de Runge-Kutta d'ordre 2 (RK2). Ainsi, en notant $\Delta t$ l'écart temporel entre les instants $t^n$ et $t^{n+1}$, les positions des marqueurs à l'instant $t^{n+1}$ se déterminent de la manière suivante :

$$ \begin{aligned} \bar{x}_{e,i}^{(*)} &= \bar{x}_{e,i}^{(n)} + \Delta t \cdot \bar{v}_{e,i}^{(n)} \\ \bar{x}_{e,i}^{(n+1)} &= \bar{x}_{e,i}^{(n)}/2 + \left( \bar{x}_{e,i}^{(*)} + \Delta t \cdot \bar{v}_{e,i}^{(n+1)} \right)/2 \end{aligned} \quad (2.8) $$

$\bar{v}_{e,i}^{(n+1)}$ est obtenue par une méthode d'extrapolation temporelle d'ordre 2 d'Adams-Bashforth :

$$ \bar{v}_{e,i}^{(n+1)} = 2\bar{v}_{e,i}^{(n)} - \bar{v}_{e,i}^{(n-1)} \quad (2.9) $$

La vitesse des marqueurs à un instant donné $k$ est calculée avec la méthode PERM (Parabolic Edge Reconstruction Method), comme définie dans les travaux de McDERMOTT et al. [74]. Elle est déterminée à partir des vitesses adjacentes du maillage eulérien :

<figure><img src="image_15.png" /><figcaption>FIGURE 2.5 – Détermination de la vitesse des marqueurs par la méthode PERM. Elle est obtenue à partir des vitesses adjacentes issues du maillage eulérien.</figcaption></figure>

## 2.3 Enrichissement

Afin de suivre la dynamique de l'interface, il est nécessaire d'avoir suffisamment de triangles pour décrire le mouvement et les déformations de celle-ci. Ainsi, nous définissons l'enrichissement comme étant l'ajout d'éléments et de marqueurs sur le maillage lagrangien, dans l'optique de capter au mieux les changements topologiques de l'interface. En effet, avec une discrétisation initiale $D^{(0)}$, si $D^{(n)}$ n'évolue pas au cours du temps, nous pouvons arriver à des situations où
le maillage lagrangien ne décrit plus l'interface et n'arrive plus à capter toutes les déformations topologiques de l'interface. Un cas typique est illustré sur la figure 2.6. Nous voyons la déformation maximale d'une sphère de rayon 0.15, initialement placée dans un cube unitaire à la position (0.35, 0.35, 0.35). Elle est soumise à un champ incompressible de vitesse, combinant des déformations dans les planes $(x - y)$ et $(x - z)$. Ce champ de vitesse, initialement proposé par LEVEQUE [65], et repris par ENRIGHT et al. [30], est construit selon l'équation (2.10) :

$$ 
\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{2.10} \\ 
w(x, y, z) &= -\sin(2\pi x)\sin(2\pi y)\sin^2(\pi z) \nonumber 
\end{align} 
$$

Il est modulé en temps par une fonction sinusoïdale $2T$-périodique $\cos(\pi t/T)$. Ainsi, la défor-

<figure><img src="image_16.png" /><figcaption>FIGURE 2.6 – Sphère au maximum de déformation ($t = T/2$) dans un champ incompressible de vitesse (Eq. (2.10)). Les deux sphères ont le même maillage initial. A gauche (a), nous n'appliquons pas d'enrichissement. A droite (b), l'enrichissement est activé.</figcaption></figure>

mation maximale est obtenue pour $t = T/2$. Comme nous pouvons le voir sur la figure 2.6a, il est clair que l'absence d'enrichissement peut conduire à des aberrations, telles que des triangles qui se recouvrent, s'entrecoupent, ou encore la présence de triangles dans des zones où il ne devrait pas en avoir. Ce sont des phénomènes que nous évitons en activant l'enrichissement (Fig. 2.6b).

La mise en œuvre de l'enrichissement passe par la détermination de critères pour raffiner ou subdiviser le maillage de l'interface. Différentes approches existent telles que celles proposées par ROGHAIR et al. [99]. Par la suite, nous présentons les critères que nous avons pour effectuer l'enrichissement du maillage de l'interface.

**Longueur des arêtes** : les arêtes de l'élément $e$ sont scindés en deux si leur longueur $d_{e,i}$ est supérieure à la taille $h$ des cellules de la grille eulérienne.

$$ \text{si } d_{e,i} \geq h \text{ alors enrichir} \tag{2.11} $$

L'enrichissement est effectué si l'équation (2.11) est satisfaite. Il est constitué de deux étapes, à savoir, d'une part **la scission des arêtes** et donc l'ajout de nouveaux marqueurs sur le triangle $e$, et d'autre part la **création de nouveaux triangles** à partir de $e$. Cet ensemble de processus conduit à un maillage interfacial transitoirement non-conforme (voir section 2.1, page 46) comme illustré sur la figure 2.3. Cela se traduit par :

$$ \exists e \in \mathcal{E}^{(n)} \mid n_e > 3 $$

A la fin du processus d'enrichissement, le maillage redevient conforme. En d'autres termes :

$$ \forall e \in \mathcal{E}^{(n)}, \quad n_e = 3 $$

**Scission des arêtes** : tous les éléments $e$ de $\mathcal{E}^{(n)}$, hormis les nouveaux formés au cours de l'enrichissement, sont analysés tour à tour. En notant $\mathcal{E}_{new}^{(n)}$ l'ensemble de éléments nouvellement
formés lors de l'enrichissement, nous pouvons synthétiser la phase de scission des arêtes comme suit :

$$ \forall e \in \mathcal{E}^{(n)} \setminus \mathcal{E}_{\text{new}}^{(n)}, \text{ pour } A = 1 + (n_e^{(n)} - 3) \text{ fixé, } \forall i \in [A, 3] $$

$$ \text{si } d_{e,i}^{(n)} \geq h \text{ alors } \begin{cases} -n_e^{(n)} = n_e^{(n)} + 1 \\ - \text{Création du marqueur } m_{[e,n_e^{(n)}]} \text{ tel que } \bar{x}_{e,n_e^{(n)}} = (\bar{x}_{e,i} + \bar{x}_{e,i_0})/2 \quad (2.12) \\ - N_M^{(n)} = N_M^{(n)} + 1 \end{cases} $$

En d'autres termes, nous regardons d'abord le nombre $n_e^{(n)}$ de segments (ou de marqueurs) qui constitue l'élément $e$. Si $n_e^{(n)} = 3$, alors $A = 1$ et nous vérifions donc la longueur de ses trois arêtes ($i \in [1, 3]$). Dans le cas $n_e^{(n)} \in [4, 5]$, une ou deux arêtes ont déjà été scindées par l'effet de ses voisins. Ainsi, pour $n_e^{(n)} = 5$, $A = 3$ et $i \in [3, 3]$. Seule l'arête 3 est analysée. Dans l'éventualité où $n_e^{(n)} = 6$, aucune analyse supplémentaire des arêtes n'est à effectuer. Ainsi, quel que soit le cas de figure, si les conditions sont remplies, de nouveaux marqueurs sont rajoutés sur les arêtes de $e$, et son nombre de segments $n_e^{(n)}$ est mis à jour, de même que le nombre total $N_M^{(n)}$ de marqueurs du maillage de l'interface.

Pour la création des nouveaux triangles, trois cas de figures sont possibles, selon les valeurs de $n_e^{(n)}$. En effet, $(n_e^{(n)} - 3)$ représente le nombre de nouveaux éléments à créer.
2.3.1 Un nouvel élément : $(n_e^{(n)} - 3) = 1$

L'élément *e* possède une seule grande arête qui a été scindée. Il est ainsi formé de 4 segments (Fig. 2.7a). Ses voisins $e_{[1]}$ et $e_{[4]}$ peuvent être différents ou non. En effet, les pointillés sur la figure 2.7 représente le fait que le voisin de *e* a été coupé ou non. S'il a été coupé, $e_{[1]} \neq e_{[4]}$. Dans le cas contraire, son voisin est l'unique élément $e_{[1]} = e_{[4]}$. Deux nouveau éléments $e'$ et $e''$ sont créés en construisant une arête entre le nouveau marqueur $m_{[e,4]}$ et le marqueur opposé $m_{[e',2]}$. Les connectivités sont ensuite mises à jour, tandis que le triangle *e* est supprimé (Fig. 2.7b). Ainsi, un seul élément est effectivement ajouté au maillage. L'ensemble des opérations à mener lors de cet enrichissement est synthétisé dans l'algorithme 2.1.

<figure><img src="image_17.png" /><figcaption>FIGURE 2.7 – L'élément *e* est constitué de 4 segments (a). Deux nouveaux éléments $e'$ et $e''$ sont formés et l'élément *e* est supprimé. Les connectivités sont mises à jour (b).</figcaption></figure>

$m_{[e',1]} = m_{[e,4]}$ ; $m_{[e',2]} = m_{[e,2]}$ ; $m_{[e',3]} = m_{[e,3]}$ ! Création de l'élément $e'$

1

$m_{[e'',1]} = m_{[e,1]}$ ; $m_{[e'',2]} = m_{[e,2]}$ ; $m_{[e'',3]} = m_{[e,4]}$ ! Création de l'élément $e''$

2

$e'_{[1]} = e_{[4]}$ ; $e'_{[2]} = e''$ ; $e'_{[3]} = e_{[3]}$ ! Connectivité de $e'$

3

$e''_{[1]} = e_{[1]}$ ; $e''_{[2]} = e_{[2]}$ ; $e''_{[3]} = e'$ ! Connectivité de $e''$

4

! Mise à jour de la connectivité des voisins

5

do $i = 1, n_e$

6

$\tilde{e} = e_{[i]}$

7

$k = 1$

8

do while $(\tilde{e}_{[k]} \neq e)$ ! Recherche du segment commun entre $\tilde{e}$ et $e$

9

$k = k + 1$

10

end do

11

if ($i \leq 2$) then

12

$\tilde{e}_{[k]} = e''$

13

else

14

$\tilde{e}_{[k]} = e'$

15

end if

16

end do

17

Suppression de *e*

18

$N_{\mathcal{E}}^{(n)} = N_{\mathcal{E}}^{(n)} + 1$ ; $N_{\mathcal{T}}^{(n)} = N_{\mathcal{T}}^{(n)} + 1$ ! Mise à jour des cardinaux de $\mathcal{E}^{(n)}$ et $\mathcal{T}^{(n)}$

19

<caption>Listing 2.1 – Enrichissement - Ajout effectif d'un élément</caption>
### 2.3.2 Deux nouveaux éléments : ($n_e^{(n)} - 3$) = 2

Dans le cas où le triangle $e$ est discrétisé en 5 segments, deux segments sont construits entre $m_{[e,4]}$ et $m_{[e,5]}$ d'une part, et entre $m_{[e,5]}$ et $m_{[e,3]}$ d'autre part. Deux éléments sont effectivement ajoutés au maillage lagrangien (Fig. 2.8). Nous résumons les diverses étapes qui conduisent à ce résultat dans l'algorithme 2.2.

<figure><img src="image_18.png" /><figcaption>FIGURE 2.8 – L'élément $e$ est formé de 5 segments (a). Trois nouveaux éléments $e'$, $e''$ et $e'''$ sont créés et l'élément $e$ est supprimé. Le choix de construire les trois éléments et de les disposer dans cet ordre est fixé. Les voisins sont également mis à jour (b).</figcaption></figure>

1. $m_{[e',1]} = m_{[e,4]}$ ; $m_{[e',2]} = m_{[e,5]}$ ; $m_{[e',3]} = m_{[e,3]}$ ! Création de l'élément $e'$

2. $m_{[e'',1]} = m_{[e,1]}$ ; $m_{[e'',2]} = m_{[e,5]}$ ; $m_{[e'',3]} = m_{[e,4]}$ ! Création de l'élément $e''$

3. $m_{[e'''},1] = m_{[e,2]}$ ; $m_{[e'''},2] = m_{[e,3]}$ ; $m_{[e'''},3] = m_{[e,5]}$ ! Création de l'élément $e'''$

4. $e'_{[1]} = e_{[4]}$ ; $e'_{[2]} = e''$ ; $e'_{[3]} = e'''$ ! Connectivité de $e'$

5. $e''_{[1]} = e_{[1]}$ ; $e''_{[2]} = e_{[5]}$ ; $e''_{[3]} = e'$ ! Connectivité de $e''$

6. $e'''_{[1]} = e_{[2]}$ ; $e'''_{[2]} = e_{[3]}$ ; $e'''_{[3]} = e'$ ! Connectivité de $e'''$

7. ! Mise à jour de la connectivité des voisins

8. do $i = 1, n_e$

9. $\tilde{e} = e_{[i]}$

10. $k = 1$

11. do while $(\tilde{e}_{[k]} \neq e)$ ! Recherche du segment commun entre $\tilde{e}$ et $e$

12. $k = k + 1$

13. end do

14. if $(\text{mod}(i, 4) = 1) \quad \tilde{e}_{[k]} = e''$ ! $i \in \{1, 5\}$

15. if $(2 \le i \le 3) \quad \tilde{e}_{[k]} = e'''$

16. if $(i = 4) \quad \tilde{e}_{[k]} = e'$

17. end do

18. Suppression de $e$

19. $N_{\mathcal{E}}^{(n)} = N_{\mathcal{E}}^{(n)} + 2$ ; $N_{\mathcal{T}}^{(n)} = N_{\mathcal{T}}^{(n)} + 2$ ! Mise à jour des cardinaux de $\mathcal{E}^{(n)}$ et $\mathcal{T}^{(n)}$

<caption>Listing 2.2 – Enrichissement - Ajout effectif de deux éléments</caption>
### 2.3.3 Trois nouveaux éléments : $(n_e^{(n)} - 3) = 3$

Le choix de construction des éléments est beaucoup moins compliqué que précédemment. L'élément $e$ est divisé en 4 petits éléments en construisant des segments entre les marqueurs $m_{[e,i]}$, $i \in [4, 6]$ (Fig. 2.9). La procédure est résumée dans l'algorithme 2.3.

<figure><img src="image_19.png" /><figcaption>FIGURE 2.9 – L'élément $e$ possède 6 segments (a). Quatre nouveaux éléments $e'$, $e''$, $e'''$ et $e''''$ sont créés tandis que l'élément $e$ est supprimé. Les voisins sont également mis à jour (b).</figcaption></figure>

```
1 m[e',1] = m[e,5] ; m[e',2] = m[e,6] ; m[e',3] = m[e,4] ! Création de l'élément e'
2 m[e'',1] = m[e,1] ; m[e'',2] = m[e,5] ; m[e'',3] = m[e,4] ! Création de l'élément e''
3 m[e'''',1] = m[e,2] ; m[e'''',2] = m[e,6] ; m[e'''',3] = m[e,5] ! Création de l'élément e''''
4 m[e''''',1] = m[e,3] ; m[e''''',2] = m[e,4] ; m[e''''',3] = m[e,6] ! Création de l'élément e'''''
5 e'[1] = e'' ; e'[2] = e''' ; e'[3] = e'''' ! Connectivité de e'
6 e''[1] = e[1] ; e''[2] = e[5] ; e''[3] = e' ! Connectivité de e''
7 e'''[1] = e[2] ; e'''[2] = e[6] ; e'''[3] = e' ! Connectivité de e'''
8 e''''[1] = e[3] ; e''''[2] = e[4] ; e''''[3] = e' ! Connectivité de e''''
9 ! Mise à jour de la connectivité des voisins
do i = 1, n_e
10
11   ã = e[i]
12   k = 1
13   do while (ã[k] ≠ e) ! Recherche du segment commun entre ã et e
14       k = k + 1
15   end do
16   if (mod(i, 4) = 1) then ! i ∈ {1, 5}
17       ã[k] = e''
18   else if (mod(i, 4) = 2) then ! i ∈ {2, 6}
19       ã[k] = e'''
20   else
21       ã[k] = e''''
22   end if
23 end do
24 Suppression de e
25 N_ε^{(n)} = N_ε^{(n)} + 3 ; N_T^{(n)} = N_T^{(n)} + 3 ! Mise à jour des cardinaux de E^{(n)} et T^{(n)}
```

<caption>Listing 2.3 – Enrichissement - Ajout effectif de trois éléments</caption>

Ainsi, l'ensemble de ces procédures (Eq. (2.12) et Algos. 2.1-2.3) nous permet d'avoir suffisamment d'éléments pour décrire des interfaces qui s'étirent. Cependant, les écoulements di-
phasiques peuvent présenter également des zones où l'interface se contracte, se traduisant par une diminution des longueurs des arêtes du maillage interfacial. Dans ce genre de situation, il convient de supprimer des marqueurs afin de garder une stabilité du maillage.

## 2.4 Désenrichissement

L'enrichissement densifie le maillage de l'interface. Toutefois, il n'est pas toujours nécessaire d'avoir beaucoup de triangles pour décrire la topologie de l'interface. En effet, pour des zones de forte courbure, avoir un maillage assez raffiné permet d'avoir une bonne évaluation de l'interface. Cependant, dans des zones de faible courbure, disposer d'un maillage trop fin est prohibitif. De surcroît, des instabilités non-physiques peuvent apparaître sur un maillage trop raffiné [33].

Par conséquent, nous désenrichissons plus dans les zones de faible courbure que dans celles de forte courbure. Nous présentons nos critères pour le désenrichissement, dont les deux premiers sont basés sur les travaux de ROGHAIR et al. [99].

**Longueur des arêtes** : les arêtes dont la longueur est plus petite que $h/4$ (le quart de la taille des cellules eulériennes) sont marquées pour une potentielle suppression.

$$d_{e,i}^{(n)} \leq h/4 \qquad (2.13)$$

**Rugosité de l'arête** : elle traduit l'appartenance d'une arête à une zone courbée ou non. La rugosité d'une arête est construite à partir de la rugosité des marqueurs qui le délimitent. En effet, $\forall e \in \mathcal{E}^{(n)}$, soit le segment $i$ de $e$, $i \in [1, 3]$, compris entre les marqueurs $m_{[e,i]}$ et $m_{[e,i_0]}$, avec $i_0$ défini par l'équation (2.2) (voir section 2.1, page 46). Soient les premiers voisinages associés à ces deux marqueurs, $\mathcal{N}_{[1,m_{[e,i]})}$ et $\mathcal{N}_{[1,m_{[e,i_0]}]}$ (voir section 2.1, page 47). La rugosité $r_{[e,i]}$ de l'arête $i$ de l'élément $e$ est définie comme la demi-somme de la rugosité de chaque marqueur de l'arête $i$, $r(m_{[e,i]})$ et $r(m_{[e,i_0]})$.

$$r(m_{[e,i]}) = \min_{(e,e')} [(1 + \bar{n}_e \cdot \bar{n}_{e'})/2], \quad e, e' \in \mathcal{N}_{[1,m_{[e,i]})} \\ r(m_{[e,i_0]}) = \min_{(e,e')} [(1 + \bar{n}_e \cdot \bar{n}_{e'})/2], \quad e, e' \in \mathcal{N}_{[1,m_{[e,i_0]}]}$$

Et la rugosité de l'arête $i$ est :

$$r_{[e,i]} = [r(m_{[e,i]}) + r(m_{[e,i_0]})] / 2 \qquad (2.14)$$

où $\bar{n}_e$ et $\bar{n}_{e'}$ sont les normales unitaires aux éléments $e$ et $e'$. De plus, $e$ et $e'$ sont voisins (Eq. (2.5), section 2.1, page 47). Ainsi, dans les zones planes, $\bar{n}_e \cdot \bar{n}_{e'} = 1$ et $r_{[e,i]} = 1$. Au contraire, dans les zones de forte courbure, $\bar{n}_e \cdot \bar{n}_{e'} = -1$ et $r_{[e,i]} = 0$. Notre objectif est d'avoir un nombre adéquat d'éléments pour avoir une bonne description de l'interface. Ainsi, le segment $i$ de $e$ est supprimable si :

$$r_{[e,i]} > 0.99 \qquad (2.15)$$

**Connectivité des marqueurs** : nous imposons au maillage la contrainte que chaque marqueur $m \in \mathcal{M}$ ait une connectivité (voir section 2.1, page 48) minimale de 4. Avoir une connectivité plus basse (égale à 3) entraînerait des erreurs lors de l'évaluation des propriétés géométriques de l'interface (voir Annexe B). Il est donc nécessaire que le désenrichissement ne franchise pas cette limite imposée de connectivité des marqueurs. En effet, le désenrichissement supprime une arête et, par voie de conséquence, les éléments qui la partagent. Les marqueurs opposés à cette arête ont leur connectivité baissée de 1. Ainsi, si l'arête $i$ est à supprimer, il faut :

$$\left\{ \begin{array}{l} \text{card} \left( \mathcal{N}_{[1,m_{[e,i]})} \right) + \text{card} \left( \mathcal{N}_{[1,m_{[e,i_0]}]} \right) - 4 \geq 4 \\ \text{card} \left( \mathcal{N}_{[1,m_{[e,i_2]}]} \right) > 4 \text{ avec } i \text{ tel que } e_{[i]} = e' \\ \text{card} \left( \mathcal{N}_{[1,m_{[e',i_2]}]} \right) > 4 \text{ avec } i \text{ tel que } e'_{[i]} = e \end{array} \right. \qquad (2.16)$$
En résumé, la somme des connectivités des marqueurs de l'arête à supprimer est supérieure ou égale à 8, et la connectivité de chaque marqueur opposé à cette arête est supérieure stricte à 4.

**Éviter la formation d'arêtes trop longues** : le désenrichissement ne doit pas générer de travail supplémentaire à la procédure d'enrichissement. Il ne faudrait pas qu'il se forme d'arêtes ayant une longueur supérieur à la limite autorisée $h$ lors du désenrichissement. Nous indiquons par (*) l'état hypothétique où le désenrichissement est effectué sur les critères précédents (Eqs. (2.13), (2.15), (2.16)). La suppression de l'arête $i$ de $e$ entraîne l'effondrement de ses extrémités sur son milieu $\bar{x} = (\bar{x}_{e,i} + \bar{x}_{e,i_0})/2$. En d'autres termes, $\bar{x}_{e,i} = \bar{x}$ et $\bar{x}_{e,i_0} = \bar{x}$. Connaissant le point d'effondrement de l'arête $i$, les nouvelles distances $d_{e,i}^{(*)}$ peuvent être évaluées sans pour autant modifier le maillage de l'interface. Ainsi, pour éviter la formation d'arêtes trop longues, il faut :

$$ \forall e \in (\mathcal{N}_{[1,m_{[e,i]}]} \cup \mathcal{N}_{[1,m_{[e,i_0]}]}), \forall i \in [1, n_e], d_{e,i}^{(*)} < h \quad (2.17) $$

**Éviter le renversement des triangles** : certaines configurations peuvent entraîner un retournement des triangles lors de la suppression d'arêtes, comme le montre la figure 2.10. Les éléments résultants ont leur normale pointée vers l'intérieur plutôt que vers l'extérieur. Aussi, pour empêcher de telles situations, les normales des éléments entre l'instant $t^{(n)}$ et l'état hypothétique (*) de désenrichissement sont comparées. Ces deux normales doivent pointer vers la même région :

$$ \forall e \in (\mathcal{N}_{[1,m_{[e,i]}]} \cup \mathcal{N}_{[1,m_{[e,i_0]}]}) \setminus (\mathcal{N}_{[1,m_{[e,i]}]} \cap \mathcal{N}_{[1,m_{[e,i_0]}]}), \quad \bar{n}_e^{(n)} \cdot \bar{n}_e^{(*)} \ge 0 \quad (2.18) $$

<figure><img src="image_20.png" /><figcaption>FIGURE 2.10 – L'arête 1 (rouge) de l'élément $e$ remplit les conditions pour être supprimer (a). Les marqueurs des voisins ne changent pas. Seules les coordonnées de $m_{[e',2]}$ et $m_{[e'',3]}$ sont changées. Ils ont pour coordonnées le milieu de l'arête 1 de l'ancien élément $e$. Dans cette configuration, l'élément $e'$ se retrouve renversé, situation qu'il faut éviter (b).</figcaption></figure>

Ainsi, lorsque tous ces critères sont remplis (Eqs. (2.13), (2.15), (2.16), (2.17) et (2.18)), l'arête considérée est supprimée comme illustré sur la figure 2.11.
<figure><img src="image_21.png" /><figcaption>FIGURE 2.11 – L’arête 1 de l’élément *e* est courte et remplit les conditions pour être supprimer (a). Après la procédure, *e* et *e'* sont effacés. Les voisins sont également mis à jour (b).</figcaption></figure>

Les opérations synthétisées dans l'algorithme 2.4 sont alors effectuées.

1 $\bar{x} = (\bar{x}_{e,i} + \bar{x}_{e,i_0})/2$ ! Coordonnées du milieu de l'arête i de e

2 $\bar{x}_{e,i} = \bar{x}$ ! Déplacement de $m_{[e,i]}$ au milieu de l'arête

3 ! Fusion de $m_{[e,i_0]}$ et $m_{[e,i]}$

4 do $\forall \tilde{e} \in \mathcal{N}_{[1,m_{[e,i_0]})}$

5 $k = 1$

6 do while $(m_{[\tilde{e},k]} \neq m_{[e,i_0]})$ ! Position de $m_{[e,i_0]}$ dans $\tilde{e}$

7 $k = k + 1$

8 end do

9 $m_{[\tilde{e},k]} = m_{[e,i]}$

10 end do

11 ! Mise à jour de $M^{(n)}, \mathcal{E}^{(n)}$ et $\mathcal{T}^{(n)}$ et leurs cardinaux

12 Suppression de $m_{[e,i_0]}$

13 Suppression de $e$ et $e'$

14 $N_M^{(n)} = N_M^{(n)} - 1$ ; $N_{\mathcal{E}}^{(n)} = N_{\mathcal{E}}^{(n)} - 2$ ; $N_{\mathcal{T}}^{(n)} = N_{\mathcal{T}}^{(n)} - 2$

15 ! Mise à jour de la connectivité des voisins de $e$ et $e'$

16 $\tilde{e} = e_{[i_2]}$ ; $\exists k \mid \tilde{e}_{[k]} = e$ ; $\tilde{e}_{[k]} = e_{[i_0]}$

17 $\tilde{e} = e_{[i_0]}$ ; $\exists k \mid \tilde{e}_{[k]} = e$ ; $\tilde{e}_{[k]} = e_{[i_2]}$

18 $\tilde{e} = e'_{[i_2]}$ ; $\exists k \mid \tilde{e}_{[k]} = e'$ ; $\tilde{e}_{[k]} = e'_{[i_0]}$

19 $\tilde{e} = e'_{[i_0]}$ ; $\exists k \mid \tilde{e}_{[k]} = e'$ ; $\tilde{e}_{[k]} = e'_{[i_2]}$

<caption>Listing 2.4 – Désenrichissement</caption>

## 2.5 Amélioration de la qualité du maillage

Les deux précédentes procédures visent à conserver un nombre adéquat d'éléments pour avoir une bonne description et un bon suivi des changements topologiques de l'interface. Les procédures qui suivent ont pour but d'améliorer la forme des triangles du maillage lagrangien. En effet, l'idée est d'avoir des triangles les plus équilatéraux possibles.

### 2.5.1 Échange d'arête

L'échange d'arête contribue à cet objectif en homogénéisant la connectivité des marqueurs sur l'ensemble du maillage et en améliorant la qualité des triangles, que nous définirons plus tard.
Le processus est schématisé sur la figure 2.12. Quatre critères sont nécessaires pour échanger les
arêtes.

<figure><img src="image_22.png" /><figcaption>FIGURE 2.12 – L’arête 1 (rouge) de l’élément *e* est échangée pour rendre les éléments *e* et *e'* les plus équilatéraux possibles.</figcaption></figure>

**Inclinaison des triangles** : l'échange d'arête modifie le volume Ω d'une quantité égale au volume du tétraèdre formé par les marqueurs de *e* et *e'* (Fig. 2.13). Or cette modification artificielle doit être limitée car elle va à l'encontre du principe de conservation de la masse. Pour ce faire, nous regardons les normales unitaires $\bar{n}_e$ et $\bar{n}_{e'}$ des éléments voisins *e* et *e'*. En effet, la variation de volume est liée à l'inclinaison entre les deux éléments. Ainsi, nous effectuons l'échange d'arêtes que dans des cas de faible inclinaison (2.19).

$$
\bar{n}_e \cdot \bar{n}_e' \geq \cos(\pi/18) \qquad (2.19)
$$

<figure><img src="image_23.png" /><figcaption>FIGURE 2.13 – Variation de volume engendrée par l’échange d’arête entre les éléments *e* et *e'*.</figcaption></figure>

**Connectivité des marqueurs** : soit *i*, l'arête à échanger, délimitée par les marqueurs $m_{[e,i]}$ et $m_{[e,i_0]}$. L'échange d'arête baisse la connectivité de $m_{[e,i]}$ et $m_{[e,i_0]}$ de 1, et augmente celle de $m_{[e,i_2]}$ de 1 (Fig. 2.12). Or pour le calcul de la courbure de l'interface, il est préférable d'avoir pour tout $m \in \mathcal{M}$, $card(N_{[1,m]} > 4)$ (Annexe B). Il en découle donc la condition (2.20) pour l'arête *i* de *e*.

$$
\mathrm{card}(\mathcal{N}_{[1,m_{[e,i]})}) > 4 \text{ et } \mathrm{card}(\mathcal{N}_{[1,m_{[e,i_0]}]}) > 4 \qquad (2.20)
$$

**Éviter le renversement des triangles** : Comme mentionné dans la section 2.4, page 56, il faut faire attention au renversement des triangles lors des différentes procédures. En effet,
dans certaines configurations, l'échange d'arête peut entraîner un renversement des triangles.
Un exemple est illustré sur la figure 2.14. Après transformation, le triangle *e* se trouve avec sa

<figure><img src="image_24.png" /><figcaption>FIGURE 2.14 – L’arête 1 (rouge) de l’élément *e* est échangée. Le triangle *e'* voit sa normale inchangée après l’opération. Ce qui n’est pas le cas du triangle *e*, dont la normale pointe après vers l’intérieur.</figcaption></figure>

normale pointant en sens inverse de la normale extérieure. Cette situation est à éviter. Ainsi, en comparant les normales entre l'instant $t^{(n)}$ et l'état hypothétique $t^{(*)}$ après modification, nous devons avoir :

$$ \bar{n}_e^{(n)} \cdot \bar{n}_e^{(*)} > 0 \qquad (2.21) $$

**Qualité des triangles** : le dernier critère concerne la qualité $Q_e$ des triangles. Elle quantifie le fait qu'un triangle est plus ou moins équilatéral. Définie de diverses façons, nous choisissons celle proposée par ROGHAIR et al. [99]. Cette définition se base sur des grandeurs facilement accessibles, à savoir les longueurs des arêtes. Elle se calcule en faisant le rapport entre le diamètre du cercle inscrit $D_i$ et le rayon du cercle circonscrit $R_c$ :

$$ Q_e = D_i/R_c \qquad (2.22) $$

Pour l'élément $e \in \mathcal{E}$, elle peut se réécrire en se servant des longueurs $d_{e,i}$ de ses arêtes (voir section 2.1, page 47), avec $i \in [1, 3]$ :

$$ Q_e = \frac{(d_{e,1} + d_{e,2} - d_{e,3}) \times (d_{e,2} + d_{e,3} - d_{e,1}) \times (d_{e,3} + d_{e,1} - d_{e,2})}{d_{e,1} \times d_{e,2} \times d_{e,3}} \qquad (2.23) $$

Ainsi, quand le triangle *e* est équilatéral, toutes les longueurs sont égales et $Q_e = 1$. A l'inverse, dans le cas d'un triangle plat, le côté le plus long est égal à la somme des deux autres. Le numérateur s'annule donc et $Q_e = 0$. Un des objectifs de l'échange d'arête est de pouvoir améliorer cette qualité. Soient *e* et *e'*, les deux triangles devant échanger une arête. La moyenne de leur qualité après opération doit être supérieure à celle avant l'échange :

$$ \left[ \left( Q_e^{(*)} + Q_{e'}^{(*)} \right) / 2 \right]^{1/k} > \left[ \left( Q_e^{(n)} + Q_{e'}^{(n)} \right) / 2 \right]^{1/k} \qquad (2.24) $$

Notre choix s'est portée pour une moyenne harmonique ($k = -1$).

Ainsi, lorsque toutes les conditions sont réunies Eqs. (2.19), (2.20), (2.21), (2.24), nous échangeons les arêtes de *e* et *e'* et mettons à jour les voisins. Les différentes étapes sont schématisées dans l'algorithme 2.5.
1 $\hat{e} = e$ ! Création temporaire de $\hat{e}$ pour stocker $e$

2 $\hat{e}' = e'$ ! Création temporaire de $\hat{e}'$ pour stocker $e'$

3 ! Échange de l'arête $i$ entre les voisins $e$ et $e'$

4 $m_{[e,i]} = m_{[\hat{e}'{,i_2}]}$

5 $m_{[e',i]} = m_{[\hat{e}',i_2]}$

6 ! Mise à jour de la connectivité de $e$ et $e'$

7 $e_{[i]} = \hat{e}'_{[i_2]}$ ; $e_{[i_2]} = e'$

8 $e'_{[i]} = \hat{e}_{[i_2]}$ ; $e'_{[i_2]} = e$

9 ! Mise à jour de la connectivité des voisins de $e$ et $e'$

10 $\tilde{e} = \hat{e}_{[i_2]}$ ; $\exists k \mid \tilde{e}_{[k]} = \hat{e}$ ; $\tilde{e}_{[k]} = e'$

11 $\tilde{e} = \hat{e}'_{[i_2]}$ ; $\exists k \mid \tilde{e}_{[k]} = \hat{e}'$ ; $\tilde{e}_{[k]} = e$

12 ! Suppression des éléments de stockage

13 Suppression de $\hat{e}$ et $\hat{e}'$

<caption>Listing 2.5 – Échange d'arête</caption>

## 2.5.2 Lissage

Toujours pour améliorer la régularité du maillage de l'interface, nous avons mis en place une procédure de lissage. Cette dernière peut se faire point par point, arête par arête ou une combinaison des deux, comme présenté par KUPRAT et al. [59]. Nous avons opté pour un lissage laplacien, point par point. Soit $m \in \mathcal{M}$, de coordonnées $\bar{x}$. Nous considérons son premier voisinage, $\mathcal{N}_{[1,m]}$, de cardinal $card(\mathcal{N}_{[1,m]})$. Nous supposons également que $m$ occupe la même place $i$ dans tout élément $e \in \mathcal{N}_{[1,m]}$. En d'autres termes :

$$ \forall e \in \mathcal{N}_{[1,m]}, m_{[e,i]} = m $$

Cela est possible car, pour $n_e = 3$, le choix du premier marqueur de $e$ est arbitraire. La nouvelle position de $m$ après lissage est donc donnée par l'équation (2.25).

$$ \forall m \in \mathcal{M}, \bar{x} = \frac{1}{2 \times \mathrm{card}(\mathcal{N}_{[1,m]})} \sum_{e \in \mathcal{N}_{[1,m]}} \sum_{\substack{j \in [1,3] \\ j \neq i}} \bar{x}_{e,j} \quad (2.25) $$

En d'autres termes, le marqueur $m$ est placé au barycentre des marqueurs constituant son premier voisinage, comme présenté sur la figure 2.15.

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

<figure><img src="image_26.png" /><figcaption>FIGURE 2.15 – Le lissage met $m_{[e,1]}$ au barycentre des marqueurs de son premier voisinage.</figcaption></figure>
### 2.5.3 Suppression de mauvais triangles

Afin de s'assurer que le maillage de l'interface est relativement de bonne qualité, une procédure supplémentaire est implémentée. Elle se focalise essentiellement sur deux types de triangles, qui peuvent restés dans le maillage malgré l'application des précédentes procédures. Nous distinguons ainsi d'une part les triangles en aiguille et les triangles plats d'autre part. Une illustration de ces deux types de triangles est présentée sur la figure 2.16. Le traitement effectué est différent selon la nature du triangle. Nous désignons par $\theta_{e,i}$ l'angle formé au marqueur $m_{[e,i]}$ de l'élément e.

<figure><img src="image_27.png" /><figcaption>FIGURE 2.16 – Un élément en aiguille (a) et un élément plat (b) à retirer du maillage lagrangien.</figcaption></figure>

**Triangle plat** : un élément plat possède un angle largement supérieur aux deux autres, résultant ainsi en une forme de chapeau. Un élément e est considéré comme plat s'il vérifie la relation (2.26).

$$ \max_{i \in [1,3]} (\theta_{e,i}) \geq 8\pi/9 \qquad (2.26) $$

Dans ce cas, un échange d'arête est opéré (Algo. 2.5), avec seulement les critères (2.20) et (2.21) pris en compte.

**Triangle en aiguille** : un élément e est qualifié d'aiguille si son plus petit angle est négligeable devant les deux autres, formant ainsi une pointe. De plus, sa plus petite arête est négligeable face à son arête la plus longue. La relation (2.27) permet de qualifier un élément $e \in \mathcal{E}$ de triangle en aiguille.

$$ \min_{i \in [1,3]} (\theta_{e,i}) \leq 8\pi/9 \quad \text{et} \quad \min_{i \in [1,3]} (d_{e,i}) / \max_{i \in [1,3]} (d_{e,i}) \leq 0.3 \qquad (2.27) $$

L'attitude adoptée est le désenrichissement de l'élément e (Algo. 2.4) en respectant seulement les conditions définies par les relations (2.16) et (2.18).

## 2.6 Conservation du volume

L'écoulement étant incompressible, nous devons nous assurer de la conservation du volume $\Omega$ au cours du temps. Pour cela, il nous faut le calculer à partir de la connaissance du maillage de l'interface.

Le volume discret $V^{(n)}$ délimité par l'interface $\partial\Omega$ peut être décomposé en plusieurs sous-volumes $V_{e}^{(n)}$, contenus dans des tétraèdres. En effet, ces tétraèdres sont constitués des triangles de l'interface et d'un point de référence. Ce point ne doit pas être trop éloigné des marqueurs de l'interface pour avoir des grandeurs du même ordre lors du calcul des volumes. Nous utilisons
<figure><img src="image_28.png" /><figcaption>FIGURE 2.17 – Le volume se calcule en faisant la somme des volumes des tétraèdres construits à partir des triangles de l’interface et d’un point de référence.</figcaption></figure>

le barycentre G des marqueurs de l'interface, de coordonnées $\bar{x}_G$ (Eq. (2.28)), comme point de référence pour construire les tétraèdres (Fig. 2.17). Celui-ci est facile d'accès. Une alternative serait le centre géométrique ou de masse des marqueurs de l'interface.

$$ \bar{x}_G = \frac{1}{\mathrm{card}(\mathcal{M})} \sum_{m \in \mathcal{M}} \bar{x}_m \qquad (2.28) $$

Ainsi, le volume algébrique du tétraèdre, construit à partir de l'élément e et du barycentre G, est :

$$ \mathcal{V}_e^{(n)} = ((\bar{x}_{e,1}^{(n)} - \bar{x}_G) \wedge (\bar{x}_{e,2}^{(n)} - \bar{x}_G)) \cdot (\bar{x}_{e,3}^{(n)} - \bar{x}_G) / 6 \qquad (2.29) $$

Nous en déduisons le volume de la bulle à l'instant $t^n$ :

$$ \mathcal{V}^{(n)} = \sum_{e \in \mathcal{E}^{(n)}} \mathcal{V}_{e}^{(n)} \geq 0 \qquad (2.30) $$

Afin de conserver le volume, les marqueurs vont être déplacés. Ce déplacement se fait dans la direction de la normale à l'interface, proportionnellement à la vitesse de celle-ci. Ainsi, pour un marqueur $m_{[e,i]}$, de vitesse $\bar{v}_{e,i}$, de normale unitaire $\bar{n}_{e,i}$, le déplacement se fait proportionnellement au vecteur :

$$ \overline{dl}_{e,i} = |\bar{v}_{e,i} \cdot \bar{n}_{e,i}| \bar{n}_{e,i} \qquad (2.31) $$

Ainsi, en notant (*) l'état de l'interface avec la correction de volume, la nouvelle position de tout marqueur $m_{[e,i]}$ est :

$$ \bar{x}_{e,i}^{(*)} = \bar{x}_{e,i}^{(n)} + \alpha \times \overline{dl}_{e,i}^{(n)}, \text{ avec } \alpha \in \mathbb{R} \qquad (2.32) $$

Déterminer $\alpha$ permet de trouver la correction à apporter à chaque marqueur pour conserver le volume. La recherche de $\alpha$ se fait en réécrivant le volume $\mathcal{V}$ après correction de la position des marqueurs. Ainsi, en considérant l'élément e, le volume du tétraèdre après correction s'écrit :

$$ \mathcal{V}_e^{(*)} = \mathcal{V}_e^{(n)} + (A_e \alpha^3 + B_e \alpha^2 + C_e \alpha) / 6 \qquad (2.33) $$

Les coefficients $A_e, B_e, C_e$ sont égaux à :

$$ \begin{align} A_e &= (\overline{dl}_{e,1}^{(n)} \wedge \overline{dl}_{e,2}^{(n)}) \cdot \overline{dl}_{e,3}^{(n)} \nonumber \\ B_e &= (\overline{dl}_{e,2}^{(n)} \wedge \overline{dl}_{e,3}^{(n)}) \cdot (\bar{x}_{e,1}^{(n)} - \bar{x}_G) + (\overline{dl}_{e,3}^{(n)} \wedge \overline{dl}_{e,1}^{(n)}) \cdot (\bar{x}_{e,2}^{(n)} - \bar{x}_G) + \nonumber \\ & \qquad (\overline{dl}_{e,1}^{(n)} \wedge \overline{dl}_{e,2}^{(n)}) \cdot (\bar{x}_{e,3}^{(n)} - \bar{x}_G) \tag{2.34} \\ C_e &= ((\bar{x}_{e,2}^{(n)} - \bar{x}_G) \wedge (\bar{x}_{e,3}^{(n)} - \bar{x}_G)) \cdot \overline{dl}_{e,1}^{(n)} + ((\bar{x}_{e,3}^{(n)} - \bar{x}_G) \wedge (\bar{x}_{e,1}^{(n)} - \bar{x}_G)) \cdot \overline{dl}_{e,2}^{(n)} + \nonumber \\ & \qquad ((\bar{x}_{e,1}^{(n)} - \bar{x}_G) \wedge (\bar{x}_{e,2}^{(n)} - \bar{x}_G)) \cdot \overline{dl}_{e,3}^{(n)} \nonumber \end{align} $$
Le volume délimité par l'interface est alors :

$$ \nu^{(*)} = \nu^{(n)} + (A\alpha^3 + B\alpha^2 + C\alpha) / 6 \quad (2.35) $$

avec

$$ A = \sum_{e \in \mathcal{E}^{(n)}} A_e, \quad B = \sum_{e \in \mathcal{E}^{(n)}} B_e, \quad C = \sum_{e \in \mathcal{E}^{(n)}} C_e \quad (2.36) $$

La conservation du volume s'écrit donc :

$$ \nu^{(*)} = \nu^{(0)} \iff A\alpha^3 + B\alpha^2 + C\alpha + 6\nu^{(n)} - 6\nu^{(0)} = 0 \quad (2.37) $$

L'inconnue $\alpha$ est donc solution d'une équation polynomiale de degré 3. Cette équation non-linéaire peut être résolue de manière itérative par une méthode de Newton ou de manière directe par la méthode de Cardan (voir Annexe A). Après la résolution de l'équation, trois solutions se présentent à nous. La solution $\alpha_0$ que nous retenons est la solution réelle la plus petite en valeur absolue. Ainsi, grâce à cette méthode, nous pouvons contrôler de manière exacte le volume discret délimité par le maillage de l'interface.

## 2.7 Validation - Interface sphérique dans un champ tourbillonnaire

Ce test a été proposé par LEVEQUE [65] et repris par ENRIGHT et al. [30]. Toutes les longueurs sont adimensionnées par la taille du domaine de simulation. Ainsi, dans une boîte cubique de côté $[0, 1]^3$ est placé une interface sphérique de rayon 0.15 et de centre (0.35, 0.35, 0.35). Elle est soumise à un champ de vitesse à divergence nulle définie comme suit :

$$ \begin{align} u(x, y, z, t) &= 2 \sin^2(\pi x) \sin(2\pi y) \sin(2\pi z) \cos(\pi t/T) \nonumber \\ v(x, y, z, t) &= -\sin(2\pi x) \sin^2(\pi y) \sin(2\pi z) \cos(\pi t/T) \tag{2.38} \\ w(x, y, z, t) &= -\sin(2\pi x) \sin(2\pi y) \sin^2(\pi z) \cos(\pi t/T) \nonumber \end{align} $$

avec, $T = 3$. Cet écoulement engendre une déformation à la fois dans les plans $(x; y)$ et $(x; z)$. Ce champ est temporellement périodique et réversible. Ainsi, à $t = T/2$, le champ de vitesse s'annule et l'interface est au maximum de sa déformation. Le champ de vitesse change alors de sens, ramenant l'interface à sa position initiale à $t = T$.

D'abord, nous effectuons des tests de nature qualitative. En effet, nous essayons dans un premier temps d'avoir des interfaces graphiquement semblables à celles présentes dans la littérature [30, 44]. Sur la figure 2.18, différents instants du maillage de l'interface au cours de son transport sont illustrés. La boîte cubique possède 80 cellules par direction, résultant en un pas d'espace $h = 1/80$. Le pas de temps est $\Delta t = 10^{-3}$. Graphiquement, nous obtenus des résultats très satisfaisant, quand nous les comparons à ceux de ENRIGHT et al. [30], qui utilisent la méthode 'semi-Lagrangian particle level set', avec un pas d'espace effectif $h = 1/512$ (Fig.2.19).

La suite des tests s'est portée sur un volet beaucoup plus quantitatif. En effet, nous nous sommes intéressés à la convergence spatiale de la position des marqueurs de l'interface en fonction de la discrétisation de la boîte cubique. En effet, nous rappelons que la longueur $d$ des arêtes des éléments du maillage de l'interface est liée à la taille $h$ des cellules de la boîte cubique (voir relations (2.11) et (2.13)).

$$ h/4 < d < h $$
<figure><img src="image_29.png" /><figcaption>FIGURE 2.18 – Évolution d'une interface sphérique dans un champ tourbillonnaire, avec $\Delta t = 10^{-3}$ et $h = 1/80$. De gauche à droite, l'interface est représentée aux instants $0, T/10, T/5$ et $T/2$, le maximum de déformation.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE 2.19 – Différents instants de l'interface à $0, T/15, 2T/15, 3T/15, 4T/15, 5T/15, 6T/15, 7T/15, 9T/15, 11T/15, 13T/15$ et $T$, tirés de [30].</figcaption></figure>

Pour cela, il nous faut mesurer l'écart entre notre solution numérique discrète et la solution de référence. N'ayant pas de solution analytique de référence pour le calcul des erreurs, nous avons donc construits une solution de référence numérique, obtenue pour un pas d'espace $h_{ref} = 1/640$ et un pas de temps $\Delta t_{ref} = 10^{-4}$. Ainsi, pour tout marqueur $m \in \mathcal{M}$, cet écart mesure le rayon de la plus petite boule centrée en $m$ et qui touche l'interface de référence (Fig. 2.20). Nous notons $I_m$ le point de contact associé au marqueur $m \in \mathcal{M}$, de coordonnées $\bar{x}_{I_m}$. Ce point vérifie la relation suivante :

$$I_m = \underset{P \in \partial \Omega_{ref}}{\operatorname{argmin}} \| \bar{x}_m - \bar{x}_P \|$$

Nous pouvons donc écrire que cet écart, noté $r_m$, est égal à :

$$r_m = \| \bar{x}_m - \bar{x}_I \| \qquad (2.39)$$
<figure><img src="image_31.png" /><figcaption>FIGURE 2.20 – Écart entre un marqueur du maillage grossier et l'interface du maillage de référence.</figcaption></figure>

Plus concrètement, $r_m$ est la plus petite des distances euclidiennes entre $m \in \mathcal{M}$ et tout élément $e \in \mathcal{E}_{ref}$. Trouver ce point de contact $I_m$ pour évaluer $r_m$ se fait en deux étapes.

**Trouver les éléments** $e \in \mathcal{E}_{ref}$ les plus proches de $m \in \mathcal{M}$ : Parcourir tous les éléments $e \in \mathcal{E}_{ref}$ pour trouver $I_m$ associé à $m \in \mathcal{M}$ est prohibitif. Soit le marqueur $m \in \mathcal{M}$ localisé dans la cellule d'indices $(i_0, j_0, k_0)$. Nous recherchons de manière itérative les couronnes, autour de la cellule initiale, qui contiennent des éléments de $\mathcal{E}_{ref}$. Ainsi, nous regardons dans les cellules d'indices $(i_0 \pm i, j_0 \pm j, k_0 \pm k)$, avec $(i, j, k) \in \mathbb{N}^3$. Tant qu'aucun élément $e \in \mathcal{E}_{ref}$ n'est rencontré, $i, j$ et $k$ sont incrémentés de 1. Lorsque la première couronne contenant des éléments $e \in \mathcal{E}_{ref}$ est trouver, nous prenons également la couronne suivante pour nous assurer d'avoir la bonne distance. En effet, comme le montre la figure 2.21, il peut se trouver des situations où l'écart minimal se localise sur la couronne d'après.

<figure><img src="image_32.png" /><figcaption>FIGURE 2.21 – Le marqueur $m$ se trouve dans la cellule $(i_0, j_0, k_0)$. Les premiers éléments $e \in \mathcal{E}_{ref}$ rencontrés sont localisés dans la première couronne (zone blanche). Cependant, la distance $r_m$ établie avec ces éléments est plus grande que celle construite avec les éléments de la deuxième couronne (zone verte), $r'_m$.</figcaption></figure>

**Projection orthogonal du marqueur** $m$ sur le plan contenant l'élément $e \in \mathcal{E}_{ref}$ :
nous notons $H$ le projeté orthogonal du marqueur $m$ sur le plan ($\mathcal{P}_e$) contenant l'élément $e \in \mathcal{E}_{ref}$.
La position de $H \in (\mathcal{P}_e)$ par rapport à l'élément $e$ nous permet de déterminer $I_m$, et donc la
distance euclidienne entre le marqueur $m$ et l'élément $e$. Trois résultats sont possibles, comme
l'illustre la figure 2.22. Nous distinguons trois zones :
<figure><img src="image_33.png" /><figcaption>FIGURE 2.22 – Le point d’intersection entre la droite partant de $m$, de vecteur directeur $n_e$, et la plan ($P_e$) peut se trouver dans trois différentes zones. S’il est dans la zone 1 (intérieur de $e$), alors $I = H_1$. S’il est dans la zone 2, alors $I = \text{projeté orthogonal de } H_2$ sur l’arête. Dans la zone 3, $I$ est un sommet de $e$. $H_3$ est proche de $m_{[e,2]}$ donc $I = m_{[e,2]}$.</figcaption></figure>

— zone 1 : c’est l’intérieur de l’élément $e$. Dans ce cas, $I = H_1$.

— zone 2 : ce sont les 3 régions du plan, à l’extérieur de $e$, comprises entre les 2 droites perpendiculaires aux 3 arêtes, passant par leurs extrémités (Fig. 2.22). Dans cette situation, $I$ est le projeté orthogonale de $H_2$ sur l’arête la plus proche.

— zone 3 : ce sont les 3 régions du plan, comprises entre les 2 droites perpendiculaires à 2 arêtes adjacentes et passant par leur extrémité commune. Dans ces parties, $I$ correspond au sommet de $e$ le plus proche. Ainsi, sur la figure 2.22, $I = m_{[e,2]}$.

De la mesure de la distance des marqueurs $m$ à la surface de référence $\partial\Omega_{ref}$, nous pouvons construire une mesure des erreurs en norme $L_2$ et $L_\infty$, présentée dans l’équation (2.40).

$$Err_2(h) = \left( \sum_{m \in \mathcal{M}} \frac{r_m^2}{N_m} \right)^{1/2} \quad , \quad Err_\infty(h) = \max_{m \in \mathcal{M}} r_m \quad (2.40)$$

La figure 2.23 présente l’évolution de l’erreur entre l’interface $\partial\Omega(h)$ et $\partial\Omega_{ref}$ pour 2 normes, à $t = T/2$. Nous observons une convergence à l’ordre 2 pour $Err_2(h)$. La convergence tombe à 1 quand nous nous intéressons à $Err_\infty(h)$. Nous constatons qu’avec l’ensemble de nos procédures pour gérer le maillage de l’interface lors de son transport et de sa déformation, nous arrivons à avoir une bonne précision sur la position des marqueurs. De plus, l’ordre 2 constaté sur $Err_2(h)$ nous permet de pouvoir réduire l’erreur spatiale sur les marqueurs. Ainsi, en divisant par 2 la taille de la grille eulérienne, nous divisons par 4 l’erreur sur la position des marqueurs.
<figure><img src="image_34.png" /><figcaption>FIGURE 2.23 – Erreurs spatiales en norme $L_2$ (triangle) et $L_\infty$ (carré) en fonction de la taille $h$ des cellules de discrétisation du domaine.</figcaption></figure>
# Bibliographie

## Gestion du maillage de l'interface

[30] D. ENRIGHT, F. LOSASSO et R. FEDKIW. "A Fast and Accurate Semi-Lagrangian Particle Level Set Method". In: *Computers & Structures* 83.6 (fév. 2005), p. 479-490 (cf. p. 50, 63, 64, 145).

[33] M. FEBRES et D. LEGENDRE. "Enhancement of a 2D front-tracking algorithm with a non-uniform distribution of Lagrangian markers". In: *J. Comput. Phys.* 358 (2018), p. 173-200 (cf. p. 55).

[44] J. P. GOIS, A. NAKANO, L. G. NONATO et G. C. BUSCAGLIA. "Front tracking with moving-least-squares surfaces". In: *J. Comput. Phys.* 227 (2008), p. 9643-9669 (cf. p. 63).

[59] KUPRAT, ANDREW, KHAMAYSEH, AHMED, GEORGE, DENISE, LARKEY et LEVI. "Volume Conserving Smoothing for Piecewise Linear Curves, Surfaces, and Triple Lines". In: *J. Comput. Phys.* 172.1 (sept. 2001), p. 99-118 (cf. p. 60).

[65] R. J. LEVEQUE. "High-resolution conservative algorithms for advection in incompressible flows". In: *SIAM J. Numer. Anal.* 33.2 (avr. 1996), p. 627-665 (cf. p. 50, 63, 145).

[74] R. MCDERMOTT et S. POPE. "The parabolic edge reconstruction method (PERM) for Lagrangian particle advection". In: *J. Comput. Phys.* 227.11 (2008), p. 5447-5491 (cf. p. 49, 89, 92, 94).

[99] I. ROGHAIR, V. S. ANNALAND et J. KUIPERS. "An improved Front-Tracking technique for the simulation of mass transfer in dense bubbly flows". In: *Chem. Eng. Sci.* 152 (oct. 2016), p. 351-369 (cf. p. 50, 55, 59, 89).

[118] M. TAVARES, D.-A. KOFFI BI, E. CHENIER et S. VINCENT. "A Two-Dimensional Second Order Conservative Front-Tracking Method with an Original Marker Advection Approach Based on Jump Relations". In: *Communications in Computational Physics* 27.5 (2020), p. 1550-1589 (cf. p. 37, 45).

[120] D. J. TORRES et J. U. BRACKBILL. "The Point-Set Method: Front-Tracking without Connectivity". In: *J. Comput. Phys.* 165 (2000), p. 620-644 (cf. p. 29, 30, 47).
# Chapitre 3

## Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces

### Sommaire

<table><tr><td>3.1</td><td>Introduction</td><td>72</td></tr><tr><td>3.2</td><td>Num. methods for eval. κ, n and f</td><td>74</td></tr><tr><td>3.2.1</td><td>Main mathematical notations</td><td>74</td></tr><tr><td>3.2.2</td><td>Discretization of the Laplace-Beltrami Operator (LBO-methods)</td><td>75</td></tr><tr><td>3.2.3</td><td>Integral Formulation of the surface tension (IF-methods)</td><td>76</td></tr><tr><td>3.2.4</td><td>Surface Reconstruction (SR-methods)</td><td>78</td></tr><tr><td>3.2.5</td><td>Discrete surfaces and solution accuracy</td><td>79</td></tr><tr><td>3.3</td><td>Analysis of a stationary surface</td><td>80</td></tr><tr><td>3.3.1</td><td>Results</td><td>81</td></tr><tr><td>3.3.2</td><td>Comparisons</td><td>85</td></tr><tr><td>3.3.3</td><td>Summary</td><td>87</td></tr><tr><td>3.4</td><td>Analysis of advected surfaces</td><td>88</td></tr><tr><td>3.4.1</td><td>Problem presentation</td><td>88</td></tr><tr><td>3.4.2</td><td>Ligament formation – Case Q &lt; 0</td><td>90</td></tr><tr><td>3.4.3</td><td>Bubble shape – Case Q &gt; 0</td><td>93</td></tr><tr><td>3.4.4</td><td>Summary</td><td>96</td></tr><tr><td>3.5</td><td>Conclusion</td><td>96</td></tr></table>
# Désir-André Koffi-Bi<fnref n="1" />, Mathilde Tavares, Éric Chénier, Stéphane Vincent

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

## Abstract :

In this work, three classes of numerical methods are investigated to evaluate the mean curvature, the unit normal vector and the surface tension on a front tracking interface encountered in the simulation of multiphase flows with separated phases. The Laplace-Beltrami Operator discretization, the Integral Formulation of the surface tension and the Surface Reconstruction technique are well-known methods used in the literature, but whose accuracy, robustness and convergence properties are seldom studied. In a first step, different variants of these methods are presented and compared against each other on a static analytical surface to measure their sensitivity to the size and regularity of the mesh. Then, to assess the influence of the surface advection scheme and the remeshing procedures, two original and time dependent analytical surfaces have been developed, leading to a ligament formation or the birth of a drop/bubble on a flat surface/liquid film. Comparisons to such dynamical surfaces are especially useful and significant, since they highlight the sensitivity of numerical methods for the interface property calculation to the errors produced by the Lagrangian transport of the front-tracking surface. Finally, the accuracy of the different approximations are correlated to their computation time, so that the user may choose the most appropriate method according to the desired accuracy and cost.

## Keywords : Front-tracking, multiphase flow, surface tension, curvature

## 3.1 Introduction

Multiphase flows are very common in industry and natural environment. Amongst the numerous examples, it can be mentioned the heat and mass transfer for safety issues in nuclear power plants or for the control of the chemical composition in alloy manufacturing, or also the liquid atomization, either of fuel in engines to reduce pollutants generated by combustion [13], or of sea water resulting from the breaking wave to understand the coupling between the ocean and the atmosphere for meteorology purpose. For these flows, the heat and mass transfers result from the energy and momentum balances in each phase, but also from the physics at the interface which is governed by the jump conditions [21]

$$[\rho(\mathbf{u} - \mathbf{V}_i)] \cdot \mathbf{n} = \dot{m} \qquad (3.1)$$

$$[-\dot{m}\mathbf{u} - p\mathbf{n} + \mu\mathbf{n} \cdot (\nabla\mathbf{u} + (\nabla\mathbf{u})^t)] = \sigma\kappa\mathbf{n} + (\nabla\sigma) \qquad (3.2)$$

$$[\lambda \cdot \nabla T] \cdot \mathbf{n} = \dot{m}L_f \qquad (3.3)$$

where [●] stands for the difference of the quantity ● on both sides of the interface, in the direction of the unit normal vector to the interface **n**, pointing outward. It is worth noticing that the forces and fluxes acting at the interface depend on the one hand, to fluid flow parameters, such

<footnote n="1">1. desir-andre.koffi-bi@univ-eiffel.fr</footnote>
as the fluid velocities **u**, the velocity of the interface **V**ᵢ, the pressures **p**, the temperatures **T**,
the densities **ρ**, the kinematic viscosities **μ**, the surface tension **σ** and the latent heat **L**f, on the
other hand to geometrical properties such as twice the mean curvature **κ** and the unit normal
vector **n** to the interface.

Since the shape of the interfaces, allowing the calculation of κ and **n**, is one of the very
important parameters to control the heat and mass transfers, the numerical models must be able
to approximate it accurately. However, such a task can turn out to be very difficult, especially
when a large range of interfacial scales have to be modelled. For example, in the atomization
process, the cascade of interface scales must be solved, which requires a very expensive mesh
refinement : the initial breakup of the continuous liquid phase gives rise to individual liquid
parcels before producing a dispersed spray. Even though the establishment of the drop size
distribution was successfully described by macro- or meso-scale models with ELSA or CEDRE
solvers [22, 139], these descriptions are only valid once the spray is formed. Therefore, the
primary atomization is still poorly understood and needs new modelling strategies. Computing
the curvature and more generally the surface properties also arise in boundary integral methods
where volumetric problems are transformed into surface problems [136].

The understanding and the control of the heat and mass transfer in multiphase flows also
requires that we distinguish the different interface shapes. Thanks to the databases stemming
from CLS-VOF schemes [125], Massot and co-workers [13, 77] have been able to propose new
ways to analyse the interface dynamics during primary atomization. They have demonstrated
that each type of structure (spheres, ligaments or sheets) is characterized by specific distribu-
tions of the interface curvature. Therefore, an accurate numerical evaluation of the interface
curvature is essential to ensure the reliability of the results.

In our contribution, we are interested in the numerical evaluation of the topological pro-
perties of interfaces, namely the mean curvature κ and the unit normal vector **n**, but also the
surface tension **f**. Depending on the model bieng used, the capillary effects are either introduced
in the Navier-Stokes equations [108], or taken into account through the jump conditions (Eqs.
3.1-3.3) by the geometrical components κ and **n**, or the surface force **f**. In the first approach,
the surface tension contribution is modelled by a volume force depending on the volume fraction
derivatives expressed over the fluid mesh [10]. Another approach relies on the writing of the
jump conditions directly onto the interface approximation. This latter method turns out to be
more accurate because it does not a-priori suffer from a spreading of the capillary force over
the fluid mesh : calculations are carried out precisely at the surface location. The ghost fluid
method [34] and the front-tracking approach [26] belong to these techniques. However, when the
interface does not fit to the nodes where the conservation equations for velocity, temperature
and pressure are solved, the local capillary information must be distributed on the surrounding
fluid nodes, which produces a smoothing effect of the surface tension, again [34, 108].

The present work is focused on the numerical evaluation of the mean curvature κ, the unit
normal vector **n** and the surface tension **f**, in the framework of a discrete interface made of a
triangle tessellation. Coming either from the computational fluid dynamics domain, for example
in the front-tracking method framework, or from the computer aided design, only few studies
have dealt with their convergence properties and robustness [67]. Moreover, no comparison
between a large panel of numerical methods exists in the literature. Lastly, the only available
validations concern static surfaces where the nodes of the discrete surface lie exactly on the
reference surface.

In our contribution, three large families of methods are studied to get the numerical approxi-
mations. In the framework of the computational fluid dynamics, we have developed a original
time dependent analytical surface, that qualitatively gives rise to characteristic patterns met
in multiphase flows, namely the formation of ligaments or the birth of a drop/bubble on a flat surface/liquid film. These new test cases have allowed us, for the first time, to examine the effect of the advection scheme, with remeshing procedures, onto the accuracy of $\kappa$, $\mathbf{n}$ and $\mathbf{f}$ approximations.

The paper is organized as follows. The different numerical approximations of $\kappa$, $\mathbf{n}$ and $\mathbf{f}$, which are stemmed from the three families of methods, are first briefly reminded (Sec. 3.2), and then compared against each other for a complex static analytical surface (Sec. 3.3). Section 3.4 is then devoted to comparisons for the time dependent surfaces, first for an exact advection but with a remeshing procedure, and then using the numerical transport scheme. Conclusions and perspectives are finally drawn in the last section.

## 3.2 Numerical methods for evaluating the curvature $\kappa$, the unit normal vector $\mathbf{n}$ and the capillary force $\mathbf{f}$

The aim of this section is to define the numerical methods used in this paper to approximate the local curvature $\kappa$, the unit normal vector $\mathbf{n}$ and the capillary force $\mathbf{f}$. These evaluations rely on a discrete representation of the interface which consists of a set of planar triangular elements, defined by their vertices, and the knowledge of the neighbours through the edges of each triangle. This kind of representation is widely used in the framework of front tracking methods.

Three main families of approximations are identified in the literature. They are based either on

- the discretization of the Laplace-Beltrami Operator,

- an Integral Formulation to get the surface tension over a specific area,

- a reconstruction method of the interface by an analytical surface.

The different methods devoted to the calculation of the interface properties, and which belong to the three families presented here above, are detailed in the following parts. But first, some notations are introduced.

### 3.2.1 Main mathematical notations

The discrete approximation of the interface is made up of a family $\mathcal{T}$ of triangles. Let $\mathcal{X}$ be the set of vertices of the whole triangles. The vertices of any $T \in \mathcal{T}$ are simply noted $\delta^2 T \subset \mathcal{X}$. So, we refer to

- the coordinate of any vertex $x \in \mathcal{X}$ by $\mathbf{x}$.

- the surface of any triangle $T \in \mathcal{T}$ by $S_T$.

We also introduce

- the set whose elements are couples of vertices $\mathcal{X}_2 = \{(x_\alpha, x_\beta) \in \mathcal{X}^2 | \exists T \in \mathcal{T}, \{x_\alpha, x_\beta\} \subset \delta^2 T\}$.

- the first and second neighbourhood of any vertex $x \in \mathcal{X}$ by $\mathcal{N}_1(x)$ and $\mathcal{N}_2(x)$, with

$$ \mathcal{N}_1(x) = \{x' \in \mathcal{X} \mid (x, x') \in \mathcal{X}_2\} \qquad (3.4) $$

$$ \mathcal{N}_2(x) = \{x' \in \mathcal{X} \setminus \mathcal{N}_1(x) \mid \exists x'' \in \mathcal{N}_1(x), (x', x'') \in \mathcal{X}_2\} \cup \mathcal{N}_1(x) \qquad (3.5) $$

To express the numerical approximations of the curvature, normal vector and capillary force, the endpoints of the triangle edges have to be ordered. This order is closely linked to the notion of the exterior normal vector, noted $\mathbf{n}$. For a closed surface, it is easily defined. As for an open one, a side of the surface is arbitrary chosen as facing "exterior", and $\mathbf{n}$ is set pointing towards this region. Then, we define
— the orientation of any triangle $T \in \mathcal{T}$ by $\mathcal{X}_2(T) \subset \mathcal{X}$

$$ \mathcal{X}_2(T) = \{(x_\alpha, x_\beta), (x_\beta, x_\gamma), (x_\gamma, x_\alpha) | \delta^2 T = \{x_\alpha, x_\beta, x_\gamma\}, ((\mathbf{x}_\beta - \mathbf{x}_\alpha) \times (\mathbf{x}_\gamma - \mathbf{x}_\alpha)) \cdot \mathbf{n} > 0\} \quad (3.6) $$

with $\times$ the cross product.

— the unit normal vector to any triangle $T \in \mathcal{T}$ by $\mathbf{n}_T$. So, $\forall T \in \mathcal{T}$, with $\delta^2 T = \{x_\alpha, x_\beta, x_\gamma\}$ and $\mathcal{X}_2(T) = \{(x_\alpha, x_\beta), (x_\beta, x_\gamma), (x_\gamma, x_\alpha)\}$, then

$$ \mathbf{n}_T = \frac{(\mathbf{x}_\beta - \mathbf{x}_\alpha) \times (\mathbf{x}_\gamma - \mathbf{x}_\alpha)}{\|(\mathbf{x}_\beta - \mathbf{x}_\alpha) \times (\mathbf{x}_\gamma - \mathbf{x}_\alpha)\|} \quad (3.7) $$

— an approximation of the unit normal vector to the surface at any vertex $x \in \mathcal{X}$ by $\mathbf{n}_{\mathcal{X}_1}(x)$:

$$ \mathbf{n}_{\mathcal{X}_1}(x) = \frac{\sum_{T \in \mathcal{T}} |x \in \delta^2 T| S_T \mathbf{n}_T}{\left\| \sum_{T \in \mathcal{T}} |x \in \delta^2 T| S_T \mathbf{n}_T \right\|} \quad (3.8) $$

This expression is simply proportional to the linear weighted combination by the triangle surface $S_T$ of the normal vectors $\mathbf{n}_T$ of all the triangles $T$ sharing the vertex $x \in \mathcal{X}$.

### 3.2.2 Discretization of the Laplace-Beltrami Operator (LBO-methods)

The Laplace-Beltrami Operator (LBO) is the generalization of Laplacian operator from Euclidean spaces to curved manifolds. Commonly used in the domain of computer aided design, it allows to get the geometrical properties of the surfaces, but also to smooth them. Different approximations of the LBO exist in the literature (see [23, 76, 91, 116] to cite but a few). In this paper, the retained method is based on that by MEYER et al. [76]. The approximation of the mean curvature normal vector $\mathbf{K}$ is:

$$ \mathbf{K}(x_i) = \frac{1}{2\mathcal{A}_{\text{mixed}}} \sum_{(T,T') \in \mathcal{T}_i^*} (\cot \alpha_{ij} + \cot \beta_{ij}) (\mathbf{x}_j - \mathbf{x}_i) = 2\kappa_i \mathbf{n}_i \quad (3.9) $$

with $\mathcal{T}_i^* = \{T \in \mathcal{T}, T' \in \mathcal{T} | \delta^2 T = \{x_i, x_l, x_j\}, \delta^2 T' = \{x_i, x_j, x_k\}\}$ and $\alpha_{ij} = \widehat{x_i x_l x_j}$ and $\beta_{ij} = \widehat{x_j x_k x_i}$ (Fig. 3.1a), $\kappa_i$ and $\mathbf{n}_i$ are the approximations of the mean curvature and unit normal

<figure><img src="image_35.png" /><figcaption>FIGURE 3.1 – Notations and stencil for the discrete LBO methods [31, 76]</figcaption></figure>

vector to the surface at $x_i$, $\mathcal{A}_{\text{mixed}}$ is the measure of the mixed area (Fig. 3.1b). The shape of the mixed area $\mathcal{A}_{\text{mixed}}$ depends on the topology of the triangles sharing the vertex $x_i$.
<figure><img src="image_36.png" /><figcaption>FIGURE 3.2 – Construction of the mixed area $A_{mixed}$ of the discrete LBO, according to MEYER et al. [76] and ESSADKI [31]</figcaption></figure>

**LBO-Meyer method [76]** In the original method, when the angles of the element are acute, its contribution to $A_{mixed}$ is the Voronoi area (Fig. 3.2a).

It corresponds to the area of the quadrangular element surface defined by the vertex $x_i$, the middle of the two edges adjacent to $x_i$ and the circumcenter $x_c$ of the triangle. For obtuse triangles, the circumcenter is outside the triangle and two cases must be discriminated (Fig. 3.2b, left). If $\overline{x_jx_i x_k} > \pi/2$, the circumcenter is replaced by the middle of the opposite edge. Otherwise, the area is that of the triangle defined by the vertex $x_i$ and each middle of the two edges adjacent to $x_i$ (Fig. 3.2b, right).

**LBO-Essad. method [31]** An alternative approach of the mixed area $A_{mixed}$ can be used. Whereas the formalism does not change for acute triangles (Fig. 3.2a), for an obtuse angle, the local contribution is based on the Voronoi area for which the circumcenter is replaced by the barycenter $x_g$ (Fig. 3.2c).

The approximated surface tension, at vertex $x_i$, is simply expressed by $\mathbf{f}_i = \sigma \mathbf{K}(x_i)$ and the approximated local curvature and normal vector to surface $S$, at node $x_i$, write :

$$ \mathbf{n}_i = \text{sign}(\mathbf{f}_i \cdot \mathbf{n}_{N_1}(x_i)) \frac{\mathbf{f}_i}{\|\mathbf{f}_i\|} , \quad \kappa_i = \frac{1}{2\sigma} \text{sign}(\mathbf{f}_i \cdot \mathbf{n}_{N_1}(x_i)) \|\mathbf{f}_i\| \quad (3.10) $$

### 3.2.3 Integral Formulation of the surface tension (IF-methods)

In the Integral Formulation family, the approximation of an average surface tension is directly computed. Initially proposed in the front-tracking framework by POPINET et al. [94] for 2D cases, this method has been later on extended to 3D cases by TRYGGVASON et al. [122]. Assuming a constant surface tension $\sigma$, the average surface tension over $S$, $\bar{\mathbf{f}}$, is the circulation of the "pulls"
exerted along the border $∂S$ :

$$ \bar{\mathbf{f}} = \frac{\sigma}{S} \int_S \kappa \mathbf{n} dS = \frac{\sigma}{S} \oint_{\partial S} \mathbf{t} \times \mathbf{n} dl \qquad (3.11) $$

where $\mathbf{t}$ is the unit tangent vector to $∂S$ and $\mathbf{n}$ the unit normal vector to $S$.

**IF-Trygg. method** [123] There is an earlier version of this method proposed by [122] where a quadratic fitting is used to evaluate $\mathbf{t}$ and $\mathbf{n}$ in the integral formulation. However we choose the recent one as it is the most used. The discrete surface is centred around the vertex $x_i$ and relies on a set of adjacent triangles (Fig. 3.3).

<figure><img src="image_37.png" /><figcaption>FIGURE 3.3 – Stencil for the original [123] method of the Integral Formulation (IF-Trygg. method)</figcaption></figure>

The surface $S$ is a collection of quadrangular polygons, defined piecewise on each triangle, whose one of the vertices are $x_i$, the middle of the two edges adjacent to $x_i$ and the barycentre $x_g$ of the triangle. After simplification, it is possible to remove the edge nodes and the barycentre, to get a relation based on the mesh vertices belonging to $N_1(x)$:

$$ \mathbf{f}_i = \frac{\sigma}{2S} \sum_{T \in \mathcal{T}_i^*} (\mathbf{x}_k - \mathbf{x}_j) \times \mathbf{n}_T \qquad (3.12) $$

with $\mathcal{T}_i^* = \{T \in \mathcal{T} | \delta^2 T = \{x_i, x_j, x_k\}, (x_j, x_k) \in \mathcal{X}_2(T)\}$, $S = \sum_{T \in \mathcal{T}_i^*} S_T$, and $\mathbf{n}_T$ the unit normal vector to the triangle $T$.

**IF-triang. method** The main difference with the previous method dwells in the shape of the surface. In the triangular version of the Integral Formulation, $S$ covers the whole triangles having $x_i$ as a vertex (Fig. 3.4).

<figure><img src="image_38.png" /><figcaption>FIGURE 3.4 – Stencil for the ‘triangle version’ of the Integral Formulation (IF-triang. method)</figcaption></figure>
The average surface tension is then expressed by

$$ \mathbf{f}_i = \frac{\sigma}{S} \sum_{(T,T') \in \mathcal{T}_i^*} (\mathbf{x}_k - \mathbf{x}_j) \times \frac{S_T \mathbf{n}_{T'} + S_{T'} \mathbf{n}_T}{\|S_T \mathbf{n}_{T'} + S_{T'} \mathbf{n}_T\|} \quad (3.13) $$

with $\mathcal{T}_i^* = \{T \in \mathcal{T}, T' \in \mathcal{T} | x_i \in \delta^2 T, x_i \notin \delta^2 T', \{x_j, x_k\} = \delta^2 T \cap \delta^2 T', (x_j, x_k) \in \mathcal{X}_2(T)\}$, $S_T$ and $S_{T'}$ respectively the surfaces of the triangles T and T', $\mathbf{n}_T$ and $\mathbf{n}_{T'}$ their respective unit normal vectors and $S = \sum_{T \in \mathcal{T}_i^*} S_T$.

Once the surface tension is evaluated, the approximated local curvature and normal vector to the surface S at node $x_i$ are expressed in the same way as for the LBO-methods (see Eq. (3.10)).

### 3.2.4 Surface Reconstruction (SR-methods)

Knowing the analytic equation of a surface enables the exact calculation of its geometrical properties, specifically the mean curvature and the unit normal vector [45]. This equation can be constructed over the entire surface, like what was proposed by POPINET et al. [94] and TRONTIN et al. [121] to approximate 1D interfaces with cubic splines, or defined locally to approximate a part of the interface mesh [26, 138]. The local approach leads obviously to simpler methods, especially for 2D surfaces and complex interface shapes. Whereas the curvature does not depend on the local coordinates, a transition matrix must be used to proceed to the change from the local to the global basis.

For the local reconstruction methods, the local basis originates at the vertex $x_i$ where the interface properties have to be approximated. Several methods are distinguished in this paper depending on

- the expression of the approximate normal vector to the surface at $x_i$,

- the stencil employed to reconstruct the local analytical surface,

- the eventual requirement that the node $x_i$ belongs to the reconstructed surface.

Once the approximate normal vector is defined, the two other vectors of the orthogonal basis can be chosen arbitrarily. In this paper, the reconstructed analytical surface has a polynomial expression and the calculation of the coefficients is the result of a least square method. The parameters of the Surface Reconstruction method are now detailed.

**Stencil for the reconstructed surface – $\mathcal{N}_1(x_i)$ or $\mathcal{N}_2(x_i)$**.

The influence of the stencil on the reconstructed surface is studied. The least square method is defined on mesh vertices that belong to, either the first neighbourhood $\mathcal{N}_1(x_i)$, or the second neighbourhood $\mathcal{N}_2(x_i)$ of the vertex $x_i$. Furthermore, some precautions must be taken into account when using the first neighbourhood $\mathcal{N}_1(x_i)$. Indeed, if the number of vertices is not sufficient to use the least square method, a completion is made by choosing arbitrary vertices in the second neighbourhood [26]. The Appendix D.1 shows that when these additional vertices are aligned with the first neighbourhood, it deteriorates the evaluation of the curvature and normal vector based on the first neighbourhood. Therefore in this peculiar situation, the solution adopted in this work is to automatically switch to a reconstruction based on the second neighbourhood when needed.

**Vertex $x_i$ solution (or not) of the reconstructed surface – $a_{00} \neq 0$ or $a_{00} = 0$**.

The local surface is approximated by the equation $Z = f(X, Y)$, with the polynomial function $f(X, Y)$ defined by:

$$ f(X, Y) = a_{20}X^2 + a_{11}XY + a_{02}Y^2 + a_{10}X + a_{01}Y + a_{00} \quad (3.14) $$

where $a_{ij}$ are unknown real coefficients to be determined, X, Y, Z are the local coordinates relative to the local basis. Depending on the method, the coefficient $a_{00}$ will be set free or imposed to zero as proposed by Du et al. [26] et ZINCHENKO et al. [138]. In this latter case, the reconstructed surface passes exactly by vertex $x_i$.
**Construction of the local basis (e_X, e_Y, n_α), with α ∈ {c, m, iter}**

The local basis is constructed first by defining the local tangent plane (e_X, e_Y), which is characterized by its normal vector. Different approaches exist to evaluate this normal vector.

— A first choice is a linear combination of all the unit normal vector n_T to triangle T, whose vertices belong to the stencil N_1(x_i) :

$$ \mathbf{n}_c(x_i) = \frac{\sum_{T \in \mathcal{T} | x_i \in \delta^2 T} S_T \mathbf{n}_T}{\left\| \sum_{T \in \mathcal{T} | x_i \in \delta^2 T} S_T \mathbf{n}_T \right\|} \qquad (3.15) $$

It can be noted that n_c(x_i) ≡ n_N_1(x_i) defined in Eq. (3.8).

— In the works by DU et al. [26], the authors utilize the formulation proposed by MAX [73] :

$$ \mathbf{n}_m(x_i) = \frac{\sum_{T \in \mathcal{T}_i^*} \frac{(\mathbf{x}_j - \mathbf{x}_i) \times (\mathbf{x}_k - \mathbf{x}_i)}{\|\mathbf{x}_j - \mathbf{x}_i\|^2 \|\mathbf{x}_k - \mathbf{x}_i\|^2}}{\left\| \sum_{T \in \mathcal{T}_i^*} \frac{(\mathbf{x}_j - \mathbf{x}_i) \times (\mathbf{x}_k - \mathbf{x}_i)}{\|\mathbf{x}_j - \mathbf{x}_i\|^2 \|\mathbf{x}_k - \mathbf{x}_i\|^2} \right\|} \qquad (3.16) $$

with T_i^* = {T ∈ T | δ^2T = {x_i, x_j, x_k}, (x_j, x_k) ∈ X_2(T)}

— Another procedure, proposed in the paper by ZINCHENKO et al. [138], is based on an iterative calculation of the unit normal vector n_iter(x_i) in order that this vector corresponds exactly to the unit normal vector of the reconstructed interface at x_i.

To distinguish the different approaches used in the Surface Reconstruction methods, from now on, they will be called SR-N_k-n_α-β, with k ∈ {1, 2} for the stencil size, α ∈ {c, m, iter} for the unit normal vector defining the local tangent plane and β ∈ {0, a_00} to enforce or not x_i to belong to the reconstructed surface. It is worth to point out that the SR-N_1-n_m-0 and SR-N_1-n_iter-0 methods were respectively proposed by DU et al. [26] and ZINCHENKO et al. [138].

The calculations of the curvature κ and unit normal vector n are performed on the reconstructed analytic function Z = f(X, Y). The surface tension at node x_i then reads :

$$ \mathbf{f}_i = \sigma\kappa(x_i)\mathbf{n}(x_i) \qquad (3.17) $$

Considering the formulae from differential geometry, the unit normal vector and the mean curvature write [45] :

$$ \tilde{\mathbf{n}} = \frac{\nabla F}{\|\nabla F\|} \qquad (3.18) $$

and

$$ \tilde{\kappa} = -\nabla \cdot \tilde{\mathbf{n}} = \frac{\nabla F \cdot H(F) \cdot \nabla F^T - \|\nabla F\|^2 \mathrm{Trace}(H(F))}{2 \|\nabla F\|^3} \qquad (3.19) $$

where H(F) is the Hessian matrix of function F(X, Y, Z) = Z - f(X, Y), with components H_αβ(F) = ∂²F / ∂α∂β and (α, β) ∈ {X, Y, Z}².

### 3.2.5 Discrete surfaces and solution accuracy

The discrete surfaces are constructed into successive steps. First, a planar surface is covered by a set of equilateral triangles of edge size d. Let us note x̃ any vertex of this mesh, with
coordinates $\tilde{\mathbf{x}} = (x, y, 0)$. Since the mesh of any interface, for example issued from a front tracking method, is often made up of more or less irregular triangles, a disturbed mesh is build from the regular one by using a mapping $\mathcal{P}(\tilde{x}) = \tilde{\tilde{x}}$ such that

$$ \tilde{\tilde{\mathbf{x}}} = \tilde{\mathbf{x}} + r \times p \times d \times (\cos 2\pi\theta, \sin 2\pi\theta, 0) \quad (3.20) $$

In this expression, $(r, \theta)$ is a couple of random variables drawn from a $[0, 1]$-uniform distribution and $p$ is a prescribed number driving the magnitude of perturbations ($0 \le p \le 0.2$). Moreover for steady problems (static case), the vertices are then projected on the known analytical surface $z = h(x, y)$ to define the discrete counterpart of the continuous surface. It is worth pointing out that, even if $p = 0$, the discrete mesh of the surface is not any more made up of equilateral triangles, since the projection has distorted the shape of all the elements. For dynamics problems, the vertices are simply advected in time by the velocity field. When $p \ne 0$, the error on the curvature and normal vector are evaluated on 10 random draws $(r, \theta)$, both for the static and dynamic cases; the retained value is simply the average.

In order to quantify the accuracy of the numerical methods with respect to the exact solutions, two measures of the error are used. Let $\phi$ be either the scalar variable $\kappa$, the unit normal vector $\mathbf{n}$ or the surface tension $\mathbf{f}$. We denote by $\phi(x_i)$ the exact value of $\phi$ at node $x_i$, and $\phi_i$ its numerical approximation. Then, two relative errors are defined over the $N$ vertices of the mesh by:

$$ \text{Err}_2^{\text{rel}}(\phi) = \left( \frac{\sum_{j=1}^{N} (\phi_j - \phi(x_j)) \cdot (\phi_j - \phi(x_j))}{\sum_{j=1}^{N} \phi(x_j) \cdot \phi(x_j)} \right)^{1/2} \quad (3.21) $$

$$ \text{Err}_{\infty}^{\text{rel}}(\phi) = \left( \frac{\max_{j \in [1,N]} (\phi_j - \phi(x_j)) \cdot (\phi_j - \phi(x_j))}{\max_{j \in [1,N]} \phi(x_j) \cdot \phi(x_j)} \right)^{1/2} \quad (3.22) $$

where ‘·’ denotes either the usual product or the scalar product if $\phi$ is a vector. Since the studied surfaces are open, the error calculations are performed only for triangles far enough from the boundaries. Notice that this condition depends on numerical methods because they may rely on different stencil.

The numerous methods presented in this paper are mainly compared through curves which present the relative errors (3.21) (3.22) as a function of the dimensionless initial mesh size defined as $d \times \max_{j \in [1,N]}(\kappa(x_j))$, with $\kappa(x)$ the exact curvature at node $x$. This latter relation can be re-interpreted as a function of the mean curvature radius $1/\kappa(x_j)$: let us assume that $d \times \max_{j \in [1,N]}(\kappa(x_j)) = 10^{-2}$, then the edges of the initial triangles are 100 times smaller than the minimum of mean radius of curvature, $\min_{j \in [1,N]}(1/\kappa(x_j))$. To distinguish in the graphs the two norms (3.21) (3.22), the curves drawn using $L_2$ are continuous whereas those with $L_{\infty}$-norm are discontinuous.

For the sake of brevity, the error curves for the surface tension will not be presented since their behaviour are similar to that of the curvature (same order of convergence and error level).

## 3.3 Analysis of a stationary surface

In this section, the static surface proposed by [67] is studied (Fig. 3.5):
<figure><img src="image_39.png" /><figcaption>FIGURE 3.5 – Top (left) and side (right) views of the discrete mesh projected onto the surface (3.23).</figcaption></figure>

$$h(x, y) = \sin(5x) \sin(5y) \quad \text{with} \quad (x, y) \in [-\frac{\pi}{5}, \frac{\pi}{5}] \times [-\frac{\pi}{5}, \frac{\pi}{5}] \qquad (3.23)$$

The main interest of this analytical surface is that it is characterized by a large variety of curvatures, sum of the two principal curvatures. Points A, B, C and D have maximal curvatures and point O exhibits a zero mean curvature with maximal but opposite principal curvatures (saddle-node point). To avoid boundary effects, the initial flat surface is defined on a larger domain, $(x, y) \in [-2, 2] \times [-2, 2]$ and all the vertices in the sub-domain $(x, y) \in [-\pi/5, \pi/5]^2$ are taken into account for the error evaluation (3.21) (3.22).

### 3.3.1 Results

#### LBO-methods

The Laplace-Beltrami Operator discretization, LBO-Essad. [31] and LBO-Meyer [76] methods, are compared in this section as a function of the mesh size, and for a regular mesh ($p = 0$) and a disturbed one ($p = 0.2$).

Without perturbation ($p = 0$), the LBO-Meyer method provides clearly the best results for the curvature (Fig. 3.6a) : it is second order accurate in $L_2$-norm and tends to be first order

<figure><img src="image_40.png" /><figcaption>FIGURE 3.6 – LBO methods - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$</figcaption></figure>