4. Modélisation numérique

4.1. Propriétés mathématiques des équations

Nous avons vu que les équations d’Euler régissant l’écoulement adiabatique non visqueux d’un fluide parfait s’écrivent sous la forme vectorielle pour un vecteur d’état \(\overrightarrow{W}=[W_i]\)

\[\frac{\partial W_{i}}{\partial t}+div(\overrightarrow{F_{i}}(\overrightarrow{W}))=0\]

\(\overrightarrow{W}=\left[\begin{array}{ccccc} \rho & \rho u & \rho v & \rho w & \rho(e+\frac{1}{2}U^{2})\end{array}\right]\) est le vecteur d’état, et \(F_{i}(\overrightarrow{W})\) le vecteur flux associé:

\[\begin{split}\overrightarrow{F}_{\rho}=\left[\begin{array}{c} \rho u\\ \rho v\\ \rho w \end{array}\right],\,\overrightarrow{F}_{u}=\left[\begin{array}{c} \rho u^{2}+p\\ \rho uv\\ \rho uw \end{array}\right],\,\overrightarrow{F}_{v}=\left[\begin{array}{c} \rho uv\\ \rho v^{2}+p\\ \rho vw \end{array}\right]\,\end{split}\]

\[\begin{split}\overrightarrow{F}_{w}=\left[\begin{array}{c} \rho uw\\ \rho vw\\ \rho w^{2}+p \end{array}\right],\,\overrightarrow{F}_{e_{t}}=\left[\begin{array}{c} \rho uh_{t}\\ \rho vh_{t}\\ \rho wh_{t} \end{array}\right]\end{split}\]

où on a noté \(e_{t}=e+\frac{1}{2}U^{2}\) l’énergie totale par unité de masse, et \(h_{t}=e+\frac{p}{\rho}+\frac{1}{2}U^{2}\) l’enthalpie totale par unité de masse.

A ces équations on ajoute l’équation d’état des gaz parfaits:

\[e=\frac{1}{\gamma-1}\frac{p}{\rho}\]

4.1.1. Analyse dans le cas 1D

Dans le cas d’un écoulement uni-dimensionnel, elles s’écrivent:

\[\begin{split}\begin{aligned} \frac{\partial\rho}{\partial t}+\frac{\partial\rho u}{\partial x} & = & 0\\ \frac{\partial\rho u}{\partial t}+\frac{\partial(\rho u^{2}+p)}{\partial x} & = & 0\\ \frac{\partial}{\partial t}(\rho e_{t})+\frac{\partial(\rho e_{t}+p)u}{\partial x} & = & 0\end{aligned}\end{split}\]

soit sous forme matricielle:

\[\frac{\partial\overrightarrow{W}}{\partial t}+\frac{\partial\overrightarrow{F}(\overrightarrow{W})}{\partial x}=0\]

En introduisant la matrice jacobienne des flux: i.e. la dérivée du vecteur flux \(\overrightarrow{F}(\overrightarrow{W})\) par rapport au vecteur d’état \(\overrightarrow{W}=[\rho,\rho u,\rho e_{t}]\)

\[\mathcal{A}=\frac{\partial F}{\partial W}\]

cette équation s’écrit sous la forme quasi-linéaire suivante:

(4.2)\[\frac{\partial\overrightarrow{W}}{\partial t}+\mathcal{A}\frac{\partial\overrightarrow{W}}{\partial x}=0\]

Un calcul formel fournit l’expression de cette Jacobienne

\[\begin{split}\mathcal{A}=\left[\begin{array}{ccc} 0 & 1 & 0\\ \frac{\gamma-3}{2}u^{2} & (3-\gamma)u & \gamma-1\\ u\,(\frac{\gamma-2}{2}u^{2}-\frac{1}{\gamma-1}c^{2}) & \frac{3-2\gamma}{2}u^{2}+\frac{1}{\gamma-1}c^{2} & u\gamma \end{array}\right]\end{split}\]

où \(u\) est la vitesse du fluide (supposée \(>0\)), et \(c=\sqrt{\gamma\frac{p}{\rho}}\) la célérité du son.

Cette matrice possède 3 valeurs propres réelles

\[\lambda=[u-c,\,u,\,u+c]\]

ce qui traduit la présence de 3 courbes caractéristiques le long desquelles se propagent les perturbations.

Ces équations forment donc un système hyperbolique non linéaire. En notant \(\mathcal{R}\) la matrice des vecteurs propres à droite \(\overrightarrow{r_{i}}\):

\[\mathcal{A}\overrightarrow{r_{i}}=\lambda_{i}\overrightarrow{r_{i}}\]

la diagonalisation de \(\mathcal{A}\) s’écrit:

\[\mathcal{A}=\mathcal{R}\Lambda\mathcal{R}^{-1}\]

où \(\Lambda\) est la matrice diagonale des vecteurs propres \(\lambda_{i}\). En multipliant l’équation (4.2) par \(\mathcal{R}^{-1}\), on obtient:

(4.3)\[\mathcal{R}^{-1}\frac{\partial\overrightarrow{W}}{\partial t}+\Lambda\mathcal{R}^{-1}\frac{\partial\overrightarrow{W}}{\partial x}=0\]

On introduit le changement de variable

\[\overrightarrow{V}=\mathcal{R}^{-1}\overrightarrow{W}\]

et le vecteur \(\overrightarrow{V}\) est appelé « variables caractéristiques ». Le système d’équations transformé est appelé « équations caractéristiques »:

\[\mathcal{R}^{-1}\frac{\partial\mathcal{R}\overrightarrow{V}}{\partial t}+\Lambda\mathcal{R}^{-1}\frac{\partial\mathcal{R}\overrightarrow{V}}{\partial x}=0\]

Dans le cas de petites perturbations, on peut simplifier ce système en considérant que la matrice \(\mathcal{R}\) est constante ainsi que les valeurs propres \(\Lambda\):

(4.4)\[\frac{\partial\overrightarrow{V}}{\partial t}+\Lambda\frac{\partial\overrightarrow{V}}{\partial x}=0\]

Ce système caractéristique est un système de 3 équations scalaires indépendantes de type advection:

\[\frac{\partial V_{i}}{\partial t}+\lambda_{i}\frac{\partial V_{i}}{\partial x}=0\]

Cette équation traduit le fait que la composante \(V_{i}\) est constante le long de la droite caractéristique de pente \(\frac{dx}{dt}=\lambda_{i}\), i.e.:

\[V_{i}(x,t)=V_{i}(x-\lambda_{i}t)\]

La solution générale de (4.4) est donc une combinaison linéaire de ces ondes caractéristiques, soit en revenant aux variables d’état \(\overrightarrow{W}\):

\[\overrightarrow{W}=\sum_{i=1}^{3}V_{i}(x-\lambda t)\,\overrightarrow{r_{i}}\]

Les perturbations de la solution se propagent donc le long de ces directions caractéristiques.

Dans le cas non-linéaire, la matrice \(\mathcal{R}\) n’est plus constante, et l’analyse précédente n’est plus exacte. On peut cependant la considérer comme une approximation au premier ordre et en déduire les propriétés de la solution. Le système hyperbolique (4.2) possède 3 valeurs propres \(\lambda_{i}=\lambda_{i}(\overrightarrow{W})\), qui dépendent de la solution \(\overrightarrow{W}\). Associées à ces 3 valeurs propres, on a des courbes caractéristiques \(\mathcal{C}_{i}=\mathcal{C}_{i}(\overrightarrow{W})\) de pente \(\lambda_{i}\) qui sont aussi fonctions de la solution \(\overrightarrow{W}\). Le long de ces caractéristiques se propagent les perturbations de la solution sous forme d’ondes non-linéaires.

