3.3 Équation de transport

3.3.1 Problème physique: convection dans un fluide

On considère le transport par un fluide d'une quantité scalaire définie par unité de volume $u(x,t)$ . On suppose que le champ de vitesse $V$ est unidimensionnel, que le scalaire $u(x,t)$ ne diffuse pas et est uniquement transporté par le fluide. La quantité $u(x,t)$ se conserve donc le long des trajectoires,i.e.:

$$\frac{Du}{Dt}=0$$ (3.12)

où $\frac{D}{Dt}$ représente la dérivée lagrangienne, i.e. la dérivée suivant la trajectoire des particules fluides qui s'écrit en fonction des dérivées partielles en temps et en espace:

$$\frac{D()}{Dt}=\frac{\partial()}{\partial t}+V\frac{\partial()}{\partial x}$$

L'équation (3.12) est donc l'équation de convection pure suivante:

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

qui traduit la conservation de $u(x,t)$ le long des trajectoires. Si on connait les trajectoires, il suffit de connaitre la solution en un point d'une trajectoire pour connaitre la solution sur toute la trajectoire.

3.3.2 Étude de la solution analytique

Si le champ de vitesse $V$ est constant, les trajectoires dans l'espace $(x,t)$ sont des droites de pente $V^{-1}$ :

$$x-V t=cste$$

On a tracé ses trajectoires (pour le cas $V>0$ ) sur la figure (3.7). Pour déterminer la solution en un point $(x,t)$ du domaine $[0,L]*[0,\tau]$ , on détermine tout d'abord la trajectoire passant par ce point. Si cette trajectoire a une intersection avec l'axe des $x$ (cas du point 1 sur la figure 3.7), la solution en $(x,t)$ est égale à la valeur de $u$ en ce point 1, i.e. est déterminée par la condition initiale $u_{0}(x)$ à $t=0$ . Si la trajectoire a une intersection avec l'axe des $t$ (cas du point 2 sur la figure 3.7), la solution en $(x,t)$ est égale à la valeur de $u$ en ce point 2, i.e. est déterminée par la condition aux limites $u_{1}(t)$ en $x=0$ .

