# Simulation numérique d'ondes interfaciales en milieu océanique

Julien Touboul, Malek Abid, Christian Kharif

► To cite this version:

Julien Touboul, Malek Abid, Christian Kharif. Simulation numérique d'ondes interfaciales en milieu océanique. CFM 2007 - 18ème Congrès Français de Mécanique, Aug 2007, Grenoble, France. hal-03361740

HAL Id: hal-03361740

https://hal.science/hal-03361740v1

Submitted on 1 Oct 2021

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 CC BY 4.0 - Attribution - International License
# Simulation numérique d'ondes interfaciales en milieu océanique

Julien Touboul, Malek Abid & Christian Kharif

IRPHE
49 rue F. Joliot Curie,
Technopôle de Château-Gombert, B.P. 146,
Marseille cedex 13.
Julien.Touboul@irphe.univ-mrs.fr

## Résumé :

Cette étude concerne la simulation numérique directe d'ondes interfaciales séparant deux fluides incompressibles. Le modèle utilisé repose sur l'équation de Navier-Stokes et un suivi d'interface de type Volume Of Fluid (VOF). La résolution de cette équation nécessite d'inverser une équation de type Poisson, et pour cela, on a recours à une méthode multigrille. L'une des principales limitations de cette méthode est le rapport des masses volumiques des deux fluides. Nous nous proposons donc d'étudier l'influence de ce paramètre sur le comportement de l'algorithme. Le problème d'une onde interfaciale stationnaire est alors considéré pour différentes valeurs du rapport des masses volumiques.

## Abstract :

This study deals with the direct numerical simulation of interfacial waves between two incompressible fluids. The model used is based on the Navier-Stokes equation with a Volume of Fluid (VOF) method for tracking the interface's motion. The solution of this equation needs to solve a Poisson like equation with a multigrid algorithm. One of the main limitations of this method is the ratio of density of both fluids. We study herein the influence of this parameter on the behavior of the algorithm. To achieve that goal, the problem of a standing interfacial wave is considered for several values of the density ratio.

## Mots-clefs :

### Ondes interfaciales, algorithme multigrille, rapport de masses volumiques

## 1 Introduction

La connaissance des ondes interfaciales en milieu océanique constitue un enjeu majeur scientifique, notamment en ingénierie côtière. Il s'agit d'un domaine vaste s'étendant des ondes internes aux ondes de surface et leur interaction avec le vent. Parmi ces phénomènes, les vagues scélérates font référence à des vagues géantes qui apparaissent soudain à la surface de la mer. Les témoignages rapportés par Mallory (1974), Lawton (2001) et d'autres, permettent de prendre en compte la mesure de l'enjeu que ces vagues représentent. Des travaux récents ont montré que le vent pouvait jouer un rôle majeur dans la dynamique de telles vagues. En effet, Touboul et al. (2006) et Touboul & Kharif (2006) ont illustré l'importance de ce couplage dans le phénomène.

Cependant, la simulation numérique directe de tels phénomènes présente de nombreuses difficultés de part la diversité des échelles à modéliser. Pour cette raison, de nombreuses méthodes d'inversion de l'équation de type Poisson ont été développées au cours des dernières années (Press et al. (1999)). Notamment, parmi elles, l'algorithme multigrille s'est révélé l'un des meilleurs compromis en terme de vitesse de convergence pour des problèmes multi échelles. L'une des limitations principales de ces méthodes est le rapport de masse volumique des deux
fluides. L'évolution de la convergence de ces méthodes en fonction de ce paramètre paraît donc être un enjeu essentiel de la modélisation du processus physique qui nous intéresse ici.

Le travail proposé s'inscrit donc dans cette démarche. En effet, dans un premier temps, la méthode utilisée est brièvement décrite, puis son comportement global est étudié dans une seconde partie. Pour cela, nous avons recours à des ondes interfaciales stationnaires dans une configuration de stratification stable, qui sont des ondes présentant un grand cisaillement à l'interface, et qui, par conséquent, sont un cas limitant du problème qui nous intéresse.

## 2 Méthode numérique

### 2.1 Formulation mathématique

La méthode présentée dans cette partie s'attache à résoudre les équations du mouvement de deux fluides visqueux, incompressibles, non miscibles, séparés par une interface. Ce problème est un problème particulièrement difficile, dans la mesure où la position de l'interface est l'une des inconnues du problème. Le mouvement des fluides est régi par l'équation de Navier-Stokes

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = - \nabla p + \nabla \cdot (2\mu \mathbf{D}) + \rho g \quad (1) $$

et la condition d'incompressibilité

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

