1  Modélisation

Ce chapitre décrit les modèles physiques utilisés dans CERF pour simuler des écoulements (sans érosion dans la version V1.0), ainsi que la manière dont ces modèles sont implémentés dans le code. Il est divisé en trois sections. La première section présente dans le cas général la méthode des volumes finis. La deuxième section décrit le modèle d’écoulement en eaux peu profondes de Saint Venant. La troisième section présente le modèle d’écoulement bifluide eulérien utilisé dans CERF.

1.1 La méthode des volumes finis

La méthode des volumes finis est une méthode numérique utilisée pour résoudre des équations aux dérivées partielles. Cette méthode est utilisée dans CERF pour résoudre les équations de Saint Venant et les équations d’écoulement bifluide.

1.1.1 Principe de la méthode des volumes finis

Dans cette section, nous rappelons rapidement l’approximation par volumes finis, en dimension \(n_{dim}\), d’un système d’équations de conservation de la forme : \[ \frac{\partial \overrightarrow{w}}{\partial t} + \nabla \cdot \overline{\overline{F}}(\overrightarrow{w}) = \overrightarrow{S}(\overrightarrow{w}) \tag{1.1}\] où \(\overrightarrow{w}(\overrightarrow{x},t)\) est le vecteur des variables conservatives de taille \(n_{dof}\) dépendant de l’espace et du temps, \(\overline{\overline{F}}\) la matrice des flux de taille \(n_{dof} \times n_{dim}\), \(\overrightarrow{S}\) le vecteur des termes sources de taille \(n_{dof}\) et \(\nabla \cdot\) l’opérateur divergence.

Le système est résolu dans \(\mathbb{R}^{n_{dim}}\). Le domaine de calcul \(\Omega \subset \mathbb{R}^{n_{dim}}\) est divisé en un ensemble de volumes de contrôle, également appelés cellules, \(\Omega = \cup_{k} C_k\). Nous commençons par intégrer l’Equation 1.1 sur une cellule \(C_k\) au moyen de la formule de Green :

\[\int _{C_k} \frac{\partial \overrightarrow{w}}{\partial t} d \Omega + \int _{C_k} \nabla \cdot \overline{\overline{F}} d \Omega = \int _{C_k} \overrightarrow{S}(\overrightarrow{w}) d \Omega\]

\[ \int _{C_k} \frac{\partial \overrightarrow{w}}{\partial t} d \Omega + \sum _{a} \int _{\partial C_{k/a}} \overline{\overline{F}} \cdot \overrightarrow{n}_{k/a} d \partial\Omega = \int _{C_k} \overrightarrow{S}(\overrightarrow{w}) d \Omega \tag{1.2}\]

où \(\overrightarrow{n}_{k/a}\) désigne le vecteur normal unitaire sur la frontière \(\partial C_{k/a}\) entre la cellule \(C_k\) et la cellule voisine \(C_a\). La sommation est effectuée sur toutes les cellules voisines de \(C_k\).

La méthode des volumes finis consiste alors à approcher la solution par des fonctions constantes par cellule, si bien qu’en posant \[\overrightarrow{w}_k(t) \approx \frac {1}{|C_k|}\int _{C_k} \overrightarrow{w}(\overrightarrow{x},t) d \Omega\] où \(|C_k|\) représente le volume de la cellule \(C_k\), nous approchons l’Equation 1.2 par \[ |C_k| \frac{\partial \overrightarrow{w}_k(t)}{\partial t} + \sum _{a} |\partial C_{k/a}| \: \widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a} \right) = |C_k|S_k(t) \tag{1.3}\]

où \(S_k(t)\approx \frac {1}{|C_k|}\int _{C_k}S(\overrightarrow{x},t) d \Omega\) représente la valeur moyenne du vecteur sollicitation sur la cellule, \(|\partial C_{k/a}|\) la surface de la face entre les cellules \(C_k\) et \(C_a\) et \(\widehat{\overrightarrow{F}}\) une approximation du flux.

Comme les inconnues \(\overrightarrow{w}\) ne sont pas définies à l’interface entre les cellules, le flux ne peut être calculé aisément. On l’approche donc par un flux numérique dépendant de l’interface et des valeurs des variables dans les cellules adjacentes. On note :

\[\overrightarrow{F}_{k/a}(t) = \widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a} \right) \approx \frac{1}{|\partial C_{k/a}|}\int _{\partial C_{k/a}} \overline{\overline{F}} \cdot \overrightarrow{n}_{k/a} ds\]

Il existe de très nombreuses approches : Rusanov, HLLC, … Dans CERF le flux numérique, par la méthode de Godunov, est calculé en utilisant la solution exacte du problème de Riemann à l’interface \(k/a\) (pour plus de détails, voir Section 1.2.5 et Section 1.3.5).

On peut atteindre l’ordre 2 en espace en utilisant la méthode MUSCL (Monotone Upstream Scheme for Conservation Laws) pour reconstruire les variables primitives aux interfaces. À partir des valeurs des variables sur les cellules voisines, on estime le gradient des variables primitives. Soit \({w_p}_i\) la ième composante des variables primitives et \({w_p}_i^k\) sa valeur au centre de la cellule, on peut estimer le gradient sur la cellule \(C_k\) comme : \[ |C_k| \widetilde{\nabla} {w_p}_i = \int _{C_k} \nabla {w_p}_i d\Omega =\sum _{a} \int _{\partial C_{k/a}} w_{p_i} \overrightarrow{n}_{k/a} d\partial\Omega =\sum_{a} |\partial C_{k/a}| \frac{1}{2}\left( {w_p}_i^k + {w_p}_i^a \right) \overrightarrow{n}_{k/a}\] et on peut alors reconstruire les variables primitives de la cellule \(C_k\) de centre \(K\) sur les centres \(F\) des faces de la cellule telles que : \[w_{p_i}^F = w_{p_i}^k + \widetilde\nabla w_{p_i} \cdot \overrightarrow{KF}\]

Pour atteindre l’ordre 2, on remplace alors \(\widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a} \right)\) par \(\widehat{\overrightarrow{F}} \left(\overrightarrow{w}_k^F(t),\overrightarrow{w}_a^F(t),\overrightarrow{n}_{k/a} \right)\)

Note

En général il est nécessaire d’adjoindre un limiteur à la méthode MUSCL. Dans CERF le limiteur de Barth-Jespersen est utilisé : \(Min(w_{p_F},w_{p_F}) \leq w_{p_i}^F \leq Max(w_{p_F},w_{p_F})\).

1.1.2 Discrétisation temporelle

Afin de faciliter la parallélisation du code, CERF ne dispose que de schémas explicites: le schéma de Runge-Kutta d’ordre 1 ou 2. L’intervalle de temps \(\left[0,T\right]\) est discrétisé de telle façon que \(t_{n+1} = t_n + \delta t_n\) où \(\delta t_n\) représente le pas de temps au temps \(t_n\). La méthode de Runge-Kutta d’ordre 1, ou encore méthode d’Euler explicite d’ordre 1, consiste donc à réécrire l’Equation 1.3 à l’instant \(t_n\) en approchant la dérivée temporelle comme une différence finie, de telle sorte que: \[ |C_k|\frac{ \overrightarrow{w}_k(t_{n+1}) - \overrightarrow{w}_k(t_{n})}{\delta t_n} + \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}(t_n) = |C_k|\overrightarrow{S}_k(t_n)\]

En notant \(\overrightarrow{w}_k(t_{n})=\overrightarrow{w}_k^{n}\), on calcule donc de façon explicite la nouvelle valeur des variables conservatives à l’instant \(t_{n+1}\): \[ \overrightarrow{w}_k^{n+1} = \overrightarrow{w}_k^{n}+\frac{\delta t_n}{|C_k|} \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}^n - \delta t_n \overrightarrow{S}_k^n\]

La méthode de Runge-Kutta d’ordre 2 consiste à calculer la solution à l’instant \(t_{n+1}\) en deux étapes. La première étape consiste à calculer la solution à l’instant \(t_{n+\frac{1}{2}}\) en utilisant le schéma d’Euler explicite d’ordre 1. La deuxième étape consiste à calculer la solution à l’instant \(t_{n+1}\) en utilisant le schéma d’Euler explicite d’ordre 1 avec la solution obtenue à l’instant \(t_{n+\frac{1}{2}}\). \[ \overrightarrow{w}_k^{n+\frac{1}{2}} = \overrightarrow{w}_k^{n}+\frac{\delta t_n}{2 |C_k|} \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}^n - \frac{\delta t_n}{2} \overrightarrow{S}_k^n \] \[ \overrightarrow{w}_k^{n+1} = \overrightarrow{w}_k^{n}+\frac{\delta t_n}{|C_k|} \sum _{a} |\partial C_{k/a}| \: \overrightarrow{F}_{k/a}^{n+\frac{1}{2}} - \delta t_n\overrightarrow{S}_k^{n+\frac{1}{2}} \]

Important 1.1: Condition de stabilité de Courant-Friedrichs-Lewy

Les méthodes explicites sont simples à mettre en oeuvre mais elles sont conditionnées par la stabilité de Courant-Friedrichs-Lewy (CFL). La condition CFL est une condition nécessaire pour garantir la stabilité du schéma numérique. Cette condition est dépendante de la vitesse de propagation des ondes dans le milieu et de la taille des cellules. En pratique, elle est déterminée pour chaque modèle hyperbolique traité. Dans CERF la taille des cellules \(d_k\) est déterminée par le ratio entre le volume de la cellule et la surface de sa frontière: \[ d_k = \frac{|C_k|}{\sum_{a}|\partial C_{k/a}|}\]
Dans le cas 1D ou 2D, les faces de symétrie ne sont pas considérées dans le calcul de la surface des frontières. La condition CFL est alors donnée par : \[ \delta t_n \leq \alpha_{CFL} \frac{\min_{k}d_k}{\max_{k} \lambda_k}\]
où \(\alpha_{CFL}\) est un coefficient à définir par l’utilisateur et \(\lambda_k\) est la vitesse de propagation des ondes dans la cellule \(C_k\) dépendant du modèle utilisé (Note 1.1 pour Saint-Venant et Note 1.2 pour l’écoulement bifluide).