$$u(x,t)=\left\{ \begin{array}{c} u_{0}(x-V t) \mbox{si }x-V t>0 u_{1}(t-\frac{x}{V}) \mbox{ si }x-V t<0\end{array}\right.$$ (3.14)

Figure 3.7: trajectoires dans l'espace $(x,t)$

La solution générale de l'équation de convection pure (3.13) est donc déterminée par des conditions aux intersections des frontières du domaine avec le pied des trajectoires. Si $V>0$ , cela correspond à la condition initiale $u_{0}$ et à la condition aux limites $u_{1}$ en $x=0$ . Si $V<0$ , cela correspond à la condition initiale $u_{0}$ et à la condition aux limites $u_{1}$ en $x=L$ .

Nous étudierons plus particulièrement deux cas:

• cas 1: convection d'une tache gaussienne:
  
  $$u_{0}(x)=e^{-\left(\frac{x-x_{0}}{\sigma}\right)^{2}}$$

• cas 2: convection d'un sinus:
  
  $$u_{0}(x)=\sin\frac{m\pi x}{L}$$

Dans le premier cas, le temps de simulation $\tau$ est tel que la tache reste dans le domaine, i.e. $\tau<\frac{L}{V}$ . La valeur de la solution sur les frontières $x=0$ et $x=L$ est alors égale à 0. Dans le second cas, on impose une condition de périodicité de la solution.

3.3.3 Discrétisation par différences finies

Nous allons étudier différents schémas pour la résolution numérique de l'équation (3.13). Les schémas étudiés s'écrivent sous la forme générique suivante:

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

qui conduit à la résolution à chaque itération en temps d'un système tri-diagonal:

$$\mathcal{A}[U^{n+1}]=\mathcal{C}[U^{n}]$$ (3.16)

Pour chaque schéma nous étudierons la stabilité et la consistance en utilisant le programme Maple suivant (3.3.3). Pour l'utiliser on rentre à la ligne 6 les coefficients $a$ et $b$ du schéma.

programme Maple 3.3.3: Analyse de la stabilité et de la consistance d'un schéma D.F. pour l'équation de convection pure 3.13

> restart;
# Etude de l'équation de convection pure
> diff(U(t,x),t)+V*diff(U(t,x),x)=0;eq:=%:
# Tests de Schema D.F.
# Définition des coefficients du schéma
> a:=[0,1/dt,0]; c:=[V/dx/2,1/dt,-V/dx/2];
> a[1]*U[n+1,i-1]+a[2]*U[n+1,i]+a[3]*U[n+1,i+1]=
     c[1]*U[n,i-1]+c[2]*U[n,i]+c[3]*U[n,i+1];
> eqh:=lhs(%)-rhs(%)=0:
# Stabilite
> Up:=(n,i)->Psi[n]*exp(I*omega*i*dx);
> subs(U[n+1,i]=Up(n+1,i),U[n+1,i+1]=Up(n+1,i+1),
   U[n+1,i-1]=Up(n+1,i-1),U[n,i]=Up(n,i),U[n,i-1]=Up(n,i-1),
   U[n,i+1]=Up(n,i+1),U[n,i-2]=Up(n,i-2),eqh);
> expand(simplify((%)*dt*exp(-I*omega*i*dx)));rel1:=%:
> simplify(subs(Psi[n+1]=G*Psi[n],rel1/Psi[n])):
> G=solve(%,G);rel2:=rhs(%):
# Coefficient d'amplification et le carré du module
> G:=simplify(subs(omega*dx=y,V=CFL*dx/dt,evalc(rel2)),trig);
> G2:=(evalc(Re(G))^2+evalc(Im(G))^2);
# Consistance
> Uex:=(p,q)->U(t+(p-n)*dt,x+(q-i)*dx);
> subs(U[n,i]=Uex(n,i),U[n,i+1]=Uex(n,i+1),U[n,i-1]=Uex(n,i-1),
      U[n+1,i]=Uex(n+1,i),U[n+1,i+1]=Uex(n+1,i+1),
      U[n+1,i-1]=Uex(n+1,i-1),lhs(eqh)); rel3:=%:
> expand(simplify(rel3-lhs(eq)));rel4:=%:
# Developpement en serie de Taylor
> U(t+dt,x)=convert(mtaylor(U(t+dt,x),[dt],4),diff);S1:=%:
> U(t+dt,x+dx)=convert(mtaylor(U(t+dt,x+dx),[dt,dx],4),diff):S2:=%:
> U(t+dt,x-dx)=convert(mtaylor(U(t+dt,x-dx),[dt,dx],4),diff):S3:=%:
> U(t,x+dx)=convert(mtaylor(U(t,x+dx),[dx],4),diff):S4:=%:
> U(t,x-dx)=convert(mtaylor(U(t,x-dx),[dx],4),diff):S5:=%:
# Erreur de troncature
> simplify(subs(S1,S2,S3,S4,S5,rel4)):
> ErrT:=collect(%,dt);

L'analyse de stabilité utilise la même approche que dans le programme (3.2.4) du paragraphe 3.2.4 . De même l'analyse de la consistance suit la même démarche que dans le programme (3.2.5) du paragraphe 3.2.5.

Pour les expérimentations numériques, nous utiliserons le programme Matlab (3.3.3) qui résout numériquement par différences finies l'équation de convection pure (3.13) sur le domaine $[0,1]$ avec un maillage de $N=201$ points, $V=1$ et une condition initiale gaussienne (cas 1 avec $x_{0}=0.2$ , $\delta=0.05$ ):

$$u_{0}(x)=e^{-\left(\frac{x-x_{0}}{\delta}\right)^{2}}$$

On calcule la solution au temps $t_{1}=0.5$ que l'on compare avec la solution exacte:

$$u(x,t_{1})=e^{-\left(\frac{x-V t_{1}-x_{0}}{\delta}\right)^{2}}$$

Pour utiliser le programme on rentre sur les lignes 13 et 14, les coefficients $a$ et $c$ du schéma (en utilisant les notations 3.15). Le temps $t_{1}$ est telle que la solution ne sort pas du domaine de calcul, et la solution en $x=0$ et $x=L$ reste non perturbée, i.e $u=0$ . Bien que la solution exacte ne dépende que de la valeur en $x=0$ (car $V>0$ ), on a besoin d'une condition aux limites numérique en $x=L$ , i.e pour $i=N$ . En effet le schèma différences finis 3.15 n'est pas applicable pour $i=N$ puisqu'il fait intervenir des valeurs hors maillage $u_{N+1}$ . Pour simplifier, nous supposerons que la solution en $x=L$ reste non perturbée et nous imposerons une condition aux limites $u=0$ aux deux extrémités $x=0$ et $x=L$ (lignes 20 à et 21).

programme Matlab 3.3.3: Simulation par D.F. de l'équation de convection pure 3.13

% resolution de l'equation de convection pure
clear; clf;
% parametres
V=1.0; L=1.0;
N=201; CFL=1/2
dx=L/(N-1);dt=CFL*dx/V;
X=[0:dx:L]; 
Tf=0.5;
% champ initial
x0=0.2; delta=0.05;
U0=exp(-((X-x0)/delta).^2);
% definition du schema D.F. operateur
a=[0 1 0]
c=[CFL/2 1 -CFL/2]
% construction des matrices tridiag
e=ones(1,N);
A=[a(1)*e; a(2)*e; a(3)*e];
C=[c(1)*e; c(2)*e; c(3)*e];
% C.L dirichlet en 1 et N
A(:,1)=0; A(2,1)=1;  C(:,1)=0; C(2,1)=1;
A(:,N)=0; A(2,N)=1;  C(:,N)=0; C(2,N)=1;
% iteration en temps
nit=round(Tf/dt)
Un1=U0;
for it=1:nit
  Un=Un1;
  B=C(1,1:N).*[0 Un(1:N-1)]+C(2,1:N).*Un(1:N)+C(3,1:N).*[Un(2:N) 0];
  Un1=tridiag(A,B);
end;
% trace solution finale
t1=nit*dt;
Ue=exp(-((X-V*t1-x0)/delta).^2);
clf; figure(1);
plot(X,U0,X,Un1,X,Ue);

3.3.4 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$$ (3.17)

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}$ comme indiqué sur le diagramme (3.8).