Fig. 4.27 propagation d’une perturbation

Considérons une perturbation située en \(x=0\) à l’instant \(t\). Notons \(s\) la vitesse de propagation de la perturbation, \(W_{L}\) l’état du fluide à gauche (left) et \(W_{R}\) l’état à droite (right). La solution \(\overrightarrow{W}\) vérifie les équations d’Euler

\[\frac{\partial\overrightarrow{W}}{\partial t}+\frac{\partial\overrightarrow{F}(\overrightarrow{W})}{\partial x}=0\]

En intégrant cette équation sur un domaine \([-x_{0},+x_{0}]\) autour de la perturbation, on obtient l’équation:

\[\frac{d}{dt}\int_{-x_{0}}^{+x_{0}}\overrightarrow{W}\,dx=\overrightarrow{F}(\overrightarrow{W_{L}})-\overrightarrow{F}(\overrightarrow{W_{R}})\]

D’autre part si \(x_{0}\) est grand devant \(s\,dt\), la variation temporelle de \(\overrightarrow{W}\) dans \([-x_{0},+x_{0}]\) s’écrit:

\[\int_{-x_{0}}^{+x_{0}}\overrightarrow{W}(t+dt)\,dx-\int_{-x_{0}}^{+x_{0}}\overrightarrow{W}(t)\,dx=\left[\overrightarrow{W_{L}}(x_{0}+sdt)+\overrightarrow{W_{R}}(x_{0}-sdt)\right]-\left[\overrightarrow{W_{L}}x_{0}+\overrightarrow{W_{R}}x_{0}\right]\]

d’où:

(4.5)\[\frac{d}{dt}\int_{-x_{0}}^{+x_{0}}\overrightarrow{W}\,dx=s\left[\overrightarrow{W_{L}}-\overrightarrow{W_{R}}\right]\]

On retrouve les relations de saut d’Hugoniot:

(4.6)\[\left[\overrightarrow{F}(\overrightarrow{W_{L}})-\overrightarrow{F}(\overrightarrow{W_{R}})\right]=s\left[\overrightarrow{W_{L}}-\overrightarrow{W_{R}}\right]\]

qui relie la vitesse de propagation de la perturbation \(s\) à l’état de part et d’autre de la perturbation \(\overrightarrow{W_{L}}\) et \(\overrightarrow{W_{R}}\). Cette vitesse de propagation \(s\) n’est donc pas quelconque. En effectuant une analyse de perturbation au premier ordre, i.e. en développant \(\overrightarrow{F}(\overrightarrow{W_{L}})\) au premier ordre au voisinage de \(\overrightarrow{W_{R}}\):

\[\begin{split}\begin{aligned} \overrightarrow{F}(\overrightarrow{W_{L}}) & = & \overrightarrow{F}(\overrightarrow{W_{R}})+\frac{\partial\overrightarrow{F}(\overrightarrow{W})}{\partial\overrightarrow{W}}(\overrightarrow{W_{L}}-\overrightarrow{W_{R}})\\ & = & \overrightarrow{F}(\overrightarrow{W_{R}})+\mathcal{A}(\overrightarrow{W_{L}}-\overrightarrow{W_{R}})\end{aligned}\end{split}\]

ce qui fournit la relation de saut pour une petite perturbation

\[\mathcal{A}(\overrightarrow{W_{L}}-\overrightarrow{W_{R}})=s\left[\overrightarrow{W_{L}}-\overrightarrow{W_{R}}\right]\]

Clairement cette relation montre que le saut \(\left[\overrightarrow{W_{L}}-\overrightarrow{W_{R}}\right]\) doit être proportionnel à un vecteur propre \(\overrightarrow{r_{i}}\) de \(\mathcal{A}\), et que la vitesse de propagation \(s\) est la valeur propre \(\lambda_{i}\) de \(\mathcal{A}\) associée. Les valeurs possibles de \(s\) sont donc \([u-c,u,u+c]\).

Les relations de saut d’Hugoniot (4.6) s’écrivent:

(4.7)\[\begin{split}\begin{aligned} \rho_{l}u_{l}-\rho_{r}u_{r} & = & s\,(\rho_{l}-\rho_{r})\nonumber \\ \rho_{l}u_{l}^{2}+p_{l}-\rho_{r}u_{r}^{2}-p_{r} & = & s\,(\rho_{l}u_{l}-\rho_{r}u_{r})\\ u_{l}(\rho_{l}et_{l}+p_{l})-u_{r}(\rho_{r}et_{r}+p_{r}) & = & s\,(\rho_{l}et_{l}-\rho_{r}et_{r})\nonumber \end{aligned}\end{split}\]

Ces relations traduisent le bilan de la masse, la quantité de mouvement et de l’énergie dans un repère lié au choc, i.e. se déplaçant à la vitesse \(s\) de propagation de la perturbation. Pour \(s=0\), on retrouve les équations du choc droit stationnaire, qui correspond à la valeur propre \(u-c=0\) (i.e. une vitesse sonique dans le choc).

Il existe 3 types de discontinuité possible dans ce système d’équations d’Euler:

1. discontinuité de contact: qui est associée à la valeur propre \(s=u\). Dans ce cas le système (4.7) conduit à \(u_{l}=u_{r}\)et \(p_{r}=p_{l}\), i.e. il n’y a pas de discontinuité sur la vitesse et la pression. La seule discontinuité possible est une discontinuité sur la masse volumique \(\rho_{l}\neq\rho_{r}\). Dans ce cas le fluide à gauche ne se mélange pas avec le fluide à droite (puisqu’on a pas pris en compte de mécanisme de diffusion), et la séparation entre les deux se propage à la vitesse du fluide.

Fig. 4.28 discontinuité de contact

2. choc de compression: qui est associé aux valeurs propres \(s=u+c\) ou \(s=u-c\). Dans ce cas on a une discontinuité de vitesse, de pression et de masse volumique. D’après le second principe de la thermodynamique, seul un choc de compression est possible: i.e la pression du fluide augmente à travers un choc. Le cas d’un choc de détente est théoriquement possible dans le système (4.7), mais il viole le second principe de la thermodynamique (i.e. dans ce cas l’entropie \(S\) diminue).

Fig. 4.29 choc droit de compression

3. ondes de détente (ou de raréfaction): qui sont associées aux valeurs propres \(s=u+c\) ou \(s=u-c\). Le choc de détente étant impossible physiquement, la détente (i.e. une diminution de pression) ne s’opère pas à travers une discontinuité (choc), mais à travers une série d’ondes divergentes qui assure la continuité de la vitesse, pression et masse volumique (on une famille de courbes caractéristiques divergentes).

Fig. 4.30 onde de détente

Dans ce cas les caractéristiques issues de la discontinuité en \(x=0\) sont les courbes \(x=st\) (\(s\) variant de \(s_{l}\) à \(s_{r}\)). Dans cette zone, on cherche donc une solution de similarité sous la forme \(\overrightarrow{W}=\overrightarrow{W}(\xi)\) avec \(\xi=\frac{x}{t}\). Cette solution vérifie les équations d’Euler:

\[\frac{\partial\overrightarrow{W}}{\partial t}+\frac{\partial\overrightarrow{F}(\overrightarrow{W})}{\partial x}=0\]

qui après le changement de variable en \(\xi\) s’écrivent:

\[\xi\frac{d\overrightarrow{W}}{d\xi}=\frac{\partial\overrightarrow{F}}{\partial\overrightarrow{W}}\frac{d\overrightarrow{W}}{d\xi}\]

ce qui traduit le fait que \(\frac{d\overrightarrow{W}}{d\xi}\) est proportionnel à un vecteur propre de la matrice jacobienne des flux \(\mathcal{A}\), et que \(\xi\) est une valeur propre de \(\mathcal{A}\).

Pour les équations d’Euler en 2D où 3D, il existe un autre type de discontinuité lié à une discontinuité de la vitesse dans la direction perpendiculaire à l’écoulement.

1. ligne de glissement: c’est la discontinuité de la vitesse, que l’on observe sur les parois où le fluide possède une vitesse de glissement par rapport à la paroi (car on a négligé les effets de viscosité).

4.2. Tube à choc

On considère un tube de longueur \(L\), séparé en son milieu par une membrane avec d’un coté un gaz à haute pression (\(p_{0},\rho_{0}\)) et de l’autre un gaz à basse pression (\(p_{1},\rho_{1}\)).

Fig. 4.31 tube à choc

On enlève la membrane à l’instant \(t=0\). On introduit donc une discontinuité de pression, masse volumique et température dans le tube.

Fig. 4.32 solution initiale dans le tube à choc

Due à la différence de pression, le gaz de la chambre haute pression va se déplacer dans la chambre basse pression. Une zone entre les 2 gaz se met en mouvement avec une vitesse \(u_{3}>0\) et une pression \(p_{3}\): \(p_{1}>p_{3}>p_{2}\) avec en amont la propagation d’une onde de choc avec une célérité \(u+c=\frac{1}{2}(u_{3}+c_{3}+c_{2})\). En arrière de cette zone se développent des ondes de détente de pente \(u-c\): \(-c_{1}<u-c<u_{3}-c_{3}\). Enfin, si on néglige la diffusion, les deux gaz ne se mélangent pas, et la séparation entre les deux correspond à une discontinuité de contact qui se propage avec la vitesse \(u=u_{3}\).

On retrouve les 3 ondes pouvant se développer dans les équations d’Euler:

1. une onde de choc de célérité \(u+c\)

2. des ondes de détente de célérité \(u-c\)

3. une discontinuité de contact de célérité \(u\)

La solution des équations d’Euler pour ce problème est donnée sur les figures suivantes à l’instant \(t=\frac{L}{4c_{1}}\), pour un rapport de pression et de masse volumique de \(3\).

Fig. 4.33 solution à \(t=\frac{L}{4c_{1}}\) dans le tube à choc

Sur l’animation suivante, on a tracé l’évolution de la trajectoire des particules fluides dans un tube à choc.

4.3. Méthodes numériques

Avant d’étudier les méthodes numériques pour les équations d’Euler, nous ferons quelques rappels sur la simulation numérique de problème d’advection linéaire.

4.3.1. Problème linéaire

Considérons l’équation d’advection linéaire suivante:

(4.8)\[\frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x}=0\]

Pour calculer une solution approchée \(u_{i}^{n}\) de la solution exacte \(u(x,t)\), on introduit un maillage en espace de pas \(dx\) et une discrétisation en temps de pas \(dt\). Classiquement, on calcule la solution approchée \(u_{i}^{n}\) par discrétisation différences finies de l’équation exacte. La solution approchée \(u_{i}^{n}\) est alors interprétée comme l’approximation de la solution exacte au point de maillage \(x=idx\), \(t=ndt\):

\[u_{i}^{n}\approx u(idx,ndt)\]

On peut cependant donner une autre interprétation de \(u_{i}^{n}\) (dans un cadre volumes finis). La solution approchée est une approximation de la valeur moyenne de la solution exacte dans la cellule entourant le noeud \(i\):

\[u_{i}^{n}\approx\frac{1}{dx}\int_{x_{i}-\frac{dx}{2}}^{x_{i}+\frac{dx}{2}}u(x,t=ndt)dx\]

C’est cette interprétation que nous utiliserons lors de l’application des schémas à des équations de bilan (équation d’Euler), puisqu’elle va permettre de caractériser les propriétés de conservation du schéma numérique.

Pour étudier la convergence de la solution numérique \(u_{i}^{n}\) du schéma aux différences finies vers la solution exacte \(u(x,t)\) de l’équation, on utilise un résultat d’analyse du à Lax (Richtmyer et Norton 1967):

Théorème de Lax

pour un problème linéaire aux valeurs initiales, la solution numérique d’un schéma itératif en temps aux différences finies converge vers la solution exacte si le schéma est consistant et stable.

4.3.1.1. Schéma explicite centré

Le schéma le plus simple est le schéma explicite centré, qui utilise une discrétisation centrée de la dérivée spatiale. Il s’écrit :

(4.9)\[\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+V\,\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2dx}=0\]

C’est un schéma explicite qui calcule la valeur inconnue \(u_{i}^{n+1}\) à l’étape \(t^{n+1}\) en fonction des valeurs connues \(\{u_{i}^{n},u_{i+1}^{n},u_{i-1}^{n}\}\) à l’étape \(t^{n}\).

Pour étudier la stabilité, on calcule le facteur d’amplification \(G\) du schéma. En introduisant le paramètre sans dimension \(C_{FL}\) ou nombre de Courant Fredrich Lewis, dit aussi nombre de Courant:

(4.10)\[C_{FL}=\frac{V\,dt}{dx}\]

le carré du module du facteur d’amplification \(G\) du schéma vérifie:

\[\left|G\right|^{2}=1+C_{FL}^{2}\sin^{2}y\,>1\,\,\,\forall y\]

Le module de \(G\) est donc toujours plus grand que 1 et le schéma est inconditionnellement instable.

Avertissement

Le schéma centré explicite (4.9) est donc inutilisable pour résoudre l’équation de convection pure (4.8).

Pour comprendre la raison de cette instabilité, nous allons faire l’étude de la consistance de ce schéma. L’erreur de troncature du schéma est la suivante:

\[ErrT=\frac{1}{2}\frac{\partial^{2}u}{\partial t^{2}}\,dt+\frac{1}{6}V\,\frac{\partial^{3}u}{\partial x^{3}}\,dx^{2}+O(dt^{2},dx^{3})\]

Le schéma est donc bien consistant à l’équation (4.8), mais instable. Pour comprendre cette instabilité, étudions le premier terme de l’erreur:

(4.11)\[\frac{1}{2}\frac{\partial^{2}u}{\partial t^{2}}\,dt\]

En dérivant l’équation (4.8) par rapport au temps, nous pouvons en déduire (en utilisant à nouveau l’équation (4.8) que la solution \(u(x,t)\) vérifie:

(4.12)\[\frac{\partial^{2}u}{\partial t^{2}}=-V\frac{\partial^{2}u}{\partial t\partial x}\,=V^{2}\frac{\partial^{2}u}{\partial x^{2}}\,\]

d’où l’expression du premier terme (4.11) de l’erreur de troncature:

(4.13)\[\frac{V^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}}\,dt\]

Par définition de l’erreur de troncature, qui est l’écart entre l’équation approchée \(Eq^{h}\) et l’équation exacte \(Eq\), calculée pour la solution exacte \(u(x,t)\), on a donc:

\[Eq^{h}-Eq=\frac{V^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}}\,dt+O(dx^{2},dt^{2})\]