1.1.3 Discrétisation spatiale

Important

Dans sa version V1.0 CERF dispose exclusivement de cellules hexaédriques et donc de faces quadrangulaires ! La généralisation est en cours de développement.

Note

CERF permet l’utilisation de maillages dynamiques, c’est à dire que régulièrement il est possible de raffiner ou déraffiner le maillage (AMR : Adaptive Mesh Refinement).

CERF est un code nativement tridimensionnel. Dans sa version initiale, CERF ne dispose que de cellules hexaédriques et donc de faces quadrangulaires. Les cellules et les faces ne sont pas nécessairement régulières. Pour autant, comme illustré Figure 1.1, il est possible d’effectuer des simulations bidimensionnelles en réalisant un maillage avec une seule cellule dans la direction \(z\) et en imposant des conditions de type “miroir” sur les faces dans le plan \(x,y\). De même, il est possible de réaliser des simulations unidimensionnelles en réalisant un maillage avec une seule cellule dans les directions \(y\) et \(z\) et en imposant des conditions de type “miroir” sur les faces dans le plan \(x,z\) et \(x,y\).

Figure 1.1: Maillages tri-bi-unidimensionnel avec des cellules hexaédriques

1.1.3.1 Stratégie de maillage

CERF permet d’utiliser des outils de raffinement et de déraffinement de maillage (AMR: Adaptive Mesh Refinement). Afin de ne pas pénaliser le temps de calcul, la modification du maillage n’est pas réalisée à chaque pas de temps, mais à des instants déterminés par l’utilisateur. Cette approche permet de conserver une bonne qualité de maillage tout en limitant le coût de calcul. Cette méthode est dénommée BB-AMR (Block Based Adaptive Mesh Refinement). Nous allons décrire ci-après les différentes étapes de la stratégie de maillage utilisée dans CERF.

Figure 1.2: Maillage initial conforme : “Master Mesh”.

Étape 1 : Pour simplifier l’explication le maillage présenté Figure 1.2 est constitué de quadrangles réguliers, qui peuvent représenter un maillage hexaédrique vu suivant \(z\) avec une seule cellule suivant \(z\). Comme illustré Figure 1.2, on part d’un maillage initial conforme constitué uniquement d’hexaèdres de forme quelconque. Dans un cas simple le maillage peut être créé par CERF, sinon il peut être importé au format GMSH V2.2. Ce maillage est dénommé le maillage maître ou Master Mesh. Il est primordial que le maillage soit conforme afin de définir aisément les faces reliant les cellules. on appellera bloc les cellules de ce maillage initial.

Figure 1.3: Numérotation du maillage initial suivant les codes de Morton.

Étape 2: On ordonne les blocs du maillage initial en les numérotant suivant le Z-order ou Code de Morton (Figure 1.3). Par cette approche, que l’on appelle Space Filling Curve (Section 1.1.3.2), la numérotation permet de parcourir les blocs, même en trois dimensions, en suivant une numérotation telle que la plupart du temps des blocs successifs dans la numérotation soient voisins dans l’espace.

Figure 1.4: Définition du niveau de raffinement des blocs.

Étape 3 : On définit ensuite le niveau de raffinement des blocs, comme illustré Figure 1.4. Initialement, le niveau de raffinement est défini par l’utilisateur. En cours de calcul, le niveau de raffinement est défini par un critère Section 1.1.3.3. Comme vérifié dans Ersoy, Golay, and Yushchenko (2013) ou Thomas Altazin and Yushchenko (2016) la différence de niveau de raffinement entre deux blocs ne doit pas excéder \(1\). Lors du raffinement le bloc est divisé par \(2\) dans chaque direction et ce pour chaque niveau nécessaire. Ainsi, si on note \(Niv\) le niveau de raffinement requis (où \(Niv=0\) représente un bloc qui n’est pas raffiné) et \(dim=1,2,3\) la dimension du problème alors le nombre de cellules dans un bloc sera : \(2^{dim \times Niv}\). Dans l’exemple Figure 1.4, le bloc \(1\) n’est pas raffiné, le bloc \(2\) est raffiné de \(1\) niveau et le bloc \(4\) de \(2\) niveaux. Ce qui conduit alors, comme illustré Figure 1.5 dans le cas où \(n_{dim}=2\), à un nombre de \(1\) cellule dans le bloc \(1\), \(4\) cellules dans le bloc \(2\) et \(16\) cellules dans le bloc \(4\).

Figure 1.5: Décomposition du maillage en domaines.

Étape 4 : Comme les schémas d’intégration en temps dans CERF sont explicites, il est aisé de paralléliser l’algorithme de résolution. Pour cela, il est nécessaire de découper le maillage en plusieurs domaines, si possible avec un nombre de cellules équilibré et des interfaces minimales afin de minimiser les échanges de données entre les domaines. Il existe de nombreux algorithmes, parfois complexes, qui effectuent cette opération de façon performante. Dans CERF nous avons choisi un compromis entre simplicité et efficacité. Nous avons opté pour un découpage piloté par les codes de Morton. Connaissant le nombre de cellules par bloc, on parcourt les blocs en suivant la numérotation de Morton et on les attribue à un domaine. On s’assure que le nombre de cellules par domaine est équilibré. Comme illustré Figure 1.5, le maillage initial est découpé en 3 domaines composés de 25, 21 et 18 cellules. Chaque domaine de calcul sera ensuite affecté à un processus de calcul

Figure 1.6: Réalisation du maillage.

Étape 5 : La dernière étape consiste alors à réaliser le maillage. Afin de communiquer les informations entre les domaines, CERF créé des cellules fantômes. Ces cellules sont des copies des cellules voisines du domaine. Elles sont utilisées pour calculer les flux numériques aux interfaces entre les domaines. Comme illustré Figure 1.6, les cellules fantômes sont représentées en gris. On notera que le maillage ainsi réalisé n’est pas conforme! Par exemple, la première cellule en bas à gauche du domaine \(1\) de Figure 1.6 est entourée de 6 faces (dans le plan \(x,y\)). Le raffinement des blocs définissant le maillage est bien sûr fonction de la dimension du problème à traiter (Figure 1.7). Enfin, si le maillage est modifié en cours de calcul, il est nécessaire d’effectuer une projection des variables conservatives sur le nouveau maillage. Dans une opération de remaillage on ne peut diminuer ou augmenter d’un niveau supérieur à \(1\), si bien que l’étape de projection est relativement simple.

Figure 1.7: Raffinement en fonction de la dimension du problème.

1.1.3.2 Code de Morton

Les courbes de remplissages ou Space Filling Curves sont des courbes qui permettent de parcourir un espace de dimension supérieure. Courbes de Peano, puis de Hilbert ou Lebesgue elles sont utilisées par exemple pour indexer des données et en particulier dans la construction de quadtree ou octree. Dans CERF, la numérotation par code de Morton ou Z-oder curve est utilisée. Cette numérotation permet d’indexer un ensemble de données (points, cellules, faces, ….) de telle sorte que des données voisines dans l’espace soient voisines dans l’indexation. Cette propriété est utilisée dans CERF pour découper le maillage en domaines de calcul en indexant les blocs du maillage maître. Considérons un point de coordonnées \((x,y,z)\) dans un espace tridimensionnel. Pour indexer ce point, on commence par convertir les coordonnées en entiers positifs. Pour ce faire on définit les coordonnées les plus petites du domaine d’étude \((x_{min},y_{min},z_{min})\) et une distance \(d_{morton}\) définissant la grille du domaine, par exemple un dixième de la taille de la plus petite cellule du domaine. On définit alors les coordonnées entières du point \((ix,iy,iz)\) par: \[ix = \frac{x-x_{min}}{d_{morton}} ,\: iy = \frac{y-y_{min}}{d_{morton}} ,\: iz = \frac{z-z_{min}}{d_{morton}} \] La numérotation de Morton consiste alors à convertir ces coordonnées entières en base \(2\) \((ix_2,iy_2,iz_2)\), puis à entrelacer les bits de chaque coordonnées afin de former un seul nombre entier en base \(2\) \(M_2\). Enfin, revenant en base \(10\), on obtient l’index de Morton \(M\) du point \((x,y,z)\).

Example 1.1 \(ix=2\), \(iy=5\), \(iz=0\).
On a alors \(ix_2=010\), \(iy_2=101\) et \(iz_2=000\).
En entrelaçant les bits on obtient \(M_2=010100010\) et en base \(10\) \(M=162\).

Note

Tous les outils d’indexation de CERF se trouvent dans le module ~/sources/UTI/mod_zorder.f90.

1.1.3.3 Critère de raffinement