Figure 3.8: schéma explicite centrée pour l'équation de convection pure

Du schéma (3.17) on déduit les coefficients $a$ et $b$ pour le programme Maple (3.3.3):

$$a:=[0,\frac{1}{dt},0]; c:=[\frac{V}{2dx}, \frac{1}{dt}, -\frac{V}{2dx}];$$

Avec ce programme, on effectue l'analyse de consistance et stabilité. 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}$$ (3.18)

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

$$\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 (3.17) est donc inutilisable pour résoudre l'équation de convection pure (3.13).

Pour comprendre la raison de cette instabilité, nous allons faire l'étude de la consistance de ce schéma. Le programme Maple (3.3.3) donne l'erreur de troncature suivante:

$$ErrT=\frac{1}{2}\frac{\partial^{2}u}{\partial... ...c{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 (3.13), mais instable. Pour comprendre cette instabilité, étudions le premier terme de l'erreur:

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

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

$$\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}}$$ (3.20)

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

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

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é 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$$ (3.22)

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

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

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

3.3.5 Schéma implicite centré

La version implicite du schéma précédent s'écrit:

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

C'est un schéma implicite qui couple les valeurs inconnues $\{u_{i}^{n+1},u_{i-1}^{n+1},u_{i+1}^{n+1}\}$ à l'étape $t^{n+1}$ à la valeur connue $u_{i}^{n}$ à l'étape $t^{n}$ comme indiqué sur le diagramme (3.9).

Figure 3.9: schéma implicite centrée pour l'équation de convection pure

A partir de (3.24), on en déduit les coefficients $a$ et $b$ pour le programme Maple (3.3.3):

$$a:=[-\frac{V}{2dx}, \frac{1}{dt}, \frac{V}{2dx}]; c:=[0,\frac{1}{dt},0];$$

Avec ce programme, on calcul le carré du module du facteur d'amplification $G$ du schéma, qui s'écrit:

$$\left\vert G\right\vert^{2}=\frac{1}{1+C_{FL}^{2}\sin^{2}y} <1 \forall y$$

C'est exactement l'inverse du facteur d'amplification du schéma explicite. Il est toujours inférieur à 1, donc le schéma implicite est inconditionnellement stable.

Le programme Maple (3.3.3) donne l'erreur de troncature suivante:

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

Le schéma implicite est donc consistant à l'équation (3.13), d'ordre $O(dt,dx^{2})$ et inconditionnellement stable.

En utilisant la propriété (3.20) de la solution exacte, le premier terme de l'erreur de troncature s'écrit:

$$-\frac{V^{2}}{2}\frac{\partial^{2}u}{\partial x^{2}} dt$$

C'est l'opposé du premier de l'erreur de troncature du schéma explicite. Le schéma implicite est donc équivalent à l'équation de convection diffusion suivante:

$$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 $\kappa=v^{2}\frac{dt}{2}$ positif. Le schéma implicite est donc stable, mais diffusif. La solution numérique est convectée, mais décroît exponentiellement proportionnellement à un coefficient de “diffusion numérique”:

$$\kappa=V^{2}\frac{dt}{2}$$ (3.25)

On le vérifie numériquement en utilisant le programme Matlab (3.3.3) avec les paramètres suivants:

$$a=[-CFL/2 1 CFL/2]; c=[0 1 0];$$

Sur la figure (3.10), 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, en particulier pour les grandes valeurs du nombre de Courant $C_{FL}$ .

Figure 3.10: solution du schéma implicite (3.24) pour différents CFL

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

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

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

$$\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$$ (3.26)

C'est un schéma de Lax-Wendroff (c31laxwend). 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}$ comme indiqué sur le diagramme (3.11).

Figure 3.11: schéma explicite d'ordre 2 pour l'équation de convection pure

Les coefficients $a$ et $b$ pour le programme Maple (3.3.3) valent:

$$a:=[0,\frac{1}{dt},0] c=[\frac{V}{2dx... ...t}-2\frac{V^{2}dt}{2dx^{2}}\„-\frac{V}{2dx}+\frac{V^{2}dt}{2dx^{2}} ];$$

Avec ce programme, le carré du module du facteur d'amplification $G$ du schéma s'écrit 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$$ (3.27)

Le programme Maple (3.3.3) donne l'erreur de troncature, qui s'écrit, en tenant compte de la propriété (3.20) :

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

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

On le vérifie numériquement en utilisant le programme Matlab (3.3.3) avec les paramètres suivants:

$$a=[0 1 0]; c=[\frac{CFL}{2}+\frac... ...rac{V^{2}dt^{2}}{2dx^{2}}, -\frac{CFL}{2}+\frac{V^{2}dt^{2}}{2dx^{2}}];$$

Sur la figure (3.12a), 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}$ . Pour faire une analyse plus fine, on a tracé sur la figure (3.12b) 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 3.12: Solution numérique avec le schéma explicite d'ordre 2 (3.26)

[solution numérique]

[erreur à $t=0.5$ ]

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

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

Une autre façon de stabilisée le schéma explicite (3.24), 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}$$  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}$$  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$$ (3.29)

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}$ comme indiqué sur le diagramme (3.13). On note que la valeur avale $u_{i+1}^{n}$ n'intervient pas dans ce cas pour le calcul de $u_{i}^{n+1}$ .

Figure 3.13: schéma explicite décentrée pour l'équation de convection pure

Les coefficients $a$ et $b$ pour le programme Maple (3.3.3) s'écrivent:

$$a:=[0,\frac{1}{dt},0] c=[\frac{V}{dx}\„ \frac{1}{dt}-\frac{V}{dx}\„ 0];$$

Avec ce programme, on calcule 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$$ (3.30)

Le programme Maple (3.3.3) fournit l'erreur de troncature, qui s'écrit:

$$ErrT=\frac{1}{2}\frac{\partial^{2}u}{\partial... ...c{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é (3.20) :

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

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

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

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$ (3.31) si la condition de stabilité (3.30) est vérifiée. Si la condition de stabilité (3.30) 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 en utilisant le programme Matlab (3.3.3) avec les paramètres suivants:

$$a=[0 1 0]; c=[C_{FL}, 1-C_{FL}, 0];$$

Sur la figure (3.14), 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 3.14: Solution numérique avec le schéma explicite décentré (3.29)

3.3.8 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é 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 2 de la convection d'une onde $u_{0}(x)=e^{I\omega x}$ . La solution exacte (3.14) de l'équation de convection pure (3.13) s'écrit:

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

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}$$ (3.33)

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

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

Figure 3.15: Diffusion numérique du schéma d'ordre 1 (3.29) et d'ordre 2 (3.26)

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 avec le programme Matlab (3.3.8)

programme Matlab 3.3.8: Convection d'une onde par D.F.

% resolution de l'equation de convection pure
% condition de periodicite schema explicite
%clear; clf;
% parametres
V=1.0; L=pi;
N=101; CFL=0.9
dx=L/(N-1);dt=CFL*dx/V
X=[0:dx:L]; 
% champ initial
m=10;
Ui=sin(m*pi*X/L);
% definition du schema D.F. operateur
a=1;
cdiff=V^2/2*dt^2/dx^2
c=[CFL/2+cdiff 1-2*cdiff -CFL/2+cdiff]
% construction des matrices tridiag
e=ones(1,N);
C=[c(1)*e; c(2)*e; c(3)*e];
% iteration en temps
nit=200
Un1=Ui; Ue=Ui;
figure(1); clf; hold on;
p1=plot(X,Ue,'EraseMode','background');
p2=plot(X,Un1,'EraseMode','background');
axis([0 L -1 1]);
for it=1:nit
  Un=Un1;
  B=C(1,1:N).*[0 Un(1:N-1)]+C(2,1:N).*Un(1:N)+C(3,1:N).*[Un(2:N) 0];
  % periodicite
  B(1)=c(1)*Un(N-1)+c(2)*Un(1)+c(3)*Un(2);
  B(N)=B(1);
  Un1=B/a;
  % trace de la solution numerique et exacte
  t1=it*dt; Ue=sin(m*pi*(X-V*t1)/L);
  pause(0.01);
  set(p1,'YData',Ue);set(p2,'YData',Un1);
  drawnow;
end;
grid on;
hold off;

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}$$ (3.34)