Avec ce schéma explicite centré, on ne résout donc pas l’équation exacte, mais l’équation équivalente suivante (en notant \(\kappa=-v^{2}\frac{dt}{2}\)):

(4.14)\[Eq^{h}\approx\frac{\partial u}{\partial t}+V\,\frac{\partial u}{\partial x}-\kappa\frac{\partial^{2}u}{\partial x^{2}}=0\]

C’est une équation de convection-diffusion classique, mais avec un coefficient de diffusion négatif \((\kappa<0)\), qui au lieu de diffuser l’énergie de la solution, apporte indéfiniment de l’énergie au système. En effet, en décomposant la solution initiale en série de Fourier:

\[u_{0}(x)=\sum_{k}C_{k}e^{I\omega_{k}x}\]

on vérifie que la solution de l’équation (4.14)s’écrit:

(4.15)\[u(x,t)=\sum_{k}C_{k}e^{I\omega_{k}(x-V\,t)}e^{-\kappa\omega^{2}t}\]

Pour \(\kappa=0\), on retrouve la solution générale de l’équation de convection pure (4.8), qui est convectée sans déformation. Si \(\kappa>0\), la solution est convectée mais décroît exponentiellement (problème de convection diffusion classique). Par contre si \(\kappa<0\) la solution est convectée mais croit exponentiellement. Ce dernier cas correspond à la solution numérique du schéma explicite centré (4.9), qui croît exponentiellement et est instable.

4.3.1.2. Schéma explicite d’ordre 2 avec diffusion numérique

Le schéma explicite (4.9) étant instable parce qu’insuffisamment diffusif, on peut le rendre stable en ajoutant un terme de diffusion numérique qui va compenser le premier terme de l’erreur de troncature (4.13), responsable de l’instabilité. Pour cela, il suffit de retrancher à l’équation (4.9) la discrétisation du premier terme de l’erreur de troncature (4.13):

\[\frac{V^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}}\,dt\approx\frac{V^{2}dt}{2}\,\frac{u_{i-1}^{n}-2u_{i}^{n}+u_{i+1}^{n}}{dx^{2}}\]

Le schéma ainsi obtenu s’écrit:

(4.16)\[\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+V\,\frac{u_{i+1}^{n}-u_{i-1}^{n}}{2dx}-\frac{V^{2}dt}{2}\,\frac{u_{i-1}^{n}-2u_{i}^{n}+u_{i+1}^{n}}{dx^{2}}=0\]

C’est un schéma de Lax-Wendroff. Ce schéma est à la base de nombreux schémas utilisés en mécanique des fluides pour traiter les problèmes à convection dominante.

C’est un schéma explicite qui calcule la valeur inconnue \(u_{i}^{n+1}\) à l’étape \(t^{n+1}\) en fonction des valeurs connues \(\{u_{i}^{n},u_{i+1}^{n},u_{i-1}^{n}\}\) à l’étape \(t^{n}\)

L’analyse de stabilité fournit le carré du module du facteur d’amplification \(G\) du schéma, en posant \(z=\cos y\):

\[\left|G\right|^{2}=(C_{FL}^{4}-C_{FL}^{2})(z-1)^{2}+1\]

C’est un polynôme de degré 2 en \(z\), que l’on étudie pour \(z\in[-1,1]\). En \(z=1\), ce polynôme vaut \(1\) et sa dérivée s’annule. La condition de stabilité \(\left|G\right|^{2}<1\,\,\forall z\in[-1,1]\) impose donc que sa valeur en \(z=-1\) soit plus petite que \(1\):

\[\left|G\right|_{z=-1}^{2}=4(C_{FL}^{4}-C_{FL}^{2})+1=(1-2C_{FL}^{2})^{2}<1\]

ce qui impose la condition de stabilité classique sur le nombre de Courant:

(4.17)\[C_{FL}<1\]

L’erreur de troncature, qui s’écrit, en tenant compte de la propriété (4.12) :

(4.18)\[ErrT=\left(\frac{1}{6}\frac{\partial^{3}u}{\partial t^{3}}\right)dt^{2}+\frac{1}{6}V\,\frac{\partial^{3}u}{\partial x^{3}}\,dx^{2}+O(dt^{3},dx^{3})\]

Astuce

Le schéma de Lax-Wendroff (4.16) est donc consistant à l’équation de convection pure (4.8), d’ordre \(O(dt^{2},dx^{2})\) et conditionnellement stable, avec une condition de stabilité donnée par (4.17).

On le vérifie numériquement en calculant la convection d’une tache gaussienne

\[u_{0}(x)=e^{-\left(\frac{x-x_{0}}{\sigma}\right)^{2}}\]

convectée par un champ de vitesse constant \(V=1\). La solution exacte est \(u(x,t)=u_{0}(x-V\,t)\).

Fig. 4.34 solution numérique Lax Wendroff

Fig. 4.35 Erreur solution Lax Wendroff

Sur la figure Fig. 4.34, on a tracé la solution numérique pour différents \(C_{FL}\) à \(t=0.5\). On constate la très bonne coïncidence entre la solution numérique et la solution exacte pour tous les \(C_{FL}\leq1\). Pour faire une analyse plus fine, on a tracé sur la figure Fig. 4.35 l’erreur entre la solution exacte et la solution numérique. On constate que l’erreur est très faible (\(<6\)%). Cette erreur croît lorsque le \(C_{FL}\) décroît pour atteindre une limite (obtenue dès \(C_{FL}=0.1\) puisque les courbes pour \(C_{FL}=0.1\) et \(C_{FL}=0.001\) sont indiscernables). Cette erreur est liée à l’erreur de discrétisation spatiale, et décroît lorsque l’on augmente \(N\).

Pour \(C_{FL}=1\), l’erreur est nulle, et on obtient dans ce cas la solution exacte. En effet l’équation (4.16) conduit pour \(C_{FL}=1\) à:

\[u_{i}^{n+1}=u_{i-1}^{n}\]

qui est la solution exacte car les deux points \((i,n+1)\) et \((i,n)\) se trouvent sur la même trajectoire, puisque:

\[idx-V(n+1)dt=(i-1)dx-Vndt\]

4.3.1.3. Schéma explicite décentré d’ordre 1

Une autre façon de stabiliser le schéma explicite (4.9), est d’utiliser une discrétisation décentrée du terme convectif \(V\,\frac{\partial u}{\partial x}\). Ce terme représente la convection de \(u\) par le champ de vitesse \(V\), et si \(V>0\) la convection transporte du scalaire \(u\) de la gauche vers la droite (et inversement si \(V<0\)). D’où l’idée d’utiliser une discrétisation décentrée de ce terme qui, si \(V>0\), fait intervenir les valeurs de \(u\) aux noeuds amonts \(i-1\) et \(i\) :

\[V\,\frac{\partial u}{\partial x}\approx V\,\frac{u_{i}-u_{i-1}}{dx}\,\,\mbox{\, si\,}V>0\]

et, si \(V<0\) , les valeurs de \(u\) aux noeuds avals \(i+1\) et \(i\):

\[V\,\frac{\partial u}{\partial x}\approx V\,\frac{u_{i+1}-u_{i}}{dx}\,\,\mbox{\, si\,}V<0\]

Le schéma ainsi obtenu s’écrit pour \(V>0\):

(4.19)\[\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+V\,\frac{u_{i}^{n}-u_{i-1}^{n}}{dx}=0\]