En général, un critère de raffinement est défini pour chaque modèle physique traité. Ce critère est souvent issu de considérations physiques à partir, par exemple, du gradient d’une des variables primitives ou de façon plus pertinente et plus complexe d’un calcul d’erreur a posteriori ou de la construction d’un indicateur d’erreur. Dans le cas des problèmes hyperboliques, il existe un critère plus général issu des travaux de Croisille (1990). Ce critère est basé sur la production numérique d’entropie. Il a été utilisé par Golay (2009), Ersoy, Golay, and Yushchenko (2013) et longuement étudié par Puppo (2004). L’inégalité d’entropie de Lax (1971), nous permet de dire que pour tout système hyperbolique de lois de conservation de la forme Equation 1.1, il existe une entropie mathématique \(s(\overrightarrow{w})\) et un flux d’entropie \(\overrightarrow{\psi}_s(\overrightarrow{w})\) tels que : \[\frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \leq 0 \tag{1.4}\] où le flux d’entropie \(\overrightarrow{\psi}_s(\overrightarrow{w})\) doit satisfaire une équation de compatibilité \(\nabla_w \psi_s = \nabla_w \overline{\overline{F}} \cdot \nabla_w s\). Pour des solutions régulières l’Equation 1.4 est une égalité. Dans le cas d’un choc, c’est une inégalité stricte. Or, lors d’une simulation numérique, on constate que la production numérique d’entropie n’est pas nulle! Il s’avère alors que la production numérique d’entropie est un très bon indicateur d’erreur, à l’exception des chocs où elle tend vers l’infini. On définit donc \(p^s_k\) la production numérique d’entropie sur la cellule \(C_k\) par :

\[p^s_k = \frac{1}{|C_k|} \int_{C_k} \left( \frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \right) d\Omega\].

Le critère de raffinement \(C_b\) est défini par bloc comme le ratio de la production numérique d’entropie dans le bloc par rapport à la production numérique d’entropie totale sur le domaine. \[ C_b = \frac{\int_{Bloc} \left( \frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \right) d\Omega}{\int_{\Omega} \left( \frac{\partial s(\overrightarrow{w})}{\partial t} + \nabla \cdot \overrightarrow{\psi}_s(\overrightarrow{w}) \right) d\Omega}\] Dans CERF, le critère de raffinement est défini par bloc. Si le critère est supérieur à une valeur seuil \(\alpha_{max}\), le bloc est raffiné. Si le critère est inférieur à une valeur seuil \(\alpha_{min}\), le bloc est déraffiné. Les valeurs seuils sont définies par l’utilisateur.

1.2 Modèle de Saint Venant

1.2.1 Équations constitutives

Le modèle d’écoulement en eaux peu profondes ou modèle de Saint-Venant est un modèle d’écoulement simplifié. Il est dérivé des équations de Navier-Stokes en négligeant les termes de viscosité, assumant une répartition hydrostatique de la pression et en supposant que la hauteur de l’eau est beaucoup plus petite que la longueur caractéristique de l’écoulement. Le modèle original de Barré de Saint-Venant est un modèle 1D Saint-Venant (1871), la dérivation du modèle est expliqué par exemple par Marche (2007) ou Poussel (2024). Les développements du modèle de CERF ont été initiés par Pons (2018). Dans un repère \((O,x,y,z)\), on considère un écoulement en eaux peu profondes sur un domaine où la bathymétrie est définie par \(z_b(x,y)\) avec une hauteur d’eau \(h(x,y)\), si bien que le niveau de la surface libre est défini par \(\zeta=h+z_b\), comme illustré Figure 1.8.

Figure 1.8: Notations du modèle d’écoulement en eaux peu profondes

Les équations de Saint-Venant sont les suivantes : \[ \frac{\partial h}{\partial t} + \nabla \cdot (h \overrightarrow{u}) = 0 \tag{1.5}\] \[ \frac{\partial h\overrightarrow{u}}{\partial t} + \nabla \cdot (h \overrightarrow{u} \otimes \overrightarrow{u} + g\frac{h^2}{2} \overline{\overline{I}} ) = -gh \nabla z_b -gh\overrightarrow{S_f} \tag{1.6}\]

où \(\overrightarrow{u}=<u,v>^T\) représente la vitesse de l’écoulement, \(g\) l’accélération de la pesanteur, \(S_f\) le terme de frottement et où l’on note \(\overline{\overline{I}}\) la matrice identité, \(\nabla\) l’opérateur gradient, \(\nabla \cdot\) l’opérateur divergence et \(\otimes\) le produit tensoriel. On définit les variables conservatives \(\overrightarrow{w}=<h,hu,hv>^T\) et les variables primitives \(\overrightarrow{w_p}=<h,u,v>^T\). On retrouve la forme générale l’Equation 1.1 avec : \[\overline{\overline{F}}(\overrightarrow{w}) = \left[ \begin{array}{cc} hu & hv \\ hu^2 + g\frac{h^2}{2} & h u v \\ h u v & hv^2 + g\frac{h^2}{2} \end{array}\right]\]

et
\[\overrightarrow{S} = -gh \nabla z_b -gh\overrightarrow{S_f} \]

1.2.2 Prise en compte de la bathymétrie

La prise en compte de la bathymétrie dans le schéma de résolution par volumes finis nécessite une attention particulière. En effet, le saut de bathymétrie inhérent au fait que la bathymétrie soit considérée constante par cellule, peut générer des vagues parasites. La méthode de reconstruction, proposée par Audusse et al. (2004) est donc implémentée dans CERF. Soit \(\Delta z_{k/a}\) le saut de bathymétrie à l’interface entre les cellules. À l’interface, on distingue alors le flux \(\overrightarrow{F}_k\) contribuant à la cellule \(C_k\) et le flux \(\overrightarrow{F}_a\) contribuant à la cellule \(C_a\), tels que:

\[\overrightarrow{F}_k \left( \overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a}, \Delta z_{k/a}\right)=\overrightarrow{F} \left(\overrightarrow{w}_k^*(t),\overrightarrow{w}_a^*(t),\overrightarrow{n}_{k/a} \right) + \left\{ \begin{array}{cc} 0 \\ \frac{g}{2} \left(h_k^2 - {h_k^*}^2 \right) \overrightarrow{n}_{k/a} \end{array} \right\} \] \[\overrightarrow{F}_a \left( \overrightarrow{w}_k(t),\overrightarrow{w}_a(t),\overrightarrow{n}_{k/a}, \Delta z_{k/a}\right)=\overrightarrow{F} \left(\overrightarrow{w}_k^*(t),\overrightarrow{w}_a^*(t),\overrightarrow{n}_{k/a} \right) + \left\{ \begin{array}{cc} 0 \\ \frac{g}{2} \left(h_a^2 - {h_a^*}^2 \right) \overrightarrow{n}_{k/a} \end{array} \right\} \] où \(\overrightarrow{w}_k^*=<h_k^*,\overrightarrow{u}_k>^T\) et \(\overrightarrow{w}_a^*=<h_a^*,\overrightarrow{u}_a>^T\), avec : \[h_k^* = max \left(0, h_k - max \left( 0, \Delta z_{k/a}\right) \right)\] \[h_a^* = max \left(0, h_a - max \left( 0, \Delta z_{k/a}\right) \right)\]

1.2.3 Critère de raffinement

Dans le cas du modèle de Saint-Venant, pour définir la production numérique d’entropie Equation 1.4, on utilise l’entropie définie par : \[ s(\overrightarrow{w}) = \frac{1}{2} \left( h ||\overrightarrow{u}||^2 + gh^2 + 2 g h z_b \right) \] et le flux d’entropie \[ \overrightarrow{\psi}_s(\overrightarrow{w}) = \left( s + \frac{g h^2}{2} \right) \overrightarrow{u} .\]

1.2.4 Prise en compte de la friction

Dans l’équation de conservation de la quantité de mouvement Equation 1.6, le terme de frottement est pris en compte. Dans CERF, le terme de frottement est modélisé par le terme de Manning-Strickler ou le terme de Darcy-Weisbach. Le terme de Manning-Strickler est défini par : \[ \overrightarrow{S_f} = C_f \frac{||\overrightarrow{u}|| \overrightarrow{u}}{h^{4/3}}\] où \(C_f = n^2\) est le coefficient de frottement de Manning-Strickler.

Le terme de Darcy-Weisbach est défini par : \[ \overrightarrow{S_f} = C_f \frac{||\overrightarrow{u}|| \overrightarrow{u}}{h}\]

Le terme de friction est traité par une méthode de “splitting”. On résout d’abord le problème Equation 1.1 sans le terme de friction, \[ \frac{\partial \overrightarrow{w}}{\partial t} + \nabla \cdot \overline{\overline{F}}(\overrightarrow{w}) = -gh \nabla z_b \tag{1.7}\]
puis on résout \[ \frac{\partial \overrightarrow{w}}{\partial t} = -gh\overrightarrow{S_f} \tag{1.8}\]

On fait alors l’hypothèse que la friction ne modifie ni la hauteur d’eau, ni la direction de l’écoulement, mais simplement l’intensité de la vitesse de l’écoulement. On peut alors résoudre le problème de friction sur la vitesse par un schéma implicite d’ordre \(1\) sur un pas de temps \(\delta t_n\) et “pénaliser” la vitesse comme dans Varra et al. (2024).
Soit \(\overrightarrow{u}^*\) la solution du problème Equation 1.7, on peut alors résoudre le problème Equation 1.8 et pénaliser la vitesse par : \[ \overrightarrow{u} = \frac{2\overrightarrow{u}^*}{1+\sqrt{1 + 4\beta}} \]
où
\[\beta = \delta t_n g C_f h^{-q}\]

avec \(q=4/3\) pour le modèle de Manning-Strickler et \(q=1\) pour le modèle de Darcy-Weisbach.

1.2.5 Solveur de Riemann pour le modèle de Saint-Venant

1.2.5.1 Généralités

Les flux numériques dans CERF sont calculés à l’aide d’un solveur de Riemann. Il existe de très nombreuses références concernant ce solveur comme par exemple Toro (1997). Ce chapitre décrit le solveur de Riemann utilisé pour le modèle de Saint-Venant en s’appuyant sur le cours de Gloria Faccanoni.