dont les paramètres sont définis aux lignes 13 à 15 dans les tableaux $a$ et $c$ . Le programme permet de tracer l'évolution de la solution numérique $Un1$ et exacte $Ue$ au cours du temps en utilisant les techniques d'animation de Matlab (lignes 22 à 25 et lignes 33 à 37). Le schéma numérique étant explicite, la solution s'obtient simplement à chaque itération par multiplication matricielle (ligne 28):

$$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 (3.34) fait intervenir la valeur $u_{-1}^{n}$ , que l'on calcul avec la condition de périodicité (ligne 30): $u_{-1}^{n}=u_{N-1}^{n}$ . Pour le noeud $i=N$ , on impose la condition $u_{N}^{n}=u_{0}^{n}$ (ligne 31).

On a tracé sur la figure (3.16) la solution obtenue avec le schéma décentré d'ordre 1 (3.29) et le schéma d'ordre 2 (3.26). 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ée (mais qui a une dispersion plus faible).

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

[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 (3.26), on transforme l'erreur de troncature (3.28) 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}} dx^{2}+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\... ...ac{\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 (3.13)), 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 (3.26) n'est plus dispersif, ce que l'on vérifie bien numériquement.

3.3.9 Convection d'une discontinuité

La capture de 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 (3.13):

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

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

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

La solution exacte corresponds à la propagation de cette discontinuité avec la célérité $V$ :

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

Sur la figure (3.17), nous avons comparé la solution exacte avec la solution approchée calculée avec le schéma upwind (3.29) et le schéma Lax Wendroff (3.26) sur un maillage de $101$ points ( $dx=0.01$ ).

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

[schéma upwind]

[schéma 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), mais elle présente des oscillations. Si on raffine le maillage, les mêmes conclusions peuvent être tirées, comme le montre la figure (3.18) où on a tracé la solution numérique sur un maillage 10 fois plus fin $(dx=0.001)$ .

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

[Schéma Upwind]

[Schéma 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é 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.

En effet en décomposant la condition initiale en série de Fourier

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

la solution exacte correspond à la propagation de chacune des ondes $e^{I\omega_{k}x}$ avec la même célérité $V$ , ce qui permet de reconstruire la solution à l'instant $t$ de façon identique à la solution initiale:

$$u(x,t)=\sum_{k=0}^{\infty}C_{k}e^{I\omega_{k}(x-Vt)}=u_{0}(x-V t)$$

D'après l'étude précédente, les schémas numériques étudiés convectent les ondes en les amortissants (avec un coefficient de diffusion numérique $\kappa$ ) et en les déphasants (avec un déphasage $\varphi$ ):

$$e^{-\kappa\omega_{k}^{2}t}e^{I\omega_{k}(x-(V-\varphi\omega_{k}^{2})t)}$$

La solution numérique reconstruite à l'instant $t$ s'écrit:

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

Si le schéma est dissipatif ( $\kappa\neq0$ ), mais peu dispersif ( $\varphi\approx0$ ), l'amplitude des ondes hautes fréquences est amortie, mais elles sont convectées à la bonne vitesse. Les discontinuitées dans la solution initiale $u_{0}$ sont lissées au cours du temps, et la solution numérique présente des fronts qui s'étalent. C'est le cas du schéma décentré d'ordre 1 (3.29).

Par contre si le schéma est dispersif ( $\varphi\neq0$ ), et peu dissipatif ( $\kappa\approx0$ ), les ondes hautes fréquences sont peu amorties, mais elles sont déphasées. C'est le cas du schéma de Lax Wendroff (3.26). 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 ( $\varphi>0$ ), 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 (3.17) et (3.18) (la fréquence de l'oscillation augmente avec le nombre de points du maillage).