4.4 Équation de convection-diffusion [Link]
suivant: 5. Annexe: EDP et monter: 4. Schémas différences finies précédent: 4.3 Équation
des ondes Table des matières
Sous-sections
4.4.1 Problème physique: dispersion d'un polluant
4.4.2 Étude de la solution exacte
[Link] Analyse en ordre de grandeur
[Link] Diffusion convection d'une gaussienne
[Link] Modes propres de diffusion
4.4.3 Discrétisation par différences finies
4.4.4 Méthode des directions alternées implicites
[Link] Stabilité et précision du schéma ADI
4.4.5 Expérimentation numérique avec Matlab
4.4 Équation de convection-diffusion
4.4.1 Problème physique: dispersion d'un polluant
On considère le problème de la dispersion d'un polluant à la surface d'un liquide en mouvement
(figure 4.15). Le polluant est miscible dans le liquide, mais plus léger. On suppose que la
vitesse du fluide est horizontale. On peut alors considérer que le polluant diffuse
essentiellement à la surface, et négliger les variations suivant la vertical .
Figure 4.15: diffusion d'une
tache de polluant
En notant la fraction massique de polluant et la densité du fluide, l'équation d'équilibre
pour traduit que la variation temporelle de la quantité de polluant dans un
volume élémentaire est égale à un bilan de flux de sur les facettes du volume. A
travers une facette de surface et de normale sortante , il existe un flux de diffusion
( est le coefficient de diffusion ) et un flux de convection ( est la
vitesse du fluide). L'équation d'équilibre s'écrit:
1 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
Compte tenu de l'équation de conservation de la masse du fluide:
cette équation s'écrit:
Compte tenu de l'hypothèse d'indépendance des quantités par rapport à , et en supposant en
outre que la densité et le coefficient de diffusion sont constants, il vient:
A cette équation il faut ajouter la condition initiale et des conditions aux
limites sur les frontières du domaine . On suppose que le polluant se trouve initialement à
l'intérieur du domaine et a pour dimension caractéristique (figure 4.15). On distingue les
frontières entrantes , i.e. telle que et les frontières sortantes , i.e. telle que
. Si la convection est prépondérante sur la diffusion, le flux entrant sur les frontières
est un flux de convection qui apporte du fluide non pollué dans le domaine . La condition
sur est une condition de Dirichlet . Sur les frontières , le fluide transporte le
polluant vers l'extérieur: on impose alors une condition aux limites de Neumann , qui
autorise la sortie du polluant du domaine .
Pour un domaine carré de dimension , et une vitesse avec des composantes positives
et , la frontière correspond aux deux cotés et et la frontière
aux deux autres cotés opposés et .
En notant , le problème modèle s'écrit pour un domaine carré de dimension
caractéristique :
Trouver tel que:
sur (4.26)
2 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
sur et sur
4.4.2 Étude de la solution exacte
Le problème (4.26) est un problème d'évolution parabolique caractéristique des problèmes de
mécanique des fluides, avec un terme de convection et un terme de diffusion. Nous allons tout
d'abord effectuer une analyse d'ordre de grandeur de chacun de ces termes.
[Link] Analyse en ordre de grandeur
Soit la dimension caractéristique de la tache initiale, l'ordre de grandeur des différents
termes de l'équation (4.26) s'écrit, en notant un temps caractéristique:
Si on considère uniquement la diffusion, le temps caractéristique vaut:
On retrouve le temps caractéristique de diffusion (relation c1tdiff) calculé pour l'équation de la
chaleur. Le temps caractéristique vérifie une relation du type (à rapprocher de la condition
de stabilité c1eq36 pour une équation de diffusion):
(4.27)
C'est le temps caractéristique de décroissance exponentielle des phénomènes de diffusion.
Si on considère uniquement la convection, le temps caractéristique vaut:
Ce temps correspond au temps de transport de la tache par le champ de vitesse sur une
distance égale à la dimension de la tache. Ce temps caractéristique vérifie donc une relation
du type (à rapprocher de la condition de stabilité de Courant c3eq12):
3 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
(4.28)
Pour notre problème, on peut définir un autre temps caractéristique de convection: le temps de
sortie du polluant hors du domaine . Il est basé sur la dimension caractéristique du
domaine et vérifie:
(4.29)
Enfin le rapport entre les temps caractéristiques de diffusion et de convection est le
nombre de Péclet:
(4.30)
qui caractérise l'importance relative du terme de convection par rapport au terme de diffusion.
[Link] Diffusion convection d'une gaussienne
Considérons la condition initiale suivante:
(4.31)
qui décrit une tache gaussienne d'amplitude 1 centrée en , de rayon .
S'il n'y a pas de convection, cette tache diffuse de façon auto-similaire , i.e. son amplitude
diminue et son rayon augmente en conservant une forme gaussienne:
En utilisant la conservation globale de dans tout le domaine:
on en déduit la relation entre l'amplitude et le rayon :
4 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
(l'intégrale d'une gaussienne vaut: ).
On cherche donc une solution de diffusion de l'équation (4.26) (avec ) sous la forme:
En reportant cette relation dans l'équation (4.26) on obtiens l'équation d'évolution de :
dont la solution vérifiant est:
(4.32)
La solution de diffusion de l'équation (4.26) s'écrit donc:
(4.33)
L'amplitude de cette gaussienne décroît donc suivant la loi:
(4.34)
En prenant en compte la convection par un champ de vitesse sans cisaillement, cette tache
gaussienne est transportée sans déformation et diffuse le long des trajectoires du champ de
vitesse comme précédemment. Pour un champ de vitesse constant, les trajectoires sont des
droites:
la solution de convection-diffusion de l'équation (4.26) s'écrit donc:
5 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
(4.35)
Cette solution est une solution en milieu infini et ne tiens pas compte des conditions aux limites
du problème (4.26). Elle constitue cependant une bonne approximation de la solution, si la
dimension de la tache est petite devant la dimension du domaine .
[solution a t=0 et t=0.5] [amplitude et taille ]
Figure: convection diffusion d'une gaussienne (4.35)
On a tracé l'évolution de cette solution (4.35) sur la figure (4.16) pour , et
. Pour ces valeurs des paramètres, un domaine de longueur , une position
initiale et , le temps caractéristique de diffusion vaut , celui de
convection , et le nombre de Péclet . Le problème est donc à convection
dominante. Au bout d'un temps , la tache est à la frontière du domaine
avec une amplitude qui a diminuée d'un tiers: .
[Link] Modes propres de diffusion
Pour rechercher des solutions vérifiant les conditions aux limites, on détermine tout d'abord les
modes propres de diffusion en utilisant la méthode de séparation de variable décrite au
6 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
paragraphe c1analytique. Le calcul est identique, et on montre facilement que les modes
propres sont les fonctions suivantes:
(4.36)
qui vérifient les conditions aux limites:
La solution générale de diffusion est alors une combinaison linéaire de ses modes:
[t=0] [Amplitude ]
Figure 4.17: mode propre de diffusion
On a tracé sur la figure (4.17) le mode et , ainsi que l'évolution temporelle de son
amplitude pour les mêmes paramètres que précédemment ( , ). Sur un temps
, l'amplitude de ce mode décroît de .
Si on prend en compte la convection dans le cas d'un champ de vitesse sans cisaillement, la
solution initiale est convectée sans déformation et diffusée le long des trajectoires. Par contre,
il n'existe pas de solutions analytiques simples qui vérifient les conditions aux limites de (4.26).
En considérant une taille de structure , le temps caractéristique de diffusion vaut
, le temps caractéristique de convection , et le nombre de Péclet .
4.4.3 Discrétisation par différences finies
7 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
Nous avons vu dans les chapitres précédents qu'une discrétisation précise de problème
parabolique est le schéma de Cranck Nicholson. Appliquée à l'équation (4.26), il s'écrit pour un
maillage cartésien de points suivant et points suivant , et de pas et :
(4.37)
C'est un schéma inconditionnellement stable d'ordre 2 en temps et en espace, i.e. en
.
A chaque itération en temps, on a à résoudre un système d'équations linéaire , de
inconnues . La matrice est une matrice penta-diagonale, qui a la même
structure que la matrice du laplacien au paragraphe c4matlap. Pour des très gros maillages, le
coût de résolution de ce système linéaire en utilisant les méthodes de résolution du paragraphe
c4matlap peut devenir rapidement prohibitif.
On va donc étudier dans le paragraphe suivant une méthode alternative: la méthode des
directions alternées implicites.
4.4.4 Méthode des directions alternées implicites
Le principe des méthodes des directions alternées implicites, notées ADI (ADI=Alternated
Directions Implicited est le sigle classique des directions alternées en anglais), est de
décomposer les opérateurs spatiaux suivant les directions d'espace et . On écrit l'équation
(4.26) sous la forme symbolique suivante:
(4.38)
où et sont les deux opérateurs suivants:
8 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
(4.39)
En notant , , les solutions au temps , et
, les développements limitées de peuvent s'écrire de façon symbolique:
On a utilisé dans ces relations le fait que et sont solutions de l'équation exacte (4.26)
pour remplacer en fonction de et .
En combinant ces deux équations, il vient:
(4.40)
En utilisant un développement au premier ordre des exponentielles,
on obtiens le schéma suivant, dans lequel il suffit d'inclure l'approximation spatiale des
opérateurs et :
(4.41)
C'est le schéam classique de Cranck Nicholson (4.37).
Pour les schémas ADI, on effectue tout d'abord une factorisation formelle dans (4.40):
avant le développement limité des exponentielles:
9 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
(4.42)
Formellement, on a un schéma de type Cranck Nicholson d'ordre 2 mais avec une erreur de
troncature différente. Pour résoudre, on introduit la solution intermédiaire telle que:
(4.43)
Ces deux équations sont équivalentes à l'équation (4.42). Pour s'en convaincre, il suffit de
multiplier la première par et la seconde par et de les combiner. L'intérêt de
cette procédure par rapport à Cranck Nicholson classique est que dans la première équation
(4.43), on est implicite suivant (i.e. suivant ) et explicite suivant (i.e. suivant ), et
vice-versa dans la seconde. La résolution de ces deux équations sera donc plus facile que la
résolution du schéma de Cranck Nicholson (4.41), dans lequel on est implicite suivant les 2
directions et .
On discrétise ensuite les opérateurs et avec des différences finies centrées sur un
maillage cartésien de points suivant et points suivant , et de pas et :
On obtiens le schéma ADI suivant pour l'équation (4.26):
(4.44)
10 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
(4.45)
Figure 4.18: schéma ADI4.44 et 4.45
Ces équations ADI correspondent à deux discrétisations de l'équation (4.26) à avec un pas
en temps . La première correspond à une discrétisation entre et avec une
implicitation dans la direction et la seconde à une discrétisation entre et avec une
implicitation dans la direction (figure 4.18). La variable intermédiaire corresponds donc
à une approximation de .
La première équation (4.44) couple les valeurs inconnues par ligne (figure 4.18), i.e. les
valeurs inconnues sur la ligne sont indépendantes des valeurs sur les autres
lignes. Elles sont solutions du système linéaire tri-diagonal de dimension
suivant:
(4.46)
11 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
avec , , , et , ,
Le second membre peut s'écrire sous la forme d'un produit matrice vecteur
avec une matrice tridiagonale:
Pour déterminer les valeurs inconnues , il faut donc résoudre systèmes linéaires tri-
diagonaux de dimension .
De même la seconde équation (4.45) couple les valeurs inconnues par colonne (figure
4.18), i.e. les valeurs inconnues sur la colonne sont indépendantes des valeurs
sur les autres colonnes. Elles sont solutions du système linéaire tri-diagonal de
dimension suivant:
(4.47)
avec , , , et , ,
Le second membre peut encore s'écrire sous la forme d'un produit matrice vecteur
:
12 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
Pour déterminer valeurs inconnues , il faut donc résoudre systèmes linéaires tri-
diagonaux de dimension .
A chaque itération en temps, on résoud systèmes tri-diagonaux de rang et systèmes
tri-diagonaux de rang , ce qui est beaucoup plus efficace que la résolution d'un seul système
linéaire de rang .
Pour les conditions aux limites, il faut modifier la première et la dernière ligne de ces systèmes
linéaires. Pour les conditions de Dirichlet sur , on modifie la première ligne de et :
et (4.48)
Cette condition fixe la valeur de la première ligne de : et la première
colonne de : .
Pour les conditions de Neumann sur , on utilise une condition miroir, qui modifie la dernière
ligne de et :
et (4.49)
De même le second membre de l'étape 1 pour la ligne est modifié:
(4.50)
ainsi que le second membre de l'étape 2 pour la colonne :
(4.51)
[Link] Stabilité et précision du schéma ADI
13 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
[Link].1 Étude de la stabilité:
L'étude de la stabilité utilise le programme Maple 4.4.4.
programme Maple 4.4.4: Etude de la stabilité du schèma ADI (4.44)
> restart:
# Equation de convection-diffusion
> diff(U(x,y,t),t)+V1*diff(U(x,y,t),x)+V2*diff(U(x,y,t),y)=
> kappa*(diff(U(x,y,t),x$2)+diff(U(x,y,t),y$2));eq:=%:
# Schema ADI
> (U[ns,i,j]-U[n,i,j])/(dt/2)+V1*(U[ns,i+1,j]-U[ns,i-1,j])/
(2*dx)+V2*(U[n,i,j+1]-U[n,i,j-1])/(2*dy)=kappa*((U[ns,i+1,j]
-2*U[ns,i,j]+U[ns,i-1,j])/dx^2+(U[n,i,j+1]-2*U[n,i,j]+
U[n,i,j-1])/dy^2); eqh1:=%:
> (U[n+1,i,j]-U[ns,i,j])/(dt/2)+V1*(U[ns,i+1,j]-U[ns,i-1,j])/
(2*dx)+V2*(U[n+1,i,j+1]-U[n+1,i,j-1])/(2*dy)=kappa*((U[ns,i+1,j]
-2*U[ns,i,j]+U[ns,i-1,j])/dx^2+(U[n+1,i,j+1]-2*U[n+1,i,j]+
U[n+1,i,j-1])/dy^2); eqh2:=%:
# Etude de la Stabilite
> Up:=(n,i,j)->Psi[n]*exp(I*omega[1]*i*dx)*exp(I*omega[2]*j*dy);
> subs(U[ns,i,j]=Up(ns,i,j),U[n,i,j]=Up(n,i,j),
U[ns,i+1,j]=Up(ns,i+1,j),U[ns,i-1,j]=Up(ns,i-1,j),
U[n,i,j+1]=Up(n,i,j+1),U[n,i,j-1]=Up(n,i,j-1),eqh1):
> rel1:=simplify(%*exp(-I*omega[1]*i*dx)*exp(-I*omega[2]*j*dy)):
> Psi[ns]/Psi[n]=solve(subs(Psi[ns]=G*Psi[n],rel1),G);rel11:=%:
#
> subs(U[ns,i,j]=Up(ns,i,j),U[n+1,i,j]=Up(n+1,i,j),
U[ns,i+1,j]=Up(ns,i+1,j),U[ns,i-1,j]=Up(ns,i-1,j),
U[n+1,i,j+1]=Up(n+1,i,j+1),U[n+1,i,j-1]=Up(n+1,i,j-1),eqh2):
> rel2:=simplify(%*exp(-I*omega[1]*i*dx)*exp(-I*omega[2]*j*dy)):
> Psi[n+1]/Psi[ns]=solve(subs(Psi[n+1]=G*Psi[ns],rel2),G);
rel22:=%:
# Facteur d'amplification
> G:=rhs(rel11)*rhs(rel22);
# Etude de chaque terme
> A1=r1*(1-cos(omega[1]*dx));A2=r2*(1-cos(omega[2]*dy));
> rel1:={%,%%}: rel11:={r1=kappa*dt/dx^2,r2=kappa*dt/dy^2}:
> B1=CFL1*sin(omega[1]*dx)/2;B2=CFL2*sin(omega[2]*dy)/2;
> rel2:={%,%%}: rel22:={CFL1=V1*dt/dx,CFL2=V2*dt/dy}:
> 'G'=((A1-1)+I*B1)/((A1+1)+I*B1)*((A2-1)+I*B2)/((A2+1)+I*B2);
GG:=rhs(%):
# Verification
> G=subs(rel1,rel2,rel11,rel22,GG):simplify(%):rhs(%)-lhs(%);
# Calcul du carre du module de G
> GG;
> G1:=((A1-1)+I*B1)/((A1+1)+I*B1);
G1M=((A1-1)^2+B1^2)/((A1+1)^2+B1^2);
> G2:=((A2-1)+I*B2)/((A2+1)+I*B2);
G2M=((A2-1)^2+B2^2)/((A2+1)^2+B2^2);
# G1M et G2M sont donc <1 ==> donc stabilite
On définit les deux étapes (4.44 et 4.45) du schéma ADI (lignes 6 à 13), puis on introduit une
perturbation décomposée en mode de Fourier (ligne 15), que l'on introduit dans les deux
équations (lignes 16 et 22). On calcul l'amplification du mode pour chacune de ces équations:
i.e. pour la première (ligne 20) et pour la seconde (ligne 22). D'où l'on déduit le
facteur d'amplification global du schéma (ligne 29):
14 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
avec , , , .
Le carré du module de vaut:
Il est plus petit que 1 puisque et sont positifs.
Le schéma ADI (4.44 et 4.45) est donc inconditionnellement stable.
[Link].2 Étude de la consistance:
Pour étudier la consistance, on utilise le programme Maple 4.4.4, qui est la suite du
programme précédent 4.4.4.
programme Maple 4.4.4: Etude de la consitance du schèma ADI (4.44)
# Erreur de troncature dvt de taylor autour de t=n+1/2
# Equation equivalente
> 1/2*eqh1+1/2*eqh2:eqh12:=lhs(%)-rhs(%):
# subsitution de la solution exacte dans cette equation
> Uex:=(r,p,q)->U(x+(p-i)*dx,y+(q-j)*dy,t+(r-n-1/2)*dt);
> subs(U[n,i,j]=Uex(n,i,j),
U[n,i+1,j]=Uex(n,i+1,j),U[n,i-1,j]=Uex(n,i-1,j),
U[n,i,j+1]=Uex(n,i,j+1),U[n,i,j-1]=Uex(n,i,j-1),
U[n+1,i,j]=Uex(n+1,i,j),
U[n+1,i+1,j]=Uex(n+1,i+1,j),U[n+1,i-1,j]=Uex(n+1,i-1,j),
U[n+1,i,j+1]=Uex(n+1,i,j+1),U[n+1,i,j-1]=Uex(n+1,i,j-1),
U[ns,i,j]=Uex(n+1/2,i,j),
U[ns,i+1,j]=Uex(n+1/2,i+1,j),U[ns,i-1,j]=Uex(n+1/2,i-1,j),
U[ns,i,j+1]=Uex(n+1/2,i,j+1),U[ns,i,j-1]=Uex(n+1/2,i,j-1),
eqh12); eqh3:=%:
# Developpement de Taylor autour de t=n+1/2
> k:=6:
> U(x,y,t-dt/2)=convert(mtaylor(U(x,y,t-dt/2),[dt],k),diff):
S1:=%:
> U(x,y+dy,t-dt/2)=convert(mtaylor(U(x,y+dy,t-dt/2),[dy,dt],
k),diff): S2:=%:
> U(x,y-dy,t-dt/2)=convert(mtaylor(U(x,y-dy,t-dt/2),[dy,dt],
k),diff):S3:=%:
> U(x+dx,y,t-dt/2)=convert(mtaylor(U(x+dx,y,t-dt/2),[dx,dt],
k),diff):S4:=%:
> U(x-dx,y,t-dt/2)=convert(mtaylor(U(x-dx,y,t-dt/2),[dx,dt],
k),diff):S5:=%:
> U(x,y+dy,t)=convert(mtaylor(U(x,y+dy,t),[dy],k),diff):
S6:=%:
> U(x,y-dy,t)=convert(mtaylor(U(x,y-dy,t),[dy],k),diff):
S7:=%:
> U(x+dx,y,t)=convert(mtaylor(U(x+dx,y,t),[dx],k),diff):
S8:=%:
> U(x-dx,y,t)=convert(mtaylor(U(x-dx,y,t),[dx],k),diff):
S9:=%:
> U(x,y,t+dt/2)=convert(mtaylor(U(x,y,t+dt/2),[dt],k),diff):
S10:=%:
> U(x,y+dy,t+dt/2)=convert(mtaylor(U(x,y+dy,t+dt/2),[dy,dt],
k),diff):S11:=%:
> U(x,y-dy,t+dt/2)=convert(mtaylor(U(x,y-dy,t+dt/2),[dy,dt],
k),diff):S12:=%:
> U(x+dx,y,t+dt/2)=convert(mtaylor(U(x+dx,y,t+dt/2),[dx,dt],
k),diff):S13:=%:
> U(x-dx,y,t+dt/2)=convert(mtaylor(U(x-dx,y,t+dt/2),[dx,dt],
k),diff): S14:=%:
# substitution dans l'équation discrete - equation exacte
15 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
> subs(S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14,eqh3)-
(lhs(eq)-rhs(eq)):
> simplify(%):collect(%,{dt,dx,dy});
# Schema d'ordre 2 en dt,dx et dy
Pour cela on fait la demi somme des deux équations (4.44 et 4.45) pour obtenir une équation
discrète équivalente à l'équation exacte (4.26) (ligne 3):
(4.52)
En comparant cette équation (4.52) au schéma de Cranck Nicholson (4.37), on peut retrouver
ce dernier en remplaçant simplement dans (4.52) les valeurs de par la moyenne
. Ce schéma ADI est donc bien équivalent au schèma de Cranck Nicolson. Les
valeurs intermédiaires sont des approximations de la solution à .
Dans l'équation équivalente (4.52), on substitue la solution approchée par la solution exacte
(lignes 6 à 15), et on effectue des développements limités autour de
(lignes 18 à 45). Après soustraction de l'équation exacte, on obtient l'erreur de troncature ,
qui est en .
Le schéma ADI (4.44 et 4.45) est donc consistant avec l'équation exacte (4.26), et est d'ordre 2
en temps et en espace, i.e. en .
4.4.5 Expérimentation numérique avec Matlab
Le programme Matlab (4.4.5) implémente le schéma ADI en utilisant une programmation
matricielle.
programme Matlab 4.4.5: Résolution numérique de l'équation (4.27)
16 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
% resolution ADI
clear
% maillage
L=1; H=1;
Nx=51; dx=L/(Nx-1); X=[0:dx:L];
Ny=51; dy=H/(Ny-1); Y=[0:dy:H];
% champ vitesse
v1=1; v2=1;
% parametres
kappa=0.01; dt=0.001; Tf=0.5; nit=round(Tf/dt)
% champ initial
delta=(0.1)^2; x0=L/4; y0=H/4;
Ui=exp(-(X-x0).^2/delta)'*exp(-(Y-y0).^2/delta);
Un=Ui; Us=Ui; Un1=Ui;
% matrices 3D
I1=ones(1,Nx); I2=ones(1,Ny);
A1 =[-kappa/dx^2-v1/(2*dx); 2/dt+2*kappa/dx^2; ...
-kappa/dx^2+v1/(2*dx)]*I1;
C1=[ kappa/dy^2+v2/(2*dy); 2/dt-2*kappa/dy^2; ...
kappa/dy^2-v2/(2*dy)]*I1;
A2 =[-kappa/dy^2-v2/(2*dy); 2/dt+2*kappa/dy^2; ...
-kappa/dy^2+v2/(2*dy)]*I2;
C2=[ kappa/dx^2+v1/(2*dx); 2/dt-2*kappa/dx^2; ...
kappa/dx^2-v1/(2*dx)]*I2;
% C.L de Dirichlet en i=1 j=1
A1(:,1)=0; A1(2,1)=1; C1(:,1)=0;
A2(:,1)=0; A2(2,1)=1; C2(:,1)=0;
% C.L de Neumann sur les frontieres i=Nx j=Ny
A1(1,Nx)=A1(1,Nx)+A1(3,Nx); A1(3,Nx)=0;
A2(1,Ny)=A2(1,Ny)+A2(3,Ny); A2(3,Ny)=0;
% iterations en temps
for it=1:nit
% 1ere etape ADI
Us(1:Nx,1)=0;
for j=2:Ny-1
B1=C1(1,1:Nx)'.*Un(1:Nx,j-1)+C1(2,1:Nx)'.*Un(1:Nx,j)+...
C1(3,1:Nx)'.*Un(1:Nx,j+1);
Us(1:Nx,j)=tridiag(A1,B1);
end;
B1=C1(2,1:Nx)'.*Un(1:Nx,Ny)+(C1(1,1:Nx)+...
C1(3,1:Nx))'.*Un(1:Nx,Ny-1);
Us(1:Nx,Ny)=tridiag(A1,B1);
% 2nd etape ADI
Un1(1,1:Ny)=0;
for i=2:Nx-1
B2=C2(1,1:Ny)'.*Us(i-1,1:Ny)'+C2(2,1:Ny)'.*Us(i,1:Ny)'+...
C2(3,1:Ny)'.*Us(i+1,1:Ny)';
Un1(i,1:Ny)=tridiag(A2,B2)';
end;
B2=C2(2,1:Ny)'.*Us(Nx,1:Ny)'+...
(C2(1,1:Ny)+C2(3,1:Ny))'.*Us(Nx-1,1:Ny)';
Un1(Nx,1:Ny)=tridiag(A2,B2)';
% iteration suivante
Un=Un1;
end;
Les paramètres du calcul sont définis aux lignes 4 à 10. Les matrices tridiagonales sont
construites sur les lignes 17 à 24, puis on applique les conditions aux limites (lignes 26 à 30).
La boucle en temps (lignes 32 à 55) inclus les deux étapes ADI et utilise la fonction tridiag
(lin3D) pour la résolution des systèmes linéaires tri-diagonaux.
[Link].1 Mode propre de diffusion:
Pour valider ce programme, nous avons tout d'abord simuler la diffusion du mode propre (4.36)
et avec un maillage de points dans chaque direction et un paramètre
. La solution calculée au bout d'un temps avec est tracée sur la
figure (4.19). L'allure de la solution (figure 4.19a) coïncide bien avec la solution exacte (figure
4.17a) , ce que confirme le tracé de l'évolution temporelle de la solution au point ,
17 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
comparée à la solution exacte (figure 4.19b).
[ ] [Amplitude en x=1 et y=1]
Figure 4.19: solution de diffusion avec le schéma ADI (mode propre , )
Pour tester la précision d'intégration en temps du schéma, nous avons calculer l'erreur au point
( , ) au bout du temps (de l'ordre du temps caractéristique de diffusion), en
fonction du pas d'intégration en temps . Le résultat de la figure (4.20a) montre que pour les
pas en temps choisis l'erreur est quasiment indépendante du pas en temps , et est donc
essentiellement une erreur de discrétisation spatiale. On note que les pas en temps choisis sont
tels que le pas en temps est beaucoup plus faible que le temps caractéristique de diffusion:
[Erreur temporelle ] [erreur spatiale ]
Figure 4.20: Erreur numérique du schéma ADI (cas de diffusion)
Nous l'avons vérifié en faisant varier le nombre de points du maillage de à
avec un pas en temps fixé. L'évolution de l'erreur en fonction du pas de
discrétisation spatiale est tracée sur la figure (4.20b), et on constate que l'erreur
décroît en .
Pour cette condition initiale, nous avons aussi fait une simulation avec une vitesse de
18 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
convection non nulle . Ce cas correspond à un nombre de Péclet . Les
iso-valeurs de la solution sont tracées sur la figure (4.21).
[ ] [ ] [ ]
Figure 4.21: Iso-valeurs de la solution (cas de convection-diffusion)
On note la convection sans déformation de la solution initiale, ce qui confirme que la condition
aux limites sur autorise la sortie des structures hors du domaine. On a aussi comparé
l'évolution temporelle du maximum de la solution exacte de diffusion et du maximum de la
solution calculée, que l'on a tracé sur la figure (4.22a). On vérifie ainsi que la décroissance de
la solution est une décroissance visqueuse.
[Evolution temporelle] [Erreur fonction de ]
Figure 4.22: Erreur numérique du schéma ADI
Nous avons ensuite étudié l'influence du pas d'intégration en temps , en traçant sur la
figure (4.22b) l'écart en fonction de entre le maximum de la solution exacte de diffusion et
le maximum de la solution calculée. On constate que cet écart croît rapidement en dessus
d'une valeur . Cette valeur est justement de l'ordre de grandeur du temps
caractéristique de convection .
En conclusion sur cette simulation, on note que le choix des paramètres numériques a été fixé
19 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
par la physique du problème, et non par des conditions numériques de stabilité:
1. la maillage permet de décrire finement la condition initiale
2. le pas en temps vérifie et .
[Link].2 Convection d'une gaussienne:
Le second cas de calcul correspond à la condition initiale gaussienne (4.31), avec les mêmes
paramètres qu'au paragraphe c4gpar.
[solution numérique à et ] [Amplitude]
Figure 4.23: Convection d'une gaussienne
Pour un maillage et un pas en temps , on a tracé sur la figure (4.23a) la
solution à et à . Cette solution se compare très bien avec la solution exacte tracée
sur la figure (4.16). On a aussi comparé l'évolution temporelle de l'amplitude de la tache
gaussienne calculée avec le schéma ADI avec l'expression analytique (4.34). Ces deux courbes
coïncident, ainsi que le montre la figure (4.23b).
[isovaleurs à ] [isovaleurs à ]
[profil à ]
Figure 4.24: Convection d'une gaussienne
20 sur 21 24/01/2019 à 12:34
4.4 Équation de convection-diffusion [Link]
Nous avons effectué une seconde simulation avec un coefficient de diffusion plus petit
. Dans ce cas, la solution numérique présente des oscillations (figure 4.24). Pour
cette valeur de , le nombre de Péclet de maille:
vaut au lieu de avec la valeur de précédente. Sur le tracé des iso-valeurs
(4.24a), on constate l'apparition de légères oscillations, caractérisées par la présence de
nombreuses lignes iso-valeurs , qui n'existent pas à (figure 4.24b). Le tracé d'un
profil à (figure 4.24c) montre bien l'apparition d'une oscillation numérique au pied de
la tache gaussienne. Ces oscillations numériques sont de même nature que celles étudiées au
paragraphe c3centre du chapitre précédent. Elles apparaissent dès que le Péclet de maille
devient plus grand que 2 et indiquent que le maillage n'est plus suffisamment fin pour capter la
solution de convection avec ce schéma centré.
suivant: 5. Annexe: EDP et monter: 4. Schémas différences finies précédent: 4.3 Équation
des ondes Table des matières
Pr. Marc BUFFAT
[Link]ffat@[Link]
2008-04-07
21 sur 21 24/01/2019 à 12:34