dans lesquelles $\rho$ et $\mu$ désignent la masse volumique et la viscosité dynamique du fluide, et $g$ la force de pesanteur. D'autre part, $\mathbf{u}$ et $p$ désignent respectivement le vecteur vitesse et le champ de pression, tandis que $\mathbf{D}$ est le tenseur des taux de déformation.

Les équations (1) et (2) sont valables dans tout le domaine $\Omega$ composé des deux fluides. $\rho$ et $\mu$ sont alors des fonctions d'espace, dépendant de l'appartenance à l'un ou l'autre des fluides. En introduisant la fonction $\chi$, qui prend pour valeur 0 dans l'un des deux fluides, et 1 dans l'autre, on a

$$ \begin{aligned} \rho &= \chi \rho_1 + (1 - \chi) \rho_2 \\ \mu &= \chi \mu_1 + (1 - \chi) \mu_2 \end{aligned} \quad (3) $$

ou $\rho_1$ et $\rho_2$ désignent les masses volumiques de chacuns des fluides, et $\mu_1$ et $\mu_2$ leurs viscosités. La fonction $\chi$ est alors transportée comme un scalaire passif par l'écoulement, et $\chi$ est donc solution de l'équation

$$ \frac{\partial \chi}{\partial t} + \mathbf{u} \cdot \nabla \chi = 0. \quad (4) $$

la résolution de l'équation (4) permet donc de simuler l'évolution de la fonction $\chi$ dans tout le domaine de calcul $\Omega$.

### 2.2 Méthode de projection

Pour résoudre ce problème, l'on a recours à une méthode de projection. On introduit un champ de vitesses intermédiaire $\mathbf{u}^*$, tel que

$$ \rho \frac{\partial \mathbf{u}^*}{\partial t} = -\rho \mathbf{u}^{(n)} \cdot \nabla \mathbf{u}^{(n)} + \nabla \cdot (2\mu \mathbf{D}^{(n)}) + \rho g \quad (5) $$
ou $\mathbf{u}^{(n)}$ désigne le champ de vitesse au pas de temps $(n)$. Ainsi, le champ de vitesses $\mathbf{u}^{(n+1)}$ au pas de temps $(n+1)$ défini par

$$ \mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\delta t}{\rho} \nabla \mathbf{p} \qquad (6) $$

sera de divergence nulle dès que le champ de pression $p$ sera solution du problème de type Poisson

