5.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.

5.3.1 Problème linéaire

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

$$\frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x}=0$$ (5.7)

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.

5.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:

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

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:

$$C_{FL}=\frac{V\, dt}{dx}$$ (5.9)

le carré du module du facteur d'amplification $G$ du schéma vérifie:

$$\left\vert G\right\vert^{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.

Le schéma centré explicite (5.8) est donc inutilisable pour résoudre l'équation de convection pure (5.7).

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+\f... ...,\frac{\partial^{3}u}{\partial x^{3}}\, dx^{2}+O(dt^{2},dx^{3})$$

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

$$\frac{1}{2}\frac{\partial^{2}u}{\partial t^{2}}\, dt$$ (5.10)

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

$$\frac{\partial^{2}u}{\partial t^{2}}=-V\frac{\partial^{2}u}{... ...ial t\partial x}\,=V^{2}\frac{\partial^{2}u}{\partial x^{2}}\,$$ (5.11)

d'où l'expression du premier terme (5.10) de l'erreur de troncature:

$$\frac{V^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}}\, dt$$ (5.12)

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}$):

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

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 (5.13) s'écrit:

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

Pour $\kappa=0$, on retrouve la solution générale de l'équation de convection pure (5.7), 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é (5.8), qui croıt exponentiellement et est instable.

5.3.1.2 Schéma explicite d'ordre 2 avec diffusion numérique

Le schéma explicite (5.8) é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 (5.12), responsable de l'instabilité. Pour cela, il suffit de retrancher à l'équation (5.8) la discrétisation du premier terme de l'erreur de troncature (5.12):