C’est un schéma explicite qui donne la valeur inconnue \(u_{i}^{n+1}\) à l’étape \(t^{n+1}\) en fonction des valeurs connues \(\{u_{i}^{n},u_{i-1}^{n}\}\) à l’étape \(t^{n}\). On note que la valeur aval \(u_{i+1}^{n}\) n’intervient pas dans ce cas pour le calcul de \(u_{i}^{n+1}\).

L’analyse de stabilité fournit le carré du module du facteur d’amplification \(G\) du schéma, qui s’écrit:

\[\left|G\right|^{2}=2(C_{FL}-C_{FL}^{2})(z-1)+1\]

C’est une fonction affine de \(z\) que l’on étudie pour \(z\in[-1,1]\), qui vaut 1 en \(z=1\). La condition de stabilité \(\left|G\right|^{2}<1\,\,\forall z\in[-1,1]\) impose donc que sa valeur en \(z=-1\) soit plus petite que \(1\):

\[\left|G\right|_{z=-1}^{2}=(1-2C_{FL})^{2}<1\]

ce qui conduit à la condition de stabilité classique:

(4.20)\[C_{FL}<1\]

L’erreur de troncature s’écrit:

\[ErrT=\frac{1}{2}\frac{\partial^{2}u}{\partial t^{2}}dt-\frac{1}{2}V\,\frac{\partial^{2}u}{\partial x^{2}}\,dx+\frac{1}{6}\frac{\partial^{3}u}{\partial t^{3}}\,dt^{2}+\frac{1}{6}V\,\frac{\partial^{3}u}{\partial x^{3}}\,dx^{2}+O(dt^{3},dx^{3})\]

soit en tenant compte de la propriété (4.12) :

\[ErrT=\frac{1}{2}V^{2}\frac{\partial^{2}u}{\partial x^{2}}dt-\frac{1}{2}V\,\frac{\partial^{2}u}{\partial x^{2}}\,dx+\frac{1}{6}\frac{\partial^{3}u}{\partial t^{3}}\,dt^{2}+\frac{1}{6}V\,\frac{\partial^{3}u}{\partial x^{3}}\,dx^{2}+O(dt^{3},dx^{3})\]

Le terme d’ordre 1 de l’erreur de troncature est donc un terme de diffusion:

\[-\kappa\frac{\partial^{2}u}{\partial x^{2}}\]

Le schéma explicite décentré est donc équivalent à l’équation de convection diffusion:

\[Eq^{h}\approx\frac{\partial u}{\partial t}+V\,\frac{\partial u}{\partial x}-\kappa\frac{\partial^{2}u}{\partial x^{2}}=0\]

avec un coefficient de diffusion numérique \(\kappa\):