$$ \left\{ \begin{array}{ll} \nabla \left( \frac{1}{\rho} \nabla p \right) = \frac{1}{\delta t} \nabla \mathbf{u}^* & \text{dans } \Omega \\ \\ \nabla p \cdot \mathbf{n} = - \frac{\rho}{\delta t} (\mathbf{u}^{n+1} - \mathbf{u}^*) \cdot \mathbf{n} = \rho g \cdot \mathbf{n} & \text{sur } \partial \Omega \end{array} \right. \qquad (7) $$

Dans les équations ci dessus, $\delta t$ est le pas de temps, $\Omega$ le domaine de calcul, $\partial \Omega$ sa frontière, et $\mathbf{n}$ la normale à cette frontière.

Le problème constitué des équations (1) et (2) se ramène donc à la résolution du système (7). Ce système est inversé au moyen d'une méthode multigrille similaire à celle présentée par Gueyffier (2000). Cette méthode est choisie pour les avantages numériques qu'elle présente, notamment en terme de temps de calcul. Cependant, une étude minutieuse de l'évolution des taux de convergence en fonction du rapport des masses volumiques $\rho_1/\rho_2$ de telles méthodes n'est pas disponible et s'avère indispensable (voir section 3).

## 2.3 Calcul d'interface affine par morceaux

Il s'agit ici de résoudre l'équation (4). En effet, une fois le champ de vitesse connu, la fonction $\chi$ peut être transportée par l'écoulement. On choisit pour cela une méthode de type "Piecewise Linear Interface Calculation" (PLIC), proposée par Li (1995). Le principe de cette méthode est de reconstruire l'interface comme une fonction affine par morceaux sur le maillage du calcul. La fonction $\chi$, une fois discrétisée, est notée $C$. Elle peut être interprétée comme la fraction volumique de l'un des fluides dans la cellule considérée. On peut ainsi définir la normale à l'interface $\mathbf{n} = \nabla C = (n_x, n_y)^t$ dans chaque cellule. Ainsi, l'interface à pour équation

$$ n_x x + n_y y = \alpha. \qquad (8) $$

$\alpha$ est un paramètre à déterminer. Pour cela, on dispose de la fraction volumique $C$. En effet, l'aire comprise sous l'interface reconstruite doit correspondre à la fraction volumique. On peut ainsi déterminer $\alpha$.

Une fois l'équation de l'interface connue, on peut procéder à une seconde étape pendant laquelle on advecte cette interface grâce à une méthode lagrangienne. En interpolant linéairement le champ de vitesses dans la cellule, l'interface advectée est une droite d'équation

$$ n'_{x}x + n'_{y}y = \alpha' \qquad (9) $$

avec

$$ \begin{aligned} n'_{x} &= n_{x} / (1 - u_{1} \delta t / dx + u_{2} \delta t / dx) \\ n'_{y} &= n_{y} / (1 - v_{1} \delta t / dy + v_{2} \delta t / dy) \\ \alpha' &= \alpha + n'_{x} u_{1} \delta t / dx + v_{1} \delta t / dy \end{aligned} \qquad (10) $$

$u_1$ et $u_2$ désignent respectivement les composantes horizontales de la vitesse à gauche et à droite de la cellule, tandis que $v_1$ et $v_2$ sont les composantes verticales en bas et en haut de la cellule. La méthode PLIC est stable et satisfait la contrainte physique $0 \le C \le 1$ dès que la condition de Courant-Friedrichs-Lewy (CFL) $||\mathbf{u}||_\infty \delta t/min(dx, dy) \le 1/2$ est vérifiée.
# 3 Vérification de la méthode

## 3.1 Tests préliminaires

Un code basé sur la méthode développée par Lafaurie et al. (1994) a été développé pour cette étude. Nous nous attachons dans un premier temps à vérifier le bon fonctionnement des algorithmes mis en œuvre. Tout d'abord, la méthode PLIC décrite dans la section précédente est utilisée dans différents champs de vitesses, notamment en translation et en rotation. La Figure 1 (à gauche) présente l'évolution d'un disque transporté par un champ de rotation $\mathbf{u} = (5 - y, x - 5)^t$. Le maillage est de $256 \times 256$. On constate que ce disque est très bien conservé. Notamment, il apparaît que la masse est conservée à environ 0.1% près au cours de toutes les simulations présentées ici. On peut considérer que l'algorithme PLIC fonctionne de manière très satisfaisante.

<figure><img src="image_2.png" /><figcaption>FIG. 1 – A gauche : évolution d'un disque advecté par la méthode PLIC dans un champ de rotation $\mathbf{u} = (5 - y, x - 5)^t$. A droite : convergence de l'algorithme multigrille vers une fonction de classe $C^\infty$ en fonction du nombre de V-cycles effectués. La courbe du haut est la norme infinie de l'erreur, celle du bas est la norme infinie du résidu.</figcaption></figure>

Nous étudions également le taux de convergence de la méthode multigrille, destinée à résoudre l'équation de type Poisson, vers la fonction de classe $C^\infty$

$$p_{\text{exacte}}(x, y) = e^{-(x-5)^2} \times e^{-(y-5)^2} \quad (11)$$

La Figure 1 (à droite) présente l'évolution de la norme infinie de l'erreur $||p - p_{\text{exacte}}||_\infty$ (courbe du haut), ainsi que la norme infinie du résidu de la méthode multigrille (courbe du bas), en fonction du nombre de V-cycles réalisés par l'algorithme multigrille. Ces courbes sont représentées pour différentes valeurs du paramètre $s = \rho_1/\rho_2$. Il est intéressant de constater que le taux de convergence de cet algorithme n'est absolument pas affecté par ce paramètre dès qu'il converge vers une fonction de classe $C^\infty$. Il est également intéressant de constater qu'une erreur de l'ordre de $10^{-6}$ en terme de résidus correspond tout de même à une erreur absolue de l'ordre de $10^{-2}$. Nous retenons donc cette erreur en terme de résidus comme erreur maximale à obtenir lors de la procédure multigrille.
## 3.2 Simulation d'ondes interfaciales stationnaires

Bien que les résultats de la section précédente soient satisfaisants, il est nécessaire de poursuivre plus avant les tests de convergence de cette méthode, dans la mesure où les problèmes qui seront abordés par la suite n'admettent pas de solutions de classe $C^\infty$. En effet, l'interface entre les deux fluides introduit une discontinuité du gradient de pression, qui résulte en une chute de la convergence de la méthode.

<figure><img src="image_3.png" /><figcaption>FIG. 2 – A gauche : évolution de la relation de dispersion en fonction du paramètre $s$. (–) évolution théorique. (o) résultat obtenu grâce à la méthode numérique. A droite : évolution de la surface libre en fonction du temps pour $s = 20$. L'enveloppe correspond à la décroissance théorique due à la viscosité.</figcaption></figure>

Nous choisissons donc d'étudier l'évolution d'une onde interfaciale stationnaire grâce à cette méthode. La Figure 2 (à gauche) représente la pulsation en fonction du paramètre $s$. La courbe en trait plein est l'évolution théorique prédite par Lamb (1932) et Plesset & Whipple (1974), à savoir, en profondeur infinie,

$$ \omega^2 = \frac{s-1}{s+1} gk. \qquad (12) $$

On constate que la méthode a convergé vers le résultat dans tous les cas correspondant à $1 \le s \le 10^3$. Toutefois, il est important de signaler que le nombre de V-cycles nécessaire à la convergence n'est absolument pas équivalent dans chacun de ces cas. Les temps de calcul ont pu varier d'un facteur 4 pour ces différentes simulations.

La figure 2 (à droite) montre l'évolution temporelle de la surface libre en $x = 0$ pour $s = 20$. L'enveloppe en pointillés représente l'atténuation théorique

$$ \frac{A(t)}{A_0} = e^{-\nu k^2 t} \qquad (13) $$

liée à la présence de viscosité ($\nu = \mu_1/\rho_1 = \mu_2/\rho_2$ représente ici la viscosité cinématique des deux fluides). Il est intéressant de noter, une fois encore, le bon comportement de la méthode. Cependant, la viscosité utilisée est importante ($\nu = 10^{-3}m^2/s$), ce qui a pour effet de stabiliser le schéma. La convergence, d'une manière générale, a tendance à se dégrader lorsque la viscosité chute. Ceci est du à la présence d'un cisaillement très important à l'interface. Il semble donc que les ondes progressives, qui présentent un meilleur comportement de ce point de vue, permettront de simuler des viscosités de l'ordre de celle de l'eau.
## 4 Conclusions

La méthode présentée ici est testée dans différents cas de figure. Notamment, nous nous attachons à déterminer l'évolution de la convergence de cette méthode lorsque le rapport de masses volumiques croît, notamment jusqu'à 1000. Les tests mis en place sont d'une part la convergence vers une fonction infiniment dérivable, et d'autre part, la convergence vers un problème d'ondes interfaciales stationnaires. Bien que la convergence se détériore dans ce dernier cas, nous constatons que le temps de calcul nécessaire à l'obtention de la convergence reste acceptable. Ces résultats sont donc très encourageants quand à la modélisation d'ondes océaniques interfaciales, incluant l'étude de l'interaction entre vent et vagues.

## Références

- Gueyffier D., 2000, Etude de l'impact de gouttes sur un film liquide mince. Développement de la corole et formation de projections. **Thèse de doctorat**, Université Pierre et Marie Curie.

- Jeffreys H., 1925, On the formation of water waves by wind. *Proc. Roy. Soc. London Ser. A*, **107**, 189-206.

- Kharif C. & Pelinovsky E., 2003, Physical mechanisms of the rogue wave phenomenon. *Eur. J. Mech. B Fluids*, **22**, 603-634.

- Lafaurie B., Scardovelli R., Zaleski S. & Zanetti G., 1994, Modelling merging and fragmentation in multiphase flows with SURFER. *J. Comp. Phys.* **113** (1), 134-147.

- Lamb H., 1932, Hydrodynamics. *Dover Publications*, New York.

- Lawton G., 2001, Monsters of the deep (the perfect wave). *New Scientist*, **170**, 28-32.

- Li J., 1995, Calcul d'intresface affine par morceaux. *C. R. Acad. Sci. Paris, Serie II b*, **320**, 391-396.

- Mallory, J.K., 1974, Abnormal waves on the South-East Africa. *Int. Hydrog. Rev.*, **51**, 89-129.

- Plesset, S. M. & Whipple, C. G., 1974, Viscous effect in Rayleigh-Taylor instability. *Phys. Fluids*, **17** (1), 1-7.

- Press, W. H., Flannery, B. P., Teukolsky, S. A. & Vetterling W. T., 1999, Numerical Recipes in Fortran 77 : The art of scientific computing. *Cambridge University Press*.

- Touboul J., Giovanangeli J. P., Kharif C. & Pelinovsky E., 2006, Freak waves under the action of wind : experiments and simulations. *Eur. J. Mech. B Fluids*, **25**, 662-676.

- Touboul J. & Kharif C., 2006, On the interaction of wind and extreme gravity waves due to modulationnal instability. *Phys. Fluids*, **18**, 108103.