$$\frac{V^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}}\, dt\app... ...{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:

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+V\,\frac{u_{i+1}^{n}-u_{i-1... ...{2}dt}{2}\,\frac{u_{i-1}^{n}-2u_{i}^{n}+u_{i+1}^{n}}{dx^{2}}=0$$ (5.15)

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\vert G\right\vert^{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\vert G\right\vert^{2}<1\,\,\forall z\in[-1,1]$ impose donc que sa valeur en $z=-1$ soit plus petite que $1$:

$$\left\vert G\right\vert _{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:

$$C_{FL}<1$$ (5.16)

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

$$ErrT=\left(\frac{1}{6}\frac{\partial^{3}u}{\partial t^{3}}\r... ...\frac{\partial^{3}u}{\partial x^{3}}\, dx^{2}+O(dt^{3},dx^{3})$$ (5.17)

Le schéma de Lax-Wendroff (5.15) est donc consistant à l'équation de convection pure (5.7), d'ordre $O(dt^{2},dx^{2})$ et conditionnellement stable, avec une condition de stabilité donnée par (5.16).

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)$.

Sur la figure (5.2a), 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 (5.2b) 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$.

Figure 5.2: Solution numérique avec le schéma explicite d'ordre 2 (5.15)

[solution numérique]

[erreur à $t=0.5$]

Pour $C_{FL}=1$, l'erreur est nulle, et on obtient dans ce cas la solution exacte. En effet l'équation (5.15) 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$$

5.3.1.3 Schéma explicite décentré d'ordre 1

Une autre façon de stabiliser le schéma explicite (5.8), 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$:

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+V\,\frac{u_{i}^{n}-u_{i-1}^{n}}{dx}=0$$ (5.18)

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\vert G\right\vert^{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\vert G\right\vert^{2}<1\,\,\forall z\in[-1,1]$ impose donc que sa valeur en $z=-1$ soit plus petite que $1$:

$$\left\vert G\right\vert _{z=-1}^{2}=(1-2C_{FL})^{2}<1$$

ce qui conduit à la condition de stabilité classique:

$$C_{FL}<1$$ (5.19)

L'erreur de troncature s'écrit:

$$ErrT=\frac{1}{2}\frac{\partial^{2}u}{\partial t^{2}}dt-\frac... ...,\frac{\partial^{3}u}{\partial x^{3}}\, dx^{2}+O(dt^{3},dx^{3})$$

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

$$ErrT=\frac{1}{2}V^{2}\frac{\partial^{2}u}{\partial x^{2}}dt-... ...,\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$:

$$\kappa=\frac{1}{2}(V\, dx-V^{2}dt)=\frac{V\, dx}{2}(1-C_{FL})$$ (5.20)

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é (5.19).

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$ (5.20) si la condition de stabilité (5.19) est vérifiée. Si la condition de stabilité (5.19) 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.

Sur la figure (5.3), 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.

Figure 5.3: Solution numérique avec le schéma explicite décentré (5.18)

5.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 (5.7) s'écrit:

$$u(x,t)=u_{0}(x-V\, t)=e^{I\,\omega x}e^{-I\,\omega V\, t}$$ (5.21)

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:

$$u_{i}^{n}=(A)^{n}e^{-In\phi}e^{I\omega idx}$$ (5.22)

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 (5.18):
   
   $$\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 (5.15):
   
   
   $$\varphi=C_{FL}y-\arctan\frac{C_{FL}\sin y}{1+C_{FL}(\cos y-1)}$$
   
   

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

Figure 5.4: Diffusion numérique du schéma d'ordre 1 (5.18) et d'ordre 2 (5.15)

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:

$$u_{i}^{n+1}=c_{1}u_{i-1}^{n}+c_{2}u_{i}^{n}+c_{3}u_{i+1}^{n}$$ (5.23)

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 (5.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 (5.5) la solution obtenue avec le schéma décentré d'ordre 1 (5.18) et le schéma d'ordre 2 (5.15). 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).

Figure 5.5: convection d'une onde ($C_{FL}=0.9)$ avec le schéma d'ordre 1 (5.18) et d'ordre 2 (5.15)

[schèma d'ordre 1]

[schèma d'ordre 2]

Pour interprèter l'apparition de cette dispersion dans le schéma d'ordre 2 (5.15), on transforme l'erreur de troncature (5.17) 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 (5.7)), 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 (5.15) n'est plus dispersif, ce que l'on vérifie bien numériquement.

5.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{eqnarray*} \frac{\partial u}{\partial t}+V\frac{\partial u}{\partial x} & = & 0\\ u(x,t=0) & = & u_{0}(x)\end{eqnarray*}$$

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

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

dont la solution exacte est

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

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

Figure 5.6: solutions numériques et exactes à $t=0.5$ pour $dx=0.01$ et $C_{FL}=0.6$

[Upwind]

[Lax Wendroff]

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 (5.7) où on a tracé la solution numérique sur un maillage 10 fois plus fin $(dx=0.001)$.

Figure 5.7: solutions numériques et exactes à $t=0.5$ pour $dx=0.001$ et $C_{FL}=0.6$

[Upwind]

[Lax Wendroff]

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})$.

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 (5.6) et (5.7) (la fréquence de l'oscillation augmente avec le nombre de points du maillage).

5.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

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0$$ (5.24)

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

$$u_{0}(0)=\left\{ \begin{array}{ll} u_{l} & x\leq x_{0}\\ u_{r} & x>x_{0}\end{array}\right.$$ (5.25)

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

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

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

$$\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}(\frac{1}{2}u^{2})=0$$ (5.27)

La relation de saut (5.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 (5.18) à l'équation de Burgers (5.24) peut s'écrire:

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+u_{i-1}^{n}\frac{u_{i}^{n}-u_{i-1}^{n}}{dx}=0$$ (5.28)

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

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}+u_{i}^{n}\frac{u_{i}^{n}-u_{i-1}^{n}}{dx}=0$$ (5.29)

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

Figure 5.8: Solutions exacte et approchée de l'équation de Burgers avec 2 schémas décentrés non conservatifs

[dx=0.01]

[dx=0.001]

Sur la figure (5.8), on a comparé la solution exacte avec la solution numérique calculée avec les 2 schémas décentrés (5.28) et (5.29) et 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 (5.28) et 20% trop grande avec le schéma upwind2 (5.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:

$$\frac{\partial u}{\partial t}+\frac{\partial F(u)}{\partial x}=0$$ (5.30)

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

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

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 (5.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{eqnarray*} \int_{x_{i-1/2}}^{x_{i+1/2}}u(x,t^{n+1})dx & = & \int_{x_{i-1/... ...(x_{i+1/2},t))dt-\int_{t^{n}}^{t^{n+1}}F(u(x_{i-1/2},t))dt\right]\end{eqnarray*}$$

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:

$$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]$$ (5.32)

Figure: Solution exacte et approchée de l'équation de Burgers avec le schéma conservatif (5.32)

[dx=0.01]

[dx=0.001]

Sur la figure (5.9), on a comparé la solution exacte de l'équation de Burgers (5.24,5.25) avec la solution numérique calculée avec le schéma conservatif décentré (5.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}$.

5.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}\,\mbo... ...t\,}}\,\,\mathcal{A}^{-}=\mathcal{R}\Lambda^{-}\mathcal{R}^{-1}$$

qui sont deux matrices vérifiant:

$$\mathcal{A}=\mathcal{A}^{+}+\mathcal{A}^{-}\,\mbox{{\, et\,}... ...eft\vert\mathcal{A}\right\vert=\mathcal{A}^{+}-\mathcal{A}^{-}$$ (5.33)

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}$:

$$F_{i+1/2}^{d}=F^{+}(W_{i})+F^{-}(W_{i+1})=\mathcal{A}^{+}W_{i}+\mathcal{A}^{-}W_{i+1}$$ (5.34)

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

$$W_{i}^{n+1}=W_{i}^{n}-\frac{dt}{dx}\left[\mathcal{A}^{+}\lef... ..._{i-1}\right)+\mathcal{A}^{-}\left(W_{i+1}-W_{i}\right)\right]$$ (5.35)

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

$$\begin{eqnarray*} F_{i+1/2}^{d} & = & \frac{1}{2}\mathcal{A}\left(W_{i+1}+W_{i}\... ...{1}{2}\mathcal{\left\vert A\right\vert}\left(W_{i+1}-W_{i}\right)\end{eqnarray*}$$

et le schéma décentré (5.35) s'écrit:

$$W_{i}^{n+1}=W_{i}^{n}-\frac{dt}{2dx}\mathcal{A}\left(W_{i+1}... ...al{A}\right\vert\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(\vert\lambda_{i}\vert)\, 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é (5.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 (5.4,5.5):

$$\begin{eqnarray*} \frac{d}{dt}\int_{-x_{0}}^{+x_{0}}\overrightarrow{W}\, dx & = ... ...overrightarrow{W_{L}})-\overrightarrow{F}(\overrightarrow{W_{R}})\end{eqnarray*}$$

5.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}$

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

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

$$\mathcal{\widehat{A}}=\left[\begin{array}{ccc} 0 & 1 & 0\\ ... ...{\gamma-1}\widehat{c}^{2} & \widehat{u}\gamma\end{array}\right]$$

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

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

avec

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

$$\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]$$

$$\mathcal{R}^{-1}=\left[\begin{array}{ccc} 1-\frac{1}{2}r_{2}... ...ac{1}{2}r_{1}-\frac{1}{2c} & \frac{1}{2v_{2}}\end{array}\right]$$

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]-... ...cal{\left\vert\widehat{A}\right\vert}\left(W_{i+1}-W_{i}\right)$$

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