(4.21)\[\kappa=\frac{1}{2}(V\,dx-V^{2}dt)=\frac{V\,dx}{2}(1-C_{FL}\]

Indication

Le schéma explicite décentré est donc consistant à l’équation de convection pure et d’ordre \(O(dt,dx)\). Il est conditionnellement stable avec une condition de stabilité {eq} c5aeq21`.

Mais c’est un schéma diffusif. La solution numérique est convectée, mais décroît exponentiellement proportionnellement à un coefficient de « diffusion numérique » \(\kappa\) (4.21) si la condition de stabilité (4.20) est vérifiée. Si la condition de stabilité (4.20) n’est pas vérifiée, la solution numérique croît exponentiellement puisque ce coefficient \(\kappa\) devient négatif, ce qui explique pourquoi le schéma devient instable.

On le vérifie numériquement sur le problème de la convection de la gaussienne.

Fig. 4.36 solution numérique schéma décentré

Sur la figure Fig. 4.36, on a tracé la solution numérique pour différents \(C_{FL}\) à \(t=0.5\). On constate l’importante diffusion numérique de la solution, qui augmente lorsque le \(C_{FL}\)diminue pour atteindre une limite (correspondant à un coefficient de diffusion numérique \(\kappa=\frac{Vdx}{2}\)). On note aussi le cas particulier \(C_{FL}=1\), où la solution numérique coïncide avec la solution exacte, puisque dans ce cas particulier le coefficient de diffusion numérique \(\kappa\) est nul.

4.3.2. Propriétés de dispersion

Nous avons vu dans les paragraphes précédents que la faible diffusion numérique d’un schéma permet d’assurer que la solution numérique est convectée sans être trop amortie. Pour les équations hyperboliques (équation de convection pure), il faut aussi que la solution numérique soit convectée à la bonne vitesse avec très peu de déphasage par rapport à la solution exacte. Pour étudier cette dispersion numérique, nous allons considérer le cas de la convection d’une onde \(u_{0}(x)=e^{I\omega x}\). La solution exacte de l’équation de convection pure (4.8) s’écrit:

\[u(x,t)=u_{0}(x-V\,t)=e^{I\,\omega x}e^{-I\,\omega V\,t}\label{c5aeq22a}\]

On note que cette solution initiale correspond justement à la perturbation \(\epsilon_{i}^{n}\) dans la méthode de stabilité de Neumann. L’étude de stabilité fournit le coefficient d’amplification \(G\) du schéma, qui permet de calculer \(\epsilon_{i}^{n}\) en fonction de la perturbation initiale \(\epsilon_{i}^{0}\):

\[\epsilon_{i}^{n}=G\,\epsilon_{i}^{n-1}=(G)^{n}\epsilon_{i}^{0}\]

La solution numérique \(u_{i}^{n}\) pour la convection d’une onde peut donc s’écrire en fonction de la solution initiale \(u_{0}\) au noeud \(i\) du maillage:

\[u_{i}^{n}=(G)^{n}u_{0}(idx)=(G)^{n}e^{I\omega idx}\]

Le facteur d’amplification \(G\) est un nombre complexe d’amplitude \(A\) et de phase \(-\phi\):

\[G=Ae^{-I\phi}\]

et la solution numérique \(u_{i}^{n}\) s’écrit sous la forme:

(4.22)\[u_{i}^{n}=(A)^{n}e^{-In\phi}e^{I\omega idx}\]

En comparant cette expression à la solution exacte au noeud \(i\):

\[u(idx,ndt)=e^{-In\omega V\,dt}e^{I\omega idx}\]

on note que \(A\), qui est le module de \(G\), nous donne l’amortissement de l’onde et la différence \(\varphi=V\omega dt-\phi\) sa dispersion. En notant que \(\phi\) dépend de \(\omega dx\) , on pose \(y=\omega dx\) et on a \(V\omega dt=C_{FL}y\). La dispersion \(\varphi\) est la fonction de \(y\) suivante:

\[\varphi=C_{FL}y-\phi(y)\]

Pour un maillage de \(N+1\) points, les modes possibles de la solution numérique correspondent à des pulsations \(\omega=\frac{k\pi}{L}\) (\(k\) de 1 à \(N-1\)) et donc \(y\) varie de \(0\) à \(\frac{N\pi}{L}dx=\pi\). On trouve pour le schéma décentré et le schéma explicite d’ordre 2, en utilisant Maple pour le calcul de \(\varphi\):

1. schéma 1: décentré d’ordre 1 (4.19):

   \[\varphi=C_{FL}y-\arctan\frac{C_{FL}\sin y}{1+C_{FL}^{2}(\cos y-1)}\]

2. schéma 2: Lax-Wendroff explicite d’ordre 2 (4.16):

   \[\varphi=C_{FL}y-\arctan\frac{C_{FL}\sin y}{1+C_{FL}(\cos y-1)}\]

Fig. 4.37 dispersion des schémas

On a tracé sur la figure Fig. 4.37 les courbes de \(\varphi\) pour différentes valeurs du \(C_{FL}\).

On constate sur cette figure que le schéma décentré est moins dispersif que le schéma d’ordre 2. Pour valider cette analyse, nous avons effectuer la simulation numérique de la convection d’une onde.

Le champ initial est une onde \(\sin\frac{m\pi x}{L}\) (lignes 10 à 11). Les schémas étudiés sont des schémas explicites du type:

(4.23)\[u_{i}^{n+1}=c_{1}u_{i-1}^{n}+c_{2}u_{i}^{n}+c_{3}u_{i+1}^{n}\]

Le schéma numérique étant explicite, la solution s’obtient simplement à chaque itération par multiplication matricielle:

\[U^{n+1}=\mathcal{C}U^{n}\]

Pour ce calcul, il faut introduire des conditions aux limites. La longueur \(L=\pi\) du domaine est un multiple de la longueur d’onde de l’onde initiale \(\frac{L}{m\pi}\). La solution vérifie donc des conditions de périodicité en \(x=0\) et \(x=L\):

\[u(0,t)=u(L,t)\,\,\,,\,\,\,\frac{\partial u}{\partial x}(0,t)=\frac{\partial u}{\partial x}(L,t)\]

Pour le noeud \(i=0\), le schéma aux différences (4.23) fait intervenir la valeur \(u_{-1}^{n}\), que l’on calcule avec la condition de périodicité : \(u_{-1}^{n}=u_{N-1}^{n}\). Pour le noeud \(i=N\), on impose la condition \(u_{N}^{n}=u_{0}^{n}\).

On a tracé sur la figure Fig. 4.37 la solution obtenue avec le schéma décentré d’ordre 1 (4.19) et le schéma d’ordre 2 (4.16). On constate bien la dispersion du schéma d’ordre 2 (décalage des courbes), mais aussi la très forte diffusion du schéma décentré (mais qui a une dispersion plus faible).

Pour interprèter l’apparition de cette dispersion dans le schéma d’ordre 2 (4.16), on transforme l’erreur de troncature (4.18) en dérivant l’équation exacte pour remplacer \(\frac{\partial^{3}u}{\partial t^{3}}=-V^{3}\frac{\partial^{3}u}{\partial x^{3}}\) :

\[ErrT=\frac{V}{6}(dx^{2}-V^{2}dt^{2})\,\frac{\partial^{3}u}{\partial x^{3}}+O(dt^{3},dx^{3})\]

L’équation aux différences pour le schéma explicite d’ordre 2 est donc équivalente à l’équation de convection dispersion suivante:

\[Eq^{h}\approx\frac{\partial u}{\partial t}+V\,\frac{\partial u}{\partial x}+\varphi\frac{\partial^{3}u}{\partial x^{3}}=0\]

avec \(\varphi=\frac{Vdx^{2}}{6}(1-C_{FL}^{2})\). La solution exacte, pour une condition initiale \(u_{0}(x)=e^{I\omega x}\), s’écrit:

\[u(x,t)=e^{I\,\omega x}e^{-I\,\omega V\,t}e^{I\varphi\omega^{3}t}=e^{I\omega(x-(V-\varphi\omega^{2})t)}\]

L’onde initiale est donc transportée avec une vitesse de convection \(V-\varphi\omega^{2}\), qui dépend de la fréquence de l’onde. Les ondes sont donc dispersées (i.e. ne se propagent pas toutes à la même vitesse comme dans l’équation exacte (4.8), avec une dispersion d’autant plus grande que l’onde est à haute fréquence. On note aussi que, pour le cas particulier \(C_{FL}=1\), le schéma (4.16) n’est plus dispersif, ce que l’on vérifie bien numériquement.

4.3.3. Convection d’une discontinuité

La capture des chocs nécessite des schémas numériques qui propagent correctement les discontinuités. Nous allons donc étudier numériquement le problème d’advection suivant:

\[\begin{split}\begin{aligned} \frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x} & = & 0\\ u(x,t=0) & = & u_{0}(x)\end{aligned}\end{split}\]

où la condition initiale \(u_{0}(x)\) est discontinue en \(x=x_{0}\)

\[\begin{split}u_{0}(x)=\left\{ \begin{array}{ll} 1 & x\leq x_{0}\\ 0 & x>x_{0} \end{array}\right.\end{split}\]

dont la solution exacte est

\[u(x,t)=u_{0}(x-V\,t)\]

Fig. 4.38 schéma d’ordre 1: convection d’une discontinuité avec \(N=101\)

Fig. 4.39 schéma d’ordre 1: convection d’une discontinuité avec \(N=1001\)

Sur la figure Fig. 4.38, nous avons comparé la solution exacte avec la solution approchée calculée sur un maillage de \(101\) points (\(dx=0.01\)).

On constate que la solution numérique avec le schéma décentré d’ordre 1 (upwind) présente une forte diffusion, et la discontinuité se trouve étalée sur une vingtaine de mailles. Avec le schéma d’ordre 2 de Lax Wendroff, la discontinuité est mieux captée (i.e. l’étalement du front est inférieure à 10 mailles). Si on raffine le maillage, les mêmes conclusions peuvent être tirées, comme le montre la figure Fig. 4.39 où on a tracé la solution numérique sur un maillage 10 fois plus fin \((dx=0.001)\).

La convergence de la solution numérique est dégradée au voisinage de la discontinuité. La méthode décentrée d’ordre 1 converge en \(O(\sqrt{dx})\) au lieu de \(O(dx),\)et la méthode d’ordre 2 converge en \(O(dx^{2/3})\) au lieu de \(O(dx^{2})\).

Fig. 4.40 schéma d’ordre 2: convection d’une discontinuité avec \(N=101\)

Fig. 4.41 schéma d’ordre 2: convection d’une discontinuité avec \(N=1001\)

Le comportement de la solution décentrée s’explique par la forte dissipation numérique du schéma upwind. Celui de la solution d’ordre 2 est lié à la dispersion numérique du schéma de Lax Wendroff, qui ralentit les ondes hautes fréquences. Or la convection d’une discontinuité correspond à la convection d’une infinité d’ondes basses et hautes fréquences (décomposition en série de Fourier), qui se propagent toutes à la même vitesse. Le schéma Lax Wendroff ralentit les ondes hautes fréquences, et donc la discontinuité n’est plus reconstituée exactement, et une oscillation haute fréquence apparaît à l’arrière de la discontinuité comme le montrent les figures Fig. 4.40 et Fig. 4.41 (la fréquence de l’oscillation augmente avec le nombre de points du maillage).

4.3.4. Problème non linéaire

En plus des difficultés numériques exposées précédemment pour une équation linéaire d’advection, la résolution numérique de problème non-linéaire se heurte à de nouvelles difficultés:

1. la solution numérique peut devenir « non-linéairement instable », i.e. même en respectant les conditions de stabilité du schéma linéarisé, la solution non-linéaire devient instable.

2. la solution numérique peut converger vers une solution qui n’est pas solution du problème exact, ou qui est physiquement inacceptable (i.e. qui viole la condition d’entropie).

A titre d’exemple, considérons l’équation non linéaire de Burgers

(4.24)\[\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0\]

avec une condition initiale \(u(x,t=0)=u_{0}(x)\) discontinue

(4.25)\[\begin{split}u_{0}(0)=\left\{ \begin{array}{ll} u_{l} & x\leq x_{0}\\ u_{r} & x>x_{0} \end{array}\right.\end{split}\]

La solution exacte est un choc qui se propage avec la célérité \(s\), et qui vérifie la relation de saut d’Hugoniot

(4.26)\[F(u_{l})-F(u_{r})=s\,[u_{l}-u_{r}]\]

où \(F(u)=\frac{1}{2}u^{2}\)est la fonction flux de l’équation de Burgers (4.24), qui s’écrit sous forme conservative:

(4.27)\[\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}(\frac{1}{2}u^{2})=0\]

La relation de saut (4.26) fournit la valeur de la célérité \(s\)

\[s=\frac{1}{2}(u_{l}+u_{r})\]

L’extension naturelle du schéma aux différences finies décentrées (4.19) à l’équation de Burgers (4.24) peut s’écrire:

(4.28)\[\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+u_{i-1}^{n}\frac{u_{i}^{n}-u_{i-1}^{n}}{dx}=0\]

en utilisant une vitesse de convection décentrée \(u_{i-1}^{n}\), ou

(4.29)\[\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+u_{i}^{n}\frac{u_{i}^{n}-u_{i-1}^{n}}{dx}=0\]

en utilisant une vitesse de convection centrée \(u_{i}^{n}\).

Fig. 4.42 solution numérique avec les 2 schémas décentrés (4.28) et (4.29) pour \(N=101\)

Fig. 4.43 solution numérique avec les 2 schémas décentrés (4.28) et (4.29) pour \(N=1001\)

Sur les figures Fig. 4.42 et Fig. 4.43 , on a comparé la solution exacte avec la solution numérique calculée avec les 2 schémas décentrés (4.28) et (4.29) pour deux maillages de \(N=101\) et \(N=1001\) points. On constate que les deux solutions numériques ne convergent pas vers la solution exacte, et donnent une célérité \(s\) 20% trop petite avec le schéma upwind1 (4.28) et 20% trop grande avec le schéma upwind2 (4.29)

Pour éviter ce comportement, il faut imposer une condition supplémentaire au schéma numérique. Il faut que le schéma soit conservatif, i.e. soit écrit sous une forme conservative.

On part de l’équation exacte sous forme conservative:

(4.30)\[\frac{\partial u}{\partial t}+\frac{\partial F(u)}{\partial x}=0\]

Un schéma conservatif par différences finies s’écrit:

(4.31)\[u_{i}^{n+1}=u_{i}^{n}-\frac{dt}{dx}\left[F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}}\right]\]

où \(F_{i+\frac{1}{2}}\) est un flux numérique en \(x=x_{i}+\frac{dx}{2}\). Cette discrétisation correspond à l’approximation par volumes finis de (4.30) obtenue en intégrant l’équation en espace sur la cellule \([x_{i}-\frac{dx}{2},x_{i}+\frac{dx}{2}]\) et en temps entre \(t^{n}\) et \(t^{n+1}\):

\[\begin{split}\begin{aligned} \int_{x_{i-1/2}}^{x_{i+1/2}}u(x,t^{n+1})dx & = & \int_{x_{i-1/2}}^{x_{i+1/2}}u(x,t^{n})dx\\ & - & \left[\int_{t^{n}}^{t^{n+1}}F(u(x_{i+1/2},t))dt-\int_{t^{n}}^{t^{n+1}}F(u(x_{i-1/2},t))dt\right]\end{aligned}\end{split}\]

La valeur approchée \(u_{i}^{n}\) est donc une approximation de la valeur moyenne sur la cellule:

\[u_{i}^{n}\approx\frac{1}{dx}\int_{x_{i-1/2}}^{x_{i+1/2}}u(x,t^{n})dx\]

et la valeur du flux numérique \(F_{i+\frac{1}{2}}\) une valeur approchée du flux physique entre \(t^{n}\) et \(t^{n+1}\):

\[F_{i+\frac{1}{2}}=\frac{1}{dt}\int_{t^{n}}^{t^{n+1}}F(u(x_{i+1/2},t)dt\]

Ainsi le schéma explicite décentré conservatif pour l’équation de Burgers s’écrit:

(4.32)\[u_{i}^{n+1}=u_{i}^{n}-\frac{dt}{dx}\left[\frac{1}{2}\left(u_{i}^{n}\right)^{2}-\frac{1}{2}\left(u_{i-1}^{n}\right)^{2}\right]\]

Fig. 4.44 solution numérique avec schéma conservatif (4.32) pour \(N=101\)

Fig. 4.45 solution numérique avec schéma conservatif (4.32) pour \(N=1001\)

Sur les figures Fig. 4.44 et Fig. 4.45, on a comparé la solution exacte de l’équation de Burgers (4.24) et (4.25) avec la solution numérique calculée avec le schéma conservatif décentré (4.32) pour deux maillages de \(N=101\) et \(N=1001\) points. On constate qu’avec ce schéma conservatif, la solution approchée converge vers la solution exacte et fournit la bonne vitesse de propagation \(s=\frac{u_{l}+u_{r}}{2}\).

4.3.5. Système d’équations hyperboliques

Soit le système de \(n\) équations hyperboliques écrit sous forme conservative:

\[\frac{\partial\overrightarrow{W}}{\partial t}+\frac{\partial\overrightarrow{F}(\overrightarrow{W})}{\partial x}=0\]

En notant \(\mathcal{A}=\frac{\partial\overrightarrow{F}}{\partial\overrightarrow{W}}\), la matrice jacobienne des flux, ce système s’écrit sous la forme non conservative:

\[\frac{\partial\overrightarrow{W}}{\partial t}+\mathcal{A}\frac{\partial\overrightarrow{W}}{\partial x}=0\]

En diagonalisant la matrice \(\mathcal{A}\), on obtient un système de n équations d’advection dont les vitesses d’advection sont les valeurs propres \(\lambda_{i}\) de \(\mathcal{A}\). Pour construire un schéma numérique décentré, il faut donc tenir compte du signe des valeurs propres de \(\mathcal{A}\) lors du calcul des flux à l’interface.

Considérons tout d’abord le cas linéaire (i.e. \(\mathcal{A}\) est indépendant de \(W\)). En notant \(\Lambda\) la matrice diagonale des valeurs propres \(\lambda_{i}\), la diagonalisation de \(\mathcal{A}\) s’écrit:

\[\mathcal{A}=\mathcal{R}\Lambda\mathcal{R}^{-1}\]

Soit \(\Lambda^{+}=max(\lambda_{i},0)\) la matrice diagonale des valeurs propres positives de \(\mathcal{A}\) et \(\Lambda^{-}=min(\lambda_{i},0)\) la matrice diagonale des valeurs propres négatives. On note

\[\mathcal{A}^{+}=\mathcal{R}\Lambda^{+}\mathcal{R}^{-1}\,\mbox{ et }\,\,\mathcal{A}^{-}=\mathcal{R}\Lambda^{-}\mathcal{R}^{-1}\]

qui sont deux matrices vérifiant:

(4.33)\[\mathcal{A}=\mathcal{A}^{+}+\mathcal{A}^{-}\,\mbox{ et }\left|\mathcal{A}\right|=\mathcal{A}^{+}-\mathcal{A}^{-}\]

Pour calculer le flux \(F(W)=\mathcal{A}W\) à l’interface \(x_{i+1/2}\), on peut utiliser un flux numérique centré:

\[F_{i+1/2}^{c}=\frac{1}{2}\left[F(W_{i})+F(W_{i+1})\right]=\frac{1}{2}\mathcal{A}(W_{i}+W_{i+1})\]

Ce flux numérique conduit à un schéma centré, qui est instable. Pour stabiliser le schéma, on décentre le flux suivant le signe des valeurs propres de \(\mathcal{A}\):

(4.34)\[F_{i+1/2}^{d}=F^{+}(W_{i})+F^{-}(W_{i+1})=\mathcal{A}^{+}W_{i}+\mathcal{A}^{-}W_{i+1}\]

Le schéma conservatif explicite décentré s’écrit alors:

(4.35)\[W_{i}^{n+1}=W_{i}^{n}-\frac{dt}{dx}\left[\mathcal{A}^{+}\left(W_{i}-W_{i-1}\right)+\mathcal{A}^{-}\left(W_{i+1}-W_{i}\right)\right]\]

En utilisant les propriétés (4.33) de \(\mathcal{A}\), le flux décentré s’écrit aussi sous la forme:

\[\begin{split}\begin{aligned} F_{i+1/2}^{d} & = & \frac{1}{2}\mathcal{A}\left(W_{i+1}+W_{i}\right)-\frac{1}{2}\left(\mathcal{A}^{+}-\mathcal{A}^{-}\right)\left(W_{i+1}-W_{i}\right)\\ & = & F_{i+1/2}^{c}-\frac{1}{2}\mathcal{\left|A\right|}\left(W_{i+1}-W_{i}\right)\end{aligned}\end{split}\]

et le schéma décentré (4.35) s’écrit:

\[W_{i}^{n+1}=W_{i}^{n}-\frac{dt}{2dx}\mathcal{A}\left(W_{i+1}^{n}-W_{i-1}^{n}\right)+\frac{dt}{2dx}\left|\mathcal{A}\right|\left(W_{i+1}^{n}-2W_{i}^{n}+W_{i-1}^{n}\right)\]

ce qui correspond à un schéma centré avec un terme supplémentaire de dissipation qui stabilise le schéma centré.

La condition de stabilité de ce schéma est une condition de type \(C_{FL}\) basée sur les vitesses d’advection du système, i.e. sur les valeurs propres de \(\mathcal{A}\):

\[C_{FL}=\frac{Max(|\lambda_{i}|)\,dt}{dx}<1\]

Pour un problème non-linéaire, la matrice \(\mathcal{A}\) dépend de la solution \(W\). L’analyse précédente doit être modifiée. Pour cela on considère un état intermédiaire \(\widehat{W}\) pour calculer la valeur \(\mathcal{\widehat{A}}\) de \(\mathcal{A}\) intervenant dans le calcul du flux décentré (4.34). Pour obtenir un schéma conservatif, la matrice \(\mathcal{\widehat{A}}\) doit vérifier:

\[\widehat{\mathcal{A}}\left(W_{i+1}-W_{i}\right)=F(W_{i+1})-F(W_{i})\]

et doit tendre vers la jacobienne exacte \(\frac{\partial F}{\partial W}\). En effet le problème linéarisé s’écrit:

\[\frac{\partial\overrightarrow{W}}{\partial t}+\mathcal{\widehat{A}}\frac{\partial\overrightarrow{W}}{\partial x}=0\]

et correspond à un flux linéarisé \(\widehat{F}(W)=\widehat{\mathcal{A}}W\). Pour être conservatif, la solution doit vérifier les relations de saut (4.5) et (4.6):

\[\begin{split}\begin{aligned} \frac{d}{dt}\int_{-x_{0}}^{+x_{0}}\overrightarrow{W}\,dx & = & \overrightarrow{F}(\overrightarrow{W_{L}})-\overrightarrow{F}(\overrightarrow{W_{R}})\\ \widehat{F}(\overrightarrow{W_{L}})-\widehat{F}(\overrightarrow{W_{R}}) & = & \overrightarrow{F}(\overrightarrow{W_{L}})-\overrightarrow{F}(\overrightarrow{W_{R}})\\ \widehat{\mathcal{A}}\left(\overrightarrow{W_{L}}-\overrightarrow{W_{R}}\right) & = & \overrightarrow{F}(\overrightarrow{W_{L}})-\overrightarrow{F}(\overrightarrow{W_{R}})\end{aligned}\end{split}\]

4.3.6. Schéma de Roe

Pour les équations d’Euler, on peut utiliser pour \(\widehat{W}\), l’état de Roe, qui est la moyenne de Roe entre \(W_{i}\) et \(W_{i+1}\)

\[\begin{split}\widehat{W}=\left[\begin{array}{l} \sqrt{\rho_{i}\rho_{i+1}}\\ \frac{\sqrt{\rho_{i}}u_{i}+\sqrt{\rho_{i+1}}u_{i+1}}{\sqrt{\rho_{i}}+\sqrt{\rho_{i+1}}}\\ \frac{\sqrt{\rho_{i}}ht_{i}+\sqrt{\rho_{i+1}}ht_{i+1}}{\sqrt{\rho_{i}}+\sqrt{\rho_{i+1}}} \end{array}\right]\end{split}\]

La matrice Jacobienne des flux de Roe s’écrit:

\[\begin{split}\mathcal{\widehat{A}}=\left[\begin{array}{ccc} 0 & 1 & 0\\ \frac{\gamma-3}{2}\widehat{u}^{2} & (3-\gamma)\widehat{u} & \gamma-1\\ \widehat{u}\,(\frac{\gamma-2}{2}\widehat{u}^{2}-\frac{1}{\gamma-1}\widehat{c}^{2}) & \frac{3-2\gamma}{2}\widehat{u}^{2}+\frac{1}{\gamma-1}\widehat{c}^{2} & \widehat{u}\gamma \end{array}\right]\end{split}\]

Sa décomposition en valeurs propres s’écrit:

\[\mathcal{\widehat{A}}=\mathcal{\widehat{R}}\widehat{\Lambda}\mathcal{\widehat{R}}^{-1}\]

avec

\[\begin{split}\Lambda=\left[\begin{array}{ccc} u & 0 & 0\\ 0 & u+c & 0\\ 0 & 0 & u-c \end{array}\right]\end{split}\]

\[\begin{split}\mathcal{R}=\left[\begin{array}{ccc} 1 & 1 & 1\\ u & u+c & u-c\\ \frac{1}{2}u^{2} & h+u+c & h-u+c \end{array}\right]\end{split}\]

\[\begin{split}\mathcal{R}^{-1}=\left[\begin{array}{ccc} 1-\frac{1}{2}r_{2} & r_{1} & -\frac{1}{v_{2}}\\ \frac{1}{4}r_{2}-\frac{1}{2}\frac{u}{c} & -\frac{1}{2}r_{1}+\frac{1}{2c} & \frac{1}{2v_{2}}\\ \frac{1}{4}r_{2}+\frac{1}{2}\frac{u}{c} & -\frac{1}{2}r_{1}-\frac{1}{2c} & \frac{1}{2v_{2}} \end{array}\right]\end{split}\]

avec \(h=\frac{1}{2}u^{2}+\frac{1}{\gamma-1}c^{2}\), \(v_{2}=\frac{1}{\gamma-1}c^{2}\), \(r_{1}=\frac{u}{v_{2}}\), \(r_{2}=\frac{u^{2}}{v_{2}}\).

Le flux de Roe s’écrit alors

\[F_{i+1/2}^{Roe}=\frac{1}{2}\left[F(W_{i+1})+F(W_{i})\right]-\frac{1}{2}\mathcal{\left|\widehat{A}\right|}\left(W_{i+1}-W_{i}\right)\]

Ce schéma de Roe est à la base de nombreux schémas numériques pour la simulation des écoulements compressibles.