On cherche donc à résoudre le problème de Saint-Venant en une dimension. On considère un écoulement en eaux peu profondes de hauteur \(h(x,t)\), de vitesse horizontale \(u(x,t)\) sur un fond plat décrit par le problème suivant : \[ \left \{ \begin{array}{ll} \displaystyle\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} = 0 \\ \displaystyle\frac{\partial (hu)}{\partial t} + \frac{\partial (hu^2 + \frac{1}{2}gh^2)}{\partial x} = 0 \end{array} \right. \tag{1.9}\] Avec les conditions initiales telles que : \[ \overrightarrow{w_p}(x,0) =\left\{ \begin{array}{ll} \overrightarrow{w_p}_L &\mbox{ si } x<0 \\ \overrightarrow{w_p}_R &\mbox{ si } x>0 \end{array} \right. \] où \(\overrightarrow{w_p}_L=<h_L,u_L>^T\) et \(\overrightarrow{w_p}_R=<h_R,u_R>^T\) sont deux états constants.

En développant l’Equation 1.9, on peut réécrire le système sous la forme suivante : \[ \frac{\partial}{\partial t} \left\{ \begin{array}{c} h \\ u \end{array} \right\} + \left[ \begin{array}{cc} u & h \\ g & u \end{array} \right] \frac{\partial}{\partial t} \left\{ \begin{array}{c} h \\ u \end{array} \right\} = \left\{ \begin{array}{c} 0 \\ 0 \end{array} \right\} \]

On cherche alors les valeurs propres \(\lambda_1\) et \(\lambda_2\) ainsi que les vecteurs propres \(\overrightarrow{r}_1\) et \(\overrightarrow{r}_2\) de la matrice ainsi mise en évidence. Un calcul rapide nous donne : \[ \lambda_1=u-\sqrt{gh} , \lambda_2=u+\sqrt{gh} \] avec \(\lambda_1 \leq \lambda_2\) et \[\overrightarrow{r}_1=\left\{ \begin{array}{c} -\sqrt{h} \\ \sqrt{g} \end{array} \right\} , \overrightarrow{r}_2=\left\{ \begin{array}{c} \sqrt{h} \\ \sqrt{g} \end{array} \right\}\].

Note 1.1

Les valeurs propres permettent de déterminer la vitesse de propagation des ondes dans la cellule Important 1.1 comme : \(\lambda_k = Max\left( |u-\sqrt{gh}| , |u+\sqrt{gh}| \right)\)

Pour déterminer la nature des deux champs caractéristiques, on calcule alors \(\nabla_{w_p}\lambda_1 \cdot \overrightarrow{r}_1 = 3/2 \neq 0\) et \(\nabla_{w_p}\lambda_2 \cdot \overrightarrow{r}_2 = 3/2 \neq 0\). Les champs caractéristiques sont donc vraiment non linéaires. Il n’y aura donc pas de discontinuité de contact et les deux ondes sont soit des chocs soit des ondes de raréfaction.

Comme il y a \(2\) ondes, on peut définir pour chaque onde \(1\) invariant de Riemann \(I_{1;1}\), \(I_{2;1}\) tels que \(\nabla_{w_p}I_{1;1} \cdot \overrightarrow{r}_1 = 0\) et \(\nabla_{w_p}I_{2;1} \cdot \overrightarrow{r}_2 = 0\). On peut aisément vérifier que \(I_{1;1}=u+2\sqrt{gh}\) et \(I_{2;1}=u-2\sqrt{gh}\) sont des invariants de Riemann.

La solution faible entropique du problème de Riemann est constituée d’au plus trois états constant \(\overrightarrow{w_p}_L\), \(\overrightarrow{w_p}_*\), \(\overrightarrow{w_p}_R\) séparés par deux discontinuités. Ces discontinuités sont soit des chocs soit des ondes de raréfaction. Par la suite, nous allons déterminer Les états \(\overrightarrow{w_p}_*\).

1.2.5.2 Étude des détentes

Considérons un état \(\overrightarrow{w_p}_g\) et un état \(\overrightarrow{w_p}_d\), reliés par une onde de raréfaction k. On doit alors satisfaire les conditions suivantes :

  • \(\lambda_k(\overrightarrow{w_p}_g)<\lambda_k(\overrightarrow{w_p}_d)\)
  • L’invariant de Riemann est conservé le long de l’onde de raréfaction \(I_{k;1}(\overrightarrow{w_p}_g)=I_{k;1}(\overrightarrow{w_p}_d)\)

Appliquons ces conditions au problème de Saint-Venant

Considérons en premier lieu une 1-détente reliant un état \(\overrightarrow{w_p}_L\) à un état \(\overrightarrow{w_p}_*\), séparés par une onde de raréfaction. Les invariants de Riemann donnent \[ I_{k;1}(\overrightarrow{w_p}_L)=I_{k;1}(\overrightarrow{w_p}_*) \rightarrow u_L+2\sqrt{g\ h_L}=u_*+2\sqrt{g\ h_*} \] soit \[ u_* = u_L + 2\sqrt{g} \left( \sqrt{h_L} - \sqrt{h_*}\right) \tag{1.10}\] De plus, \[\lambda_k(\overrightarrow{w_p}_L)<\lambda_k(\overrightarrow{w_p}_*) \rightarrow u_L-\sqrt{g\ h_L} < u_*-\sqrt{g\ h_*} \rightarrow h_* < h_L\] et \[ u_* > u_L\]

Considérons enfin une 2-détente reliant un état \(\overrightarrow{w_p}_*\) à un état \(\overrightarrow{w_p}_R\), séparés par une onde de raréfaction. Les invariants de Riemann donnent \[ I_{k;1}(\overrightarrow{w_p}_*)=I_{k;1}(\overrightarrow{w_p}_R) \rightarrow u_*-2\sqrt{g\ h_*}=u_R-2\sqrt{g\ h_R} \] soit \[ u_* = u_R + 2\sqrt{g} \left(\sqrt{h_*} - \sqrt{h_R}\right) \tag{1.11}\] De plus, \[\lambda_k(\overrightarrow{w_p}_*)<\lambda_k(\overrightarrow{w_p}_R) \rightarrow u_*+\sqrt{g\ h_*} < u_R+\sqrt{g\ h_R} \rightarrow h_* < h_R\] et \[ u_* < u_R\]

1.2.5.3 Étude des chocs

Considérons un état \(\overrightarrow{w_p}_g\) et un état \(\overrightarrow{w_p}_d\), reliés par une onde de choc \(k\). On doit alors satisfaire les conditions suivantes:

  • Les conditions de Rankine-Hugoniot où l’on note \(\dot\sigma_k\) la vitesse du choc : \(\displaystyle \dot\sigma_k=\frac{\left[\overrightarrow{F}(\overrightarrow{w})\right]}{\left[\overrightarrow{w}\right]}\)
  • Les conditions entropiques de Lax pour un champs \(k\) vraiment non linéaire : \[ \left\{ \begin{array}{l} \lambda_{k-1}(\overrightarrow{w_p}_g) < \dot\sigma_k \\ \lambda_{k}(\overrightarrow{w_p}_d) < \dot\sigma_k < \lambda_{k}(\overrightarrow{w_p}_g) \\ \dot\sigma_k < \lambda_{k+1}(\overrightarrow{w_p}_d) \end{array} \right. \]

Appliquées au problème de Saint-Venant, les conditions de Rankine-Hugoniot donnent : \[ \dot\sigma_k = \frac{h_d u_d - h_g u_g}{h_d - h_g} = \frac{h_d u_d^2 + \frac{g}{2}h_d^2 - h_g u_g^2 + \frac{g}{2}h_g^2}{h_d u_d - h_g u_g} \] on peut définir en particulier \(j_k\) tel que

\[ j_k = h_g u_g - \dot\sigma_k h_g = h_d u_d - \dot\sigma_k h_d .\]

Par élimination de la vitesse du choc, on obtient : \[\left( u_g - u_d \right)^2 = \frac{g}{2 h_d h_g} \left( h_g - h_d \right)^2 \left( h_g + h_d \right)\] et donc \[ j_k = \frac{g}{2} \frac{h_g - h_d}{u_d - u_g} \]

Considérons en premier lieu un 1-choc reliant un état \(\overrightarrow{w_p}_L\) à un état \(\overrightarrow{w_p}_*\), séparés par une onde de choc. Les conditions de Lax donnent : \[ \displaystyle\left\{ \begin{array}{l} \lambda_1(\overrightarrow{w_p}_*) < \dot\sigma_1 < \lambda_1(\overrightarrow{w_p}_L) \\ \dot\sigma_1 < \lambda_2(\overrightarrow{w_p}_*) \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_1 < u_L-\sqrt{g h_L} \\ \dot\sigma_1 < u_*+\sqrt{g h_*} \end{array} \right.\]

Les conditions de Rankine-Hugoniot donnent : \[ \displaystyle\left\{ \begin{array}{l} j_1 = h_*\left( u_*-\dot\sigma_1\right) = h_L\left( u_L-\dot\sigma_1\right) \\ \displaystyle\dot\sigma_1 = \frac{h_* u_*^2 + \frac{g}{2}h_*^2 - h_L u_L^2 - \frac{g}{2}h_L^2}{h_* u_* - h_L u_L} \end{array} \right.\]

Donc à partir des conditions de Lax et en utilisant la première condition de Rankine Hugonio, on obtient : \[ \displaystyle\left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_1 \\ u_L-\sqrt{g h_L} > \dot\sigma_1 \\ u_*+\sqrt{g h_*} > \dot\sigma_1 \end{array} \right. \rightarrow \left\{ \begin{array}{l} h_*\sqrt{g h_*} > j_1 \\ h_L\sqrt{g h_L} < j_1 \\ - h_*\sqrt{g h_*} < j_1 \end{array} \right.\] Soit \[ h_* > h_L \] et comme \(j_1=\frac{g}{2} \frac{h_L - h_*}{u_* - u_L} > 0\), on a \[ u_* < u_L \] et en utilisant l’expression de \(\dot\sigma_1\), on trouve la vitesse \[ u_*=u_L+\left(h_L-h_*\right)\sqrt{\frac{g}{2}\frac{h_L + h_*}{h_L h_*}} \tag{1.12}\]

Considérons enfin un 2-choc reliant un état \(\overrightarrow{w_p}_*\) à un état \(\overrightarrow{w_p}_R\), séparés par une onde de choc. Les conditions de Lax donnent : \[ \displaystyle\left\{ \begin{array}{l} \lambda_1(\overrightarrow{w_p}_*) < \dot\sigma_2 \\ \lambda_2(\overrightarrow{w_p}_R) < \dot\sigma_2 < \lambda_2(\overrightarrow{w_p}_*) \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_2 \\ < u_R+\sqrt{g h_R} < \dot\sigma_2 < u_*+\sqrt{g h_*} \end{array} \right.\]

Les conditions de Rankine-Hugoniot donnent : \[ j_2 = h_*\left( u_*-\dot\sigma_2 \right) = h_R \left( u_R-\dot\sigma_2 \right)= \frac{g}{2} \frac{h_* - h_R}{u_R - u_*}\]

Donc à partir des conditions de Lax et en utilisant la première condition de Rankine Hugonio, on obtient : \[ \displaystyle\left\{ \begin{array}{l} u_*-\sqrt{g h_*} < \dot\sigma_2 \\ u_R+\sqrt{g h_R} < \dot\sigma_2 \\ u_*+\sqrt{g h_*} > \dot\sigma_2 \end{array} \right. \rightarrow \left\{ \begin{array}{l} h_*\sqrt{g h_*} > j_2 \\ -h_R\sqrt{g h_R} > j_2 \\ - h_*\sqrt{g h_*} > j_2 \end{array} \right.\]

Soit \[ h_* > h_R \] et comme \(j_2 < 0\), on a \[ u_* > u_L \] et en utilisant l’expression de \(\dot\sigma_2\), on trouve la vitesse \[ u_* = u_R + \left( h_* - h_R \right) \sqrt{\frac{g}{2} \frac{h_R + h_*}{h_R h_*}} \tag{1.13}\]

1.2.5.4 Construction de la solution

Pour déterminer la solution du problème de Riemann, on doit alors déterminer l’état \(\overrightarrow{w_p}_*\) en fonction des états \(\overrightarrow{w_p}_L\) et \(\overrightarrow{w_p}_R\). On relie l’état intermédiaire à l’état gauche par une 1-onde avec Equation 1.12 et Equation 1.10 : \[ u_* = v^L(h_*) = \left\{ \begin{array}{l} u_L - \displaystyle\left( \sqrt{h_*} - \sqrt{h_L}\right) \frac{2\sqrt{g}}{\sqrt{h_L} + \sqrt{h_*}} \mbox{ si } h_* \leq h_L \mbox{ (1-détente)}\\ u_L - \displaystyle\left( \sqrt{h_*} - \sqrt{h_L}\right) \sqrt{\frac{g}{2} \frac{h_* + h_L}{h_* h_L}} \mbox{ si } h_* > h_L \mbox{ (1-choc)}\end{array} \right. \] et on relie l’état intermédiaire à l’état droit par une 2-onde avec Equation 1.11 et Equation 1.13 : \[ u_* = v^R(h_*) = \left\{ \begin{array}{l} u_R + \displaystyle\left( \sqrt{h_*} - \sqrt{h_R}\right) \frac{2\sqrt{g}}{\sqrt{h_R} + \sqrt{h_*}} \mbox{ si } h_* \leq h_R \mbox{ (2-détente)}\\ u_R + \left( \sqrt{h_*} - \sqrt{h_R}\right) \displaystyle\sqrt{\frac{g}{2} \frac{h_* + h_R}{h_* h_R}} \mbox{ si } h_* > h_R \mbox{ (2-choc)} \end{array} \right. \]

Donc pour trouver \(h_*\), on doit résoudre l’équation \(v^L(h_*)-v^R(h_*)=0\) par un schéma de Newton-Raphson.

En introduisant la fonction \(\zeta(h,\psi)=\left\{ \begin{array}{l} \displaystyle\frac{2\sqrt{g}}{\sqrt{h} + \sqrt{\psi}} \mbox{ si } h \leq \psi \\ \displaystyle\sqrt{\frac{g}{2} \frac{h + \psi}{h \psi}} \mbox{ si } h > \psi \end{array} \right.\), on notera simplement \[u_* = u_R + \left( \sqrt{h_*} - \sqrt{h_R}\right) \zeta(h_*,h_R)\]

Les calculs précédents nous permettent de définir les vitesses des ondes de choc \(\dot\sigma_1 = u_L - h_* \displaystyle\sqrt{\frac{g}{2} \frac{h_L + h_*}{h_L h_*}}\) et \(\dot\sigma_2 = u_R + h_* \displaystyle\sqrt{\frac{g}{2} \frac{h_R + h_*}{h_R h_*}}\) .

et les vitesses des ondes de raréfaction \(\lambda_1 = u_L - \sqrt{g h_L}\) , \(\lambda_1^* = u_* - \sqrt{g h_*}\), \(\lambda_2^* = u_* + \sqrt{g h_*}\) et \(\lambda_2 = u_R + \sqrt{g h_R}\).

La solution dans une zone de détente est donnée par: \[\overrightarrow{w_p}_{1det}(x,t) = \left\{ \begin{array}{l} \displaystyle \frac{1}{9g} \left( u_L+2 \sqrt{g h_L} - \frac{x}{t} \right)^2 \\ \displaystyle \frac{1}{3} \left( u_L + 2 \sqrt{g h_L} + 2 \frac{x}{t} \right) \end{array} \right\}\] et \[\overrightarrow{w_p}_{2det}(x,t) = \left\{ \begin{array}{l} \displaystyle \frac{1}{9g} \left( -u_R + 2 \sqrt{g h_R} + \frac{x}{t} \right)^2 \\ \displaystyle \frac{1}{3} \left( u_R - 2 \sqrt{g h_R} + 2 \frac{x}{t} \right) \end{array} \right\} .\]

Enfin nous pouvons écrire la solution complète du problème de Riemann.

Cas 1: \(h_* > h_L\) et \(h_* > h_R\) la solution est composée d’un 1-choc et d’un 2-choc :

\[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ si } x < \dot\sigma_1 t \\ \overrightarrow{w_p}_* & \mbox{ si } \dot\sigma_1 t < x < \dot\sigma_2 t \\ \overrightarrow{w_p}_R & \mbox{ si } x > \dot\sigma_2 t \end{array} \right.\]

Cas 2: \(h_L < h_* < h_R\) la solution est composée d’un 1-choc et d’une 2-détente : \[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ si } x < \dot\sigma_1 t \\ \overrightarrow{w_p}_* & \mbox{ si } \dot\sigma_1 t < x < \lambda_2^* t \\ \overrightarrow{w_p}_{2det}(x,t) & \mbox{ si } \lambda_2^* t < x < \lambda_2^R t \\ \overrightarrow{w_p}_R & \mbox{ si } x > \lambda_2^R t \end{array} \right.\]

Cas 3: \(h_R < h_* < h_L\) la solution est composée d’une 1-détente et d’un 2-choc : \[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ si } x < \lambda_1^L t \\ \overrightarrow{w_p}_{1det}(x,t) & \mbox{ si } \lambda_1^L t < x < \lambda_1^* t \\ \overrightarrow{w_p}_* & \mbox{ si } \lambda_1^* t < x < \dot\sigma_2 t \\ \overrightarrow{w_p}_R & \mbox{ si } x > \dot\sigma_2 t \end{array} \right.\]

Cas 4: \(h_* < h_L\) et \(h_* < h_R\) la solution est composée d’une 1-détente et d’une 1-détente : \[\overrightarrow{w_p}(x,t) = \left\{ \begin{array}{ll} \overrightarrow{w_p}_L & \mbox{ si } x < \lambda_1^L t \\ \overrightarrow{w_p}_{1det}(x,t) & \mbox{ si } \lambda_1^L t < x < \lambda_1^* t \\ \overrightarrow{w_p}_* & \mbox{ si } \lambda_1^* t < x < \lambda_2^* t \\ \overrightarrow{w_p}_{2det}(x,t) & \mbox{ si } \lambda_2^* t < x < \lambda_2^R t\\ \overrightarrow{w_p}_R & \mbox{ si } x > \lambda_2^R t \end{array} \right.\]

Note

Dans le cas où il apparaît des zones sèches l’algorithme est adapté comme explicité par exemple dans ce cours de Gloria Faccanoni ou la thèse de Pons (2018).

Note

Ce solveur de Riemann est utilisé dans CERF. Il est programmé en Fortran90 dans ~/sources/PHY/svriemann.F90.

Important 1.2

Dans la pratique, il est nécessaire de définir une valeur seuil définissant la zone “sèche”. Dans CERF, on définit donc \(h_{crit}\) tel que si \(h<h_{crit}\), alors on considère que la zone est sèche et on impose \(h =0\) et \(\overrightarrow{u} = \overrightarrow{0}\). Cette valeur est modifiable par l’utilisateur.

1.3 Modèle d’écoulement bifluide

1.3.1 Équations constitutives

Ce modèle a été développé à l’origine pour simuler des écoulements en hydrodynamique maritime, propagation et déferlement de vagues. Pour le déferlement l’approche physiquement la plus pertinente est de considérer les deux fluides air et eau avec un comportement régi par un modèle Navier-Stokes incompressible avec ou sans turbulence. La prise en compte de l’interface est par ailleurs un sujet en soi ! Quelle que soit l’approche de discrétisation utilisée, il s’agit de simulations très coûteuses en temps de calcul. De plus, du fait de dissipation numérique, il est difficile de simuler des écoulements à très grandes échelles. C’est pourquoi, pour simuler des phénomènes de propagation de vagues, on préférera souvent des modèles de type Saint-Venant, Serre-Green-Nagdhi ou écoulement irrotationnel non visqueux. Ces modèles, par contre, modélisent mal, ou pas, les phénomènes de déferlement. La littérature sur ces sujets est très vaste et nous ne prétendons pas faire un état de l’art ici. Nous nous contenterons de proposer une étude comparative partielle dans Helluy, Philippe et al. (2005).
Le modèle bifluide présenté ici est issu de la volonté de trouver un compromis entre la pertinence des modèles de type Navier-Stokes et la simplicité des modèles de type Saint-Venant. On considère donc un mélange de deux fluides, si bien qu’aucun “suivi” d’interface n’est nécessaire et aucune tension de surface ne peut être prise en compte. On fait de plus l’hypothèse, qui se vérifie dans la plupart des cas étudiés en hydrodynamique, que la viscosité peut être négligée, ce qui d’un point de vue numérique nous affranchit avantageusement de la modélisation de l’opérateur Laplacien. Enfin, comme dans Chorin (1967), on considère que les deux fluides sont très faiblement compressibles. Ainsi, on s’affranchit de la nécessité d’imposer l’incompressibilité, ce qui est coûteux en temps de calcul, par contre il devient nécessaire d’adjoindre une équation d’état pour le mélange. Après étude, d’après Golay and Helluy (2007), il semble qu’une équation d’état isotherme soit suffisante pour la plupart des cas étudiés, si bien que l’on peut également s’affranchir de la résolution de l’équation sur l’énergie. On obtient ainsi un modèle Eulérien de la forme :

\[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \overrightarrow{u}) = 0 \tag{1.14}\] \[ \frac{\partial \rho\overrightarrow{u}}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{u} \otimes \overrightarrow{u} + p \overline{\overline{I}} \right) = \rho\overrightarrow{g} \tag{1.15}\]

où \(\rho\) est la masse volumique du mélange, \(\overrightarrow{u}=<u,v,w>^T\) est la vitesse du mélange, \(\overrightarrow{g}\) est le vecteur de gravité et où l’on note \(\overline{\overline{I}}\) la matrice identité, \(\nabla \cdot\) l’opérateur divergence et \(\otimes\) le produit tensoriel. Equation 1.14 représente la conservation de la masse et Equation 1.15 la conservation de la quantité de mouvement. \(p\) est la pression du mélange qui est déterminée par l’équation d’état :
\[ p = p_0 + c_0^2 \left(\rho - \left( \phi \rho_a + (1 - \phi) \rho_w \right) \right) \tag{1.16}\]
où \(\rho_a\) représente la densité de l’air, \(\rho_w\) la densité de l’eau, \(p_0\) la pression de référence et \(c_0\) la vitesse du son dans le mélange. \(\phi\) représente la fraction volumique de l’air dans le mélange. \(\phi\) satisfait une équation de transport, que l’on prendra sous la forme non conservative (Golay and Helluy (2007)) de la forme :
\[ \frac{\partial \phi}{\partial t} + \overrightarrow{u} \cdot \nabla \phi = 0 \tag{1.17}\]

Usuellement, on prendra \(p_0 = 10^5 Pa\), \(\rho_a = 1. kg/m^3\), \(\rho_w = 1000. kg/m^3\) et \(c_0 = 20. m/s\).
En hydrodynamique, les écoulements sont de l’ordre de quelques \(m/s\), si bien qu’une vitesse du son de l’ordre de \(20 m/s\) conduit à un nombre de Mach \(M = \displaystyle \frac{||\overrightarrow{u}||}{c_0} < 0.3\). On peut alors considérer que le mélange est quasiment incompressible. Cette vitesse du son est complètement artificielle et ne correspond à aucune réalité physique. Elle est choisie pour des raisons numériques. En effet, le choix de \(c_0\) va influencer la stabilité du schéma numérique explicite utilisé pour résoudre les équations par la condition de CFL (Important 1.1).

Cependant, de façon fortuite, d’après l’équation de Wood (Figure 1.9), alors que la vitesse du son dans l’eau est de l’ordre de 1481 \(m/s\) et dans l’air de l’ordre de 343 \(m/s\), la vitesse du son dans un mélange air-eau est de l’ordre de 23 \(m/s\). Il ne s’agit nullement d’une justification du modèle utilisé ici, mais il est intéressant de noter que la vitesse du son dans un mélange air-eau est de l’ordre de la vitesse du son artificielle utilisée dans le modèle bifluide.

Figure 1.9: Vitesse du son dans un mélange air-eau d’après l’équation de Wood

Les équations Equation 1.14, Equation 1.15 et Equation 1.17, complétées par Equation 1.16 constituent le problème hyperbolique résolu par CERF pour modéliser les écoulements bifluides. On définit les variables conservatives \(\overrightarrow{w}=<\rho,\rho u,\rho v,\rho w,\phi>^T\) et les variables primitives \(\overrightarrow{w_p}=<\rho, u, v, w, p>^T\). On retrouve la forme générale Equation 1.1.

1.3.2 Critère de raffinement

Dans le cas du modèle bifluide, pour définir la production numérique d’entropie Equation 1.4, on utilise l’entropie définie par: \[ s(\overrightarrow{w}) = \frac{1}{2}\rho ||\overrightarrow{u}||^2 + c_0^2 \rho ln(\rho) - c_0^2 \left( \rho_w - \rho_a \right) \phi\] et le flux d’entropie \[ \overrightarrow{\psi}_s(\overrightarrow{w}) = \left( s + c_0^2 \rho + c_0^2 \left( \rho_w - \rho_a \right) \right) \overrightarrow{u} = \left( \frac{1}{2}\rho ||\overrightarrow{u}||^2 + c_0^2 \rho \left( ln(\rho) + 1 \right) \right) \overrightarrow{u} .\]

1.3.3 Raidissement de la surface libre

Le schéma de résolution par volumes finis du système hyperbolique Equation 1.14, Equation 1.15 et Equation 1.17, même à l’ordre 2, est connu pour être dissipatif. En conséquence l’interface air-eau diffuse et s’élargit. Pour contrer ce phénomène, on peut utiliser un modèle de raidissement de l’interface. L’idée est de rajouter un terme de source à l’équation de transport de \(\phi\) Equation 1.17 qui va pénaliser les zones de mélange avec une méthode de splitting. Sur un pas de temps fictif on calcule donc : \[ \frac{\partial \phi}{\partial \tau} = \phi^2 (1 - \phi)^2(\phi - c)\] où \(c\) est déterminé de façon à conserver la fraction volumique, c’est à dire:
\[ c = \frac{\int _{\Omega} \phi^3 (1 - \phi)^2 d\Omega}{\int _{\Omega} \phi^2 (1 - \phi)^2 d\Omega}\]

On peut remarquer que pour \(\phi=0\) ou \(\phi=1\), ce terme source est sans effet. Enfin, comme \(\phi\) est modifié, par soucis de conservation on modifie également la masse volumique et la quantité de mouvement du mélange :
\[ \frac{\partial \rho}{\partial \tau} = \phi^2 (1 - \phi)^2(\phi - c)(\rho_a - \rho_w)\] \[ \frac{\partial \left(\rho \overrightarrow{u} \right)}{\partial \tau} = \phi^2 (1 - \phi)^2(\phi - c) \overrightarrow{u}(\rho_a - \rho_w)\]

Le pas de temps fictif est calibré pour représenter une part du flux. Cette part est déterminée par un coefficient de raidissement \(\alpha_{sharp}\). Dans la pratique, ce coefficient est de l’ordre de \(10^{-3} - 10^{-2}\).

1.3.4 Notion de “fluide rigide”

En suivant le principe des domaines fictifs, on peut pénaliser une partie des cellules pour simuler un obstacle en imposant une vitesse nulle. En suivant les travaux de Coquerelle and Cottet (2008) on peut également pénaliser une partie des cellules pour obtenir un comportement de type “solide rigide”. L’extension des travaux de Coquerelle and Cottet (2008) à un modèle bifluide a été réalisé par Altazin (2017). Considérons un “solide” occupant le domaine \(\Omega_s\), de densité uniforme \(\rho_s\), de centre de gravité \(G\), de masse \(M_s\), de matrice d’inertie \(\overline{\overline{J_s}}\). La vitesse de tout point \(P\) du solide est donnée par :
\[ \overrightarrow{v}_s = \overrightarrow{v}_G + \overrightarrow{\omega}_G \wedge \overrightarrow{GP}\]
où \(\overrightarrow{v}_G\) est la vitesse du centre de gravité, \(\overrightarrow{\omega}_G\) est la vitesse angulaire du solide et \(\overrightarrow{GP}\) est le vecteur reliant le centre de gravité au point \(P\). En minimisant \(\int _{\Omega_s} || \rho\overrightarrow{u} - \rho_s\overrightarrow{v}_s||^2 d\Omega_s\), on obtient pour une densité uniforme du solide : \[\overrightarrow{v}_G=\frac{\int _{\Omega_s} \rho\overrightarrow{u} d\Omega_s}{M_s}\]
et
\[\overrightarrow{\omega}_G = \overline{\overline{J_s}}^{-1} \int _{\Omega_s} \overrightarrow{GP} \wedge \overrightarrow{u} .\]
À l’issue de chaque pas de temps, on calcule la position et la vitesse du centre de gravité du solide et la vitesse de rotation du solide. On affecte ensuite \(\rho_s\) et \(\overrightarrow{v}_s\) à toutes les cellules du “solide”. Dans CERF le solide est déterminé par son enveloppe. Cette enveloppe est définie par un maillage triangulaire. Les cellules à l’intérieur de l’enveloppe sont déterminées par une technique de “ray casting”. Dans CERF on peut considérer un solide “fixe”, à “vitesse imposée” ou “libre”.

1.3.5 Solveur de Riemann pour le modèle bifluide

1.3.5.1 Généralités

Les flux numériques dans CERF sont calculés à l’aide d’un solveur de Riemann. Il existe de très nombreuses références concernant ce solveur comme par exemple Toro (1997). Ce chapitre décrit le solveur de Riemann utilisé pour le modèle d’écoulement bifluide en s’appuyant sur Golay and Helluy (2007).

On cherche donc à résoudre le problème bifluide isotherme en une dimension. On considère un écoulement de densité \(\rho(x,t)\), de vitesse horizontale \(u(x,t)\), de pression \(p(x,t)\), de fraction volumique \(\phi(x,t)\) décrit par le problème suivant : \[ \left \{ \begin{array}{ll} \displaystyle\frac{\partial \rho}{\partial t} + \frac{\partial (\rho u)}{\partial x} = 0 \\ \displaystyle\frac{\partial (\rho u)}{\partial t} + \frac{\partial (\rho u^2 + p)}{\partial x} = 0 \\ \displaystyle\frac{\partial \rho \phi}{\partial t} + \frac{\partial \rho \phi u}{\partial x} = 0 \\ \end{array} \right. \tag{1.18}\] Avec une équation d’état de la forme : \[ p = p_0 + c_0^2 \left(\rho - \left( \phi \rho_a + (1 - \phi) \rho_w \right) \right) \] et les conditions initiales telles que : \[ \overrightarrow{w_p}(x,0) =\left\{ \begin{array}{ll} \overrightarrow{w_p}_L &\mbox{ si } x<0 \\ \overrightarrow{w_p}_R &\mbox{ si } x>0 \end{array} \right. \] où \(\overrightarrow{w_p}_L=<\rho_L, u_L, \phi_L>^T\) et \(\overrightarrow{w_p}_R=<\rho_R, u_R, \phi_R>^T\) sont deux états constants.

En développant l’Equation 1.18, on peut réécrire le système sous la forme suivante : \[ \frac{\partial}{\partial t} \left\{ \begin{array}{c} \rho \\ u \\ \phi \end{array} \right\} + \left[ \begin{array}{ccc} u & \rho & 0 \\ \displaystyle\frac{c_0^2}{\rho} & u & \displaystyle\frac{c_0^2\left( \rho_w - \rho_a \right)}{\rho} \\ 0 & 0 & u \end{array} \right] \frac{\partial}{\partial x} \left\{ \begin{array}{c} \rho \\ u \\ \phi \end{array} \right\} = \left\{ \begin{array}{c} 0 \\ 0 \\ 0\end{array} \right\} . \]

On cherche alors les valeurs propres \(\lambda_1\), \(\lambda_2\) et \(\lambda_3\) ainsi que les vecteurs propres \(\overrightarrow{r}_1\), \(\overrightarrow{r}_2\) et \(\overrightarrow{r}_3\) de la matrice ainsi mise en évidence. Un calcul rapide nous donne: \[ \lambda_1=u - c_0 , \lambda_2=u, \lambda_3=u + c_0 \] avec \(\lambda_1 \leq \lambda_2 \leq \lambda_3\) et \[\overrightarrow{r}_1=\left\{ \begin{array}{c} -\rho/c_0 \\ 1 \\0 \end{array} \right\} , \overrightarrow{r}_2=\left\{ \begin{array}{c} \left( \rho_w - \rho_a \right) \\ 0 \\ -1 \end{array} \right\}, \overrightarrow{r}_3=\left\{ \begin{array}{c} \rho/c_0 \\ 1 \\0 \end{array} \right\} .\]

Note 1.2

Les valeurs propres permettent de déterminer la vitesse de propagation des ondes dans la cellule Important 1.1 comme : \(\lambda_k = Max\left( |u - c_0| , |u|, |u + c_0| \right)\)

Pour déterminer la nature des trois champs caractéristiques, on calcule alors \(\nabla_{w_p}\lambda_1 \cdot \overrightarrow{r}_1 = 1 \neq 0\), \(\nabla_{w_p}\lambda_2 \cdot \overrightarrow{r}_2 = 0\) et \(\nabla_{w_p}\lambda_3 \cdot \overrightarrow{r}_3 = 1 \neq 0\). Les champs caractéristiques 1 et 3 sont donc vraiment non linéaires alors que le champs 2 est linéairement dégénéré. Il y aura donc une de discontinuité de contact et deux ondes qui sont soit des chocs soit des ondes de raréfaction.

Comme il y a \(3\) ondes, on peut définir pour chaque onde \(2\) invariants de Riemann. Les invariants de Riemann de l’onde \(k\) sont tels que \(\nabla_{w_p}I_{k;1} \cdot \overrightarrow{r}_k = 0\) et \(\nabla_{w_p}I_{k;2} \cdot \overrightarrow{r}_k = 0\). On peut aisément vérifier que :
\(I_{1;1} = \phi\) , \(I_{1;2} = u + c_0 ln(\rho)\),
\(I_{2;1} = u\) , \(I_{2;2} = p\),
\(I_{3;1} = \phi\) , \(I_{3;2} = u - c_0 ln(\rho)\).

La solution faible entropique du problème de Riemann est constituée d’au plus quatre états constants \(\overrightarrow{w_p}_L\), \(\overrightarrow{w_p}_I\), \(\overrightarrow{w_p}_{II}\), \(\overrightarrow{w_p}_R\) séparés par un choc ou une détente, une discontinuité de contact, puis un choc ou une détente (Figure 1.10). Par la suite, nous allons déterminer Les états \(\overrightarrow{w_p}_I\) et \(\overrightarrow{w_p}_{II}.\)

Figure 1.10: Schéma des solutions du problème de Riemann avec le modèle bi-fluide isotherme

1.3.5.2 Étude des détentes

Considérons un état \(\overrightarrow{w_p}_g\) et un état \(\overrightarrow{w_p}_d\), reliés par une onde de raréfaction k. On doit alors satisfaire les conditions suivantes:

  • \(\lambda_k(\overrightarrow{w_p}_g)<\lambda_k(\overrightarrow{w_p}_d)\)
  • Les invariants de Riemann sont conservés le long de l’onde de raréfaction \(I_{k;i}(\overrightarrow{w_p}_g)=I_{k;i}(\overrightarrow{w_p}_d)\) pour \(i=1,2\)

Appliquons ces conditions au problème bifluide isotherme.

Considérons en premier lieu une 1-détente reliant un état \(\overrightarrow{w_p}_L\) à un état \(\overrightarrow{w_p}_I\), séparés par une onde de raréfaction. Les invariants de Riemann donnent :
\[\left\{ \begin{array}{l}I_{1;1}(\overrightarrow{w_p}_L) = I_{1;1}(\overrightarrow{w_p}_I) \\ I_{1;2}(\overrightarrow{w_p}_L) = I_{1;2}(\overrightarrow{w_p}_I)\end{array} \right. \rightarrow \left\{ \begin{array}{l}\phi_L=\phi_I \\ u_L + c_0 ln(\rho_L) = u_I + c_0 ln(\rho_I)\end{array} \right.\] soit \[ u_I = u_L + c_0 ln\left( \frac{\rho_L}{\rho_I}\right) \tag{1.19}\] De plus, \[\lambda_1(\overrightarrow{w_p}_L)<\lambda_1(\overrightarrow{w_p}_I) \rightarrow u_L < u_I \] et en conséquence d’après Equation 1.19, \(\rho_L > \rho_I\).

Considérons enfin une 3-détente reliant un état \(\overrightarrow{w_p}_{II}\) à un état \(\overrightarrow{w_p}_R\), séparés par une onde de raréfaction. La condition sur les vitesses s’écrit : \[\lambda_3(\overrightarrow{w_p}_{II})<\lambda_3(\overrightarrow{w_p}_R) \rightarrow u_{II} < u_R .\]

Les invariants de Riemann donnent \[\left\{ \begin{array}{l}I_{3;1}(\overrightarrow{w_p}_{II}) = I_{3;1}(\overrightarrow{w_p}_R) \\ I_{3;2}(\overrightarrow{w_p}_{II}) = I_{3;2}(\overrightarrow{w_p}_R)\end{array} \right. \rightarrow \left\{ \begin{array}{l}\phi_{II}=\phi_R \\ u_{II} - c_0 ln(\rho_{II}) = u_R - c_0 ln(\rho_R)\end{array} \right.\] soit \[ u_{II} = u_R - c_0 ln\left( \frac{\rho_R}{\rho_{II}}\right) \tag{1.20}\]
donc comme \(u_{II} < u_R\), on a \(\rho_{II} < \rho_R\).

1.3.5.3 Étude des chocs

Considérons un état \(\overrightarrow{w_p}_g\) et un état \(\overrightarrow{w_p}_d\), reliés par une onde de choc \(k\). On doit alors satisfaire les conditions suivantes :

  • Les conditions de Rankine-Hugoniot où l’on note \(\dot\sigma_k\) la vitesse du choc : \(\displaystyle \dot\sigma_k=\frac{\left[\overrightarrow{F}(\overrightarrow{w})\right]}{\left[\overrightarrow{w}\right]}\)
  • Les condition entropiques de Lax pour un champs \(k\) vraiment non linéaire : \[ \left\{ \begin{array}{l} \lambda_{k-1}(\overrightarrow{w_p}_g) < \dot\sigma_k < \lambda_{k}(\overrightarrow{w_p}_g) \\ \lambda_{k}(\overrightarrow{w_p}_d) < \dot\sigma_k < \lambda_{k+1}(\overrightarrow{w_p}_d) \end{array} \right. \]

Appliquées au problème bifluide, les conditions de Rankine-Hugoniot donnent : \[ \dot\sigma_k = \frac{\rho_d u_d - \rho_g u_g}{\rho_d - \rho_g} = \frac{\rho_d u_d^2 + p_d - \rho_g u_g^2 - p_g}{\rho_d u_d - \rho_g u_g} = \frac{\rho_d u_d \phi_d - \rho_g u_g \phi_g}{\rho_d \phi_d - \rho_g \phi_g} \] \[ \left\{ \begin{array}{l} \rho_d \dot\sigma_k - \rho_g \dot\sigma_k = \rho_d u_d - \rho_g u_g \\ \rho_d u_d \dot\sigma_k - \rho_g u_g \dot\sigma_k = \rho_d u_d^2 + p_d - \rho_g u_g^2 - p_g \\ \rho_d \phi_d \dot\sigma_k - \rho_g \phi_g \dot\sigma_k = \rho_d u_d \phi_d - \rho_g u_g \phi_g\end{array} \right.\] \[ \left\{ \begin{array}{l} \rho_g (u_g - \dot\sigma_k) = \rho_d (u_d - \dot\sigma_k) \\ u_d \rho_d (u_d - \dot\sigma_k) -u_g \rho_g (u_g - \dot\sigma_k) + p_d - p_g = 0 \\ \phi_g \rho_g (u_g - \dot\sigma_k) = \phi_d \rho_d (u_d - \dot\sigma_k) \end{array} \right.\] Soit, si on note \(j_k = \rho_g (u_g - \dot\sigma_k) = \rho_d (u_d - \dot\sigma_k)\), on obtient : \[\phi_g = \phi_d\]
et
\[j_k = \frac{p_d - p_g}{u_g - u_d}\]

Par élimination de la vitesse du choc, on a également : \[\left( \rho_d u_d - \rho_g u_g \right)^2 = \left( \rho_d - \rho_g \right) \left( \rho_d u_d^2 - \rho_g u_g^2 + p_d - p_g \right)\] ou encore
\[\rho_d \rho_g \left( u_d - u_g \right)^2 = \left( \rho_d - \rho_g \right) \left( p_d - p_g \right)\]

Considérons en premier lieu un 1-choc reliant un état \(\overrightarrow{w_p}_L\) à un état \(\overrightarrow{w_p}_I\), séparés par une onde de choc. Les conditions d’admissibilité de Lax s’écrivent : \[ \displaystyle\left\{ \begin{array}{l} \dot\sigma_1 < \lambda_1(\overrightarrow{w_p}_L) \\ \lambda_1(\overrightarrow{w_p}_I) <\dot\sigma_1 < \lambda_2(\overrightarrow{w_p}_I) \end{array} \right. \rightarrow \left\{ \begin{array}{l} \dot\sigma_1 < u_L - c_0 \\ u_I - c_0 < \dot\sigma_1 \\ \dot\sigma_1 < u_I \end{array} \right. \rightarrow \left\{ \begin{array}{l} \rho_L c_0 < j_1 \\ j_1 < \rho_I c_0 \\ \dot\sigma_1 < u_I \end{array} \right.\]

Donc en particulier \(u_I < u_L\) et \(\rho_L < \rho_I\). Les conditions de Rankine-Hugoniot donnent : \[ \displaystyle\left\{ \begin{array}{l} \phi_L = \phi_I \\ \left( u_I - u_L \right)^2 = \displaystyle\frac{\rho_I - \rho_L}{\rho_I \rho_L}\left( p_I - p_L \right) \end{array} \right.\] d’où, comme \(u_I < u_L\) et \(\rho_L < \rho_I\) : \[ u_I=u_L -\sqrt{\frac{\rho_I - \rho_L}{\rho_I \rho_L}\left( p_I - p_L \right)} = u_L + c_0 \frac{\rho_L - \rho_I}{\sqrt{\rho_I \rho_L}} \tag{1.21}\]

Considérons enfin un 3-choc reliant un état \(\overrightarrow{w_p}_{II}\) à un état \(\overrightarrow{w_p}_R\), séparés par une onde de choc. Les conditions de Lax donnent : \[ \displaystyle\left\{ \begin{array}{l} \lambda_2(\overrightarrow{w_p}_{II}) < \dot\sigma_3 < \lambda_3(\overrightarrow{w_p}_{II}) \\ \lambda_3(\overrightarrow{w_p}_R) < \dot\sigma_3 \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_{II} < \dot\sigma_3 \\ \dot\sigma_3 < u_{II} + c_0 \\ u_R + c_0 < \dot\sigma_3 \end{array} \right. \rightarrow \left\{ \begin{array}{l} u_{II} < \dot\sigma_3 \\ - c_0 \rho_{II} < j_3 \\ j_3 < - c_0 \rho_R \end{array} \right.\]

Donc en particulier \(u_R < u_{II}\) et \(\rho_R < \rho_{II}\).Les conditions de Rankine-Hugoniot s’écrivent : \[ \displaystyle\left\{ \begin{array}{l} \phi_{II} = \phi_R \\ \left( u_R - u_{II} \right)^2 = \displaystyle\frac{\rho_R - \rho_{II}}{\rho_R \rho_{II}}\left( p_R - p_{II} \right) \end{array} \right.\] d’où , comme \(u_R < u_{II}\) et \(\rho_R < \rho_{II}\) : \[ u_{II}=u_R + \sqrt{\frac{\rho_R - \rho_{II}}{\rho_R \rho_{II}}\left( p_R - p_{II} \right)} = u_R - c_0 \frac{\rho_R - \rho_{II}}{\sqrt{\rho_R \rho_{II}}} \tag{1.22}\]

1.3.5.4 Étude de la discontinuité de contact

Considérons un état \(\overrightarrow{w_p}_g\) et un état \(\overrightarrow{w_p}_d\), reliés par une discontinuité de contact. Au travers d’une k-discontinuité de contact l’entropie et les invariants de Riemann sont conservés. La discontinuité se déplaçant à la vitesse \(\dot\sigma_k\) on a : \[\displaystyle \dot\sigma_k = \lambda_{k}(\overrightarrow{w_p}_g) = \lambda_{k}(\overrightarrow{w_p}_d) = \frac{\left[\overrightarrow{F}(\overrightarrow{w})\right]}{\left[\overrightarrow{w}\right]}\]

Considérons donc une 2-discontinuité de contact reliant l’état \(\overrightarrow{w_p}_I\) à l’état \(\overrightarrow{w_p}_{II}\). On a alors \[ \dot\sigma_2 = u_I = u_{II} = \frac{\rho_{II} u_{II} - \rho_I u_I}{\rho_{II} - \rho_I} = \frac{\rho_{II} u_{II}^2 + p_{II} - \rho_I u_I^2 - p_I}{\rho_{II} u_{II} - \rho_I u_I} = \frac{\rho_{II} u_{II} \phi_{II} - \rho_I u_I \phi_I}{\rho_{II} \phi_{II} - \rho_I \phi_I} \]

D’où l’on déduit : \[ u_I = u_{II} \mbox{ et } p_I = p_{II} \tag{1.23}\]

1.3.5.5 Construction de la solution

Pour simplifier les notations, définissons une fonction \(H(a,b)\) telle que \(H(a,b) = \displaystyle\left\{ \begin{array}{ll} c_0 \displaystyle\frac{\rho_a - \rho_b}{\sqrt{\rho_a \rho_b}} & \mbox{ si } \rho_a < \rho_b \\ c_0 ln \left(\displaystyle\frac{\rho_a}{\rho_b} \right) & \mbox{ si } \rho_a > \rho_b \end{array} \right.\).
On peut alors écrire pour la première onde, en utilisant Equation 1.19, Equation 1.21 \(u_I = u_L + H(\rho_L, \rho_I)\) et pour la troisième onde, en utilisant Equation 1.20, Equation 1.22 \(u_{II} = u_R - H(\rho_R, \rho_{II})\).

On sait par l’étude de la discontinuité de contact que : \[u_I - u_{II} = 0\] donc
\[u_L + H(\rho_L, \rho_I) - u_R + H(\rho_R, \rho_{II}) = 0\] Or, \(\rho_I=\rho_I(p_I,\phi_I)\) et \(\rho_{II}=\rho_{II}(p_{II},\phi_{II})\). On sait par ailleurs que \(p_I = p_{II}\), nous noterons alors la pression intermédiaire \(p^*\) telle que \(p^* = p_I = p_{II}\). De plus on sait \(\phi_I = \phi_L\) et \(\phi_{II} = \phi_R\). On peut alors écrire : \[u_L - u_R + H(\rho_L, \rho(p^*,\phi_L)) + H(\rho_R, \rho(p^*,\phi_R)) = 0\] C’est donc une équation non-linéaire en \(p^*\) que l’on peut résoudre par une méthode de Newton-Raphson. Connaissant \(p^*\), on peut alors déterminer les états \(\overrightarrow{w_p}_I\) et \(\overrightarrow{w_p}_{II}\).

Note

Ce solveur de Riemann est utilisé dans CERF. Il est programmé en Fortran90 dans ~/sources/PHY/riemisot.F90.