3.4 Équation de convection-diffusion

3.4.1 Problème physique

On veut déterminer la répartition de température $T(x)$ dans un fluide incompressible de vitesse $V$ et de conductivité thermique $\kappa$ entre 2 sections $x_{1}=0$ et $x_{1}=L$ où la température est maintenue constante. L'équation de conservation de l'énergie s'écrit pour un régime stationnaire:

$$V \frac{\partial T}{\partial x_{1}}=\kappa\frac{\partial^{2}T}{\partial x_{1}^{2}}$$

avec les conditions aux limites $T(0)=T_{0}$ et $T(L)=T_{1}$ . En notant $u=\frac{T-T_{0}}{T_{1}-T_{0}}$ , $x=\frac{x_{1}}{L}$ , le problème modèle s'écrit:

$$P_{e} \frac{\partial u}{\partial x}=\frac{\partial^{2}u}{\partial x^{2}} u(0)=0,u(1)=1$$ (3.35)

où $P_{e}=\frac{VL}{\kappa}$ est le nombre de Péclet, nombre sans dimension caractéristique du problème.

3.4.2 Étude de la solution analytique

La solution analytique de (3.35) est obtenue facilement avec Maple en utilisant la fonction $dsolve$ . Elle s'écrit:

$$u(x)=\frac{1-e^{P_{e}x}}{1-e^{P_{e}}}$$ (3.36)

Figure 3.19: solution exacte de l'équation de convection-diffusion (3.35)

L'allure de cette solution est tracée sur la figure (3.19) pour différentes valeurs du nombre de Péclet $P_{e}$ . Sur cette figure on constate que se développe une couche limite thermique dont l'épaisseur $\delta$ décroît lorsque le nombre de Péclet augmente. Plus précisément, en définissant l'épaisseur $\delta$ , comme étant la distance à la paroi du point où la température relative décroit de $90$ % , le rapport $\xi=\frac{\delta}{L}$ vérifie:

$$u(\xi)=\frac{1-e^{P_{e}(1-\xi)}}{1-e^{P_{e}}}=0.1$$

soit, en réarrangeant les termes:

$$1-\frac{e^{P_{e}}}{e^{P_{e}}-1}(1-e^{-P_{e}\xi})=0.1$$

Pour des grands nombres de Péclet, cette relation conduit à :

$$e^{-P_{e}\xi}=0.1\Rightarrow\xi=\frac{\log10}{P_{e}}$$

L'épaisseur de la couche limite thermique est donc inversement proportionnelle au nombre de $P_{e}$ . Le problème devient numériquement de plus en plus raide lorsque $P_{e}$ augmente.

3.4.3 Discrétisation avec un schéma centré d'ordre 2

La discrétisation par différences finies centrées de l'équation de convection diffusion (3.35) s'écrit:

$$P_{e} \frac{u_{i+1}-u_{i-1}}{2dx}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{dx^{2}}$$ (3.37)

Le problème étant stationnaire, on ne fait que l'analyse de la consistance du schéma (3.35). Ceci peut être fait avec Maple, et on obtiens l'erreur de troncature suivante:

$$ErrT=\left(\frac{1}{6}P_{e} \frac{\partial^{... ...rac{1}{12}\frac{\partial^{4}u}{\partial x^{4}}\right) dx^{2}+O(dx^{4})$$

que l'on peut transformer en dérivant 2 fois l'équation (3.35):

$$ErrT=\frac{1}{12}\frac{\partial^{4}u}{\partial x^{4}}dx^{2}+O(dx^{4})$$

Le schéma (3.37) est donc consistant à l'équation (3.35) et d'ordre 2.

Pour poursuivre l'analyse, nous allons déterminer la solution de l'équation aux différences (3.37) avec le programme Maple (3.4.3) pour un maillage de $N+1$ points.

programme Maple 3.4.3: Etude du schéma centré 3.37 pour l'équation 3.35

> restart:
# Etude de l'equation de convection diffusion
> Pe*diff(U(x),x)-diff(U(x),x$2)=0;eq:=%:
> U(0)=0,U(1)=1; CL:=%:
# determination de la solution exacte
> dsolve({eq,CL},U(x)):simplify(%):
> uex:=unapply(rhs(%),x);
# Schema de D.F. centree
> Pe*(U[i+1]-U[i-1])/(2*dx)-(U[i+1]-2*U[i]+U[i-1])/dx^2=0;eqh:=%:
# Solution exacte de l'equation D.F
> eqh;
> collect(simplify(dx^2*eqh),[U[i],U[i+1],U[i-1]]);eq1:=%:
> subs(Pe=Peh/dx,eq1);eq2:=%:
# Relation de recurrence a trois niveaux: 
# recherche des racines de l'equation caracteristique
> subs(U[i]=lambda,U[i-1]=1,U[i+1]=lambda^2,eq2);
> Lambda:=[solve(%,lambda)];
# Solution generale
> Ug:=unapply(C1*Lambda[1]^i+C2*Lambda[2]^i,i);
# Determination des constantes avec les C.L.
> Ug(0);Ug(N)=1;solve({%,%%},{C1,C2});rel:=%:
> Uh:=unapply(simplify(subs(rel,Ug(i))),i);
# Etude de la convergence vers la solution exacte
> subs(Pe=Peh*N,uex(i/N));Uexi:=%:
> Uh(i);
> subs(Peh=Pe*dx,i=k*N,N=1/dx,Uexi);subs(Peh=Pe*dx,i=k*N,N=1/dx,Uh(i));
> simplify(taylor(%%-%,dx,4));
# Donc convergence d'ordre 2

On commence par déterminer la solution exacte de l'équation (3.35) (lignes 6 à 7), puis on écrit le schéma aux différences finies (3.37) sous la forme d'une relation de récurrence à 3 niveaux (lignes 11 à 13):

$$\frac{P_{e}^{h}-2}{2} u_{i+1}+2u_{i}-\frac{P_{e}^{h}+2}{2} u_{i-1}=0$$

en introduisant le nombre de Péclet de maille $P_{e}^{h}$ :

$$P_{e}^{h}=P_{e}dx=\frac{V (Ldx)}{\kappa}$$

qui est un nombre de Péclet basé sur la longueur de la maille de calcul (ici $Ldx$ dans le domaine physique $[0,L]$ ).

La solution générale de cette relation de récurrence s'écrit en fonction des racines $\Lambda$ de l'équation caractéristique associée (calculées aux lignes 16 et 17) sous la forme (ligne 19):

$$u_{i}=C_{1}(\Lambda_{1})^{i}+C_{2}(\Lambda_{2})^{i}$$

Les constantes $C_{1}$ et $C_{2}$ sont déterminées par les conditions aux limites $u_{0}=0$ et $u_{N}=1$ (lignes 21 et 22). D'où la solution exacte de l'équation aux différences (3.37):

$$u_{i}=\frac{\left(\frac{2+P_{e}^{h}}{2-P_{e}^{h}}\right)^{i}-1}{\left(\frac{2+P_{e}^{h}}{2-P_{e}^{h}}\right)^{N}-1}$$ (3.38)

qu'il faut comparer à la solution exacte (3.36) calculée aux noeuds du maillage $x=i dx=\frac{i}{N}$ , qui s'écrit en introduisant $P_{e}^{h}=P_{e}N$ (lignes 24 et 25) :

$$u(\frac{i}{N})=\frac{e^{P_{e}^{h}i}-1}{e^{P_{e}^{h}N}-1}$$ (3.39)

Pour étudier la convergence, on se place en un point $x=idx=\frac{i}{N}$ fixé et on étudie l'écart $u(\frac{i}{N})-u_{i}$ lorsque $dx=1/N$ tends vers zéro. Pour cela on écrit tous les paramètres ( $P_{e}^{h}$ , $N$ , $i$ ) en fonction de $dx$ en notant que le rapport $\frac{i}{N}=k$ reste constant (car le point $x$ est fixé). On calcul ensuite un développement en série de Taylor de l'écart en fonction de $dx$ (lignes 26 à 27). Le résultat donné par Maple est le suivant:

$$u(\frac{i}{N})-u_{i}=O(dx^{2})$$

La solution numérique converge vers la solution exacte avec une précision d'ordre 2 (ce que nous donnait l'erreur de troncature).

L'expression de la solution numérique (3.37) montre que l'on approxime des exponentielles par des fonctions puissances. De plus, pour un maillage fixé (i.e. avec un nombre de points $N+1$ donné), on constate que le terme dans les fonctions puissances $\frac{2+P_{e}^{h}}{2-P_{e}^{h}}$ devient négatif lorsque l'on augmente le nombre de Péclet $P_{e}$ . D'un point à l'autre la solution change alors de signe et oscille. Plus précisément, pour des valeurs de $P_{e}^{h}$ telles que:

$$P_{e}^{h}=P_{e}\frac{1}{N}>2$$ (3.40)

la solution numérique présente des oscillations numériques (et non une instabilité numérique). Pour comprendre l'apparition de ces oscillations, considérons le maillage avec 3 points de la figure ci-dessous:

Les valeurs en $u_{0}$ et en $u_{2}$ sont fixées par les conditions aux limites et le seul degré de liberté est la valeur $u_{1}$ . Cette valeur $u_{1}$ doit être tel qu'au niveau discret, on vérifie:

$$P_{e} \left(\frac{\partial u}{\partial x}\right)_{i}=\left(\frac{\partial^{2}u}{\partial x^{2}}\right)_{i}$$

soit

$$P_{e}\frac{u_{2}-u_{0}}{2dx}=\frac{u_{2}-2u_{1}+u_{0}}{dx^{2}}$$ (3.41)

Avec la discrétisation choisie, le membre de gauche est constant et indépendant de $u_{1}$ . Le membre de droite est la valeur approchée de la dérivée seconde au noeud 1 de la courbe passant par les 3 points $[x_{0},x_{1},x_{2}]$ . La valeur de $u_{1}$ est donc telle que la courbure est égale à la valeur du membre de gauche. Cette valeur augmentant avec le nombre de Péclet, la courbure augmente et la valeur de $u_{1}$ s'éloigne par valeur négative de la valeur moyenne $\frac{u_{2}+u_{0}}{2}$ (qui correspond à une courbure nulle). D'après (3.41), la valeur de $u_{1}$ s'écrit:

$$u_{1}=\frac{u_{2}+u_{0}}{2}-P_{e}\frac{u_{2}-u_{0}}{4dx}$$

Cette valeur de $u_{1}$ devient inférieur à $u_{0}$ dès que $P_{e}dx>2$ , qui est justement la condition (3.40) d'apparition des oscillations numériques.

Ces oscillations numériques apparaissent dans la couche limite thermique en présence de forte courbure. Elles indiquent que le maillage est insuffisant pour décrire précisément la courbure de la couche limite.

3.4.4 Expérimentation numérique avec Matlab

programme Matlab 3.4.4: Résolution par différences finies de l'équation 3.35

% resolution de l'equation de convection diffusion
%clear;
% parametres
N=11; Peh=2.1
L=1.0; dx=1/(N-1); Pe=Peh/dx;
X=[0:dx:L];
% definition du schema D.F. operateur
a=[-(Peh+2)/2 2  (Peh-2)/2]
% construction des matrices tridiag
e=ones(1,N);
A=[a(1)*e; a(2)*e; a(3)*e];
B=zeros(1,N);
% C.L dirichlet en 1 et N
A(:,1)=0; A(2,1)=1;  B(1)=0;
A(:,N)=0; A(2,N)=1;  B(N)=1;
% resolution
U=tridiag(A,B);
% trace solution finale
clf; figure(1);
Xe=[0:1/100:L];
Ue=(exp(Pe*Xe)-1)/(exp(Pe)-1);
plot(X,U,Xe,Ue);

Nous avons programmé le schéma centré avec Matlab (programme 3.4.4). Sur la figure (3.20), on a comparée la solution numérique calculée avec $N=11$ points et la solution exacte pour différentes valeurs de $P_{e}^{h}$ . On constate que pour $P_{e}^{h}\leq2$ , la solution numérique est proche de la solution exacte et ne présente aucune oscillation. Par contre pour $P_{e}^{h}>2$ , la solution numérique présente des oscillations, qui montre que le maillage n'est plus assez fin pour capter la solution exacte.

Figure 3.20: Comparaison de la solution exacte et de la solution numérique avec le schéma centré 3.37 pour différentes valeurs de $P_{e}^{h}$

[cas $P_e^h<2$ avec $N=11$ ]

[cas $P_e^h>2$ avec $N=11$ ]

3.4.5 Discrétisation avec un schéma décentré

Pour éviter les oscillations numériques, on peut utiliser une discrétisation décentrée du terme de convection, comme au paragraphe (3.3.7). On obtiens ainsi un schéma décentré (ou upwind), qui s'écrit pour $V>0$ :

$$P_{e} \frac{u_{i}-u_{i-1}}{dx}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{dx^{2}}$$ (3.42)

L'analyse de l'erreur de troncature s'effectue avec Maple comme pour le schéma centré, et on obtiens l'erreur de troncature suivante:

$$ErrT=-\frac{P_{e}}{2}\frac{\partial^{2}u}{\partial x^{2}}dx+O(dx^{2})$$

Le schéma décentré (3.42) est donc bien consistant à l'équation (3.35) et d'ordre 1. Cette erreur de troncature montre aussi que ce schéma décentré 3.42 est équivalent à l'équation de convection diffusion suivante:

$$P_{e} \frac{\partial u}{\partial x}-(1+\frac{P_{e}^{h}}{2})\frac{\partial^{2}u}{\partial x^{2}}=0$$ (3.43)

que l'on peut reécrire sous forme dimensionnelle:

$$V \frac{\partial u}{\partial x}-(\kappa+\frac{Vdx}{2})\frac{\partial^{2}u}{\partial x^{2}}=0$$

Le schéma décentré introduit donc une diffusion numérique supplémentaire $\frac{Vdx}{2}$ par rapport à la diffusion physique $\kappa$ . Ce schéma décentré est donc sur-dissipatif. En suivant la même démarche qu'au paragraphe 3.4.3, on peut déterminer la solution analytique de l'équation aux différences (3.42). Pour cela on écrit le schéma aux différences finies (3.42) sous la forme d'une relation de récurrence à 3 niveaux:

$$-u_{i+1}+(2+P_{e}^{h})u_{i}-(1+P_{e}^{h}) u_{i-1}=0$$

La solution générale de cette relation de récurrence s'écrit en fonction des racines $\Lambda=[1,P_{e}^{h}+1]$ de l'équation caractéristique associée sous la forme:

$$u_{i}=C_{1}(\Lambda_{1})^{i}+C_{2}(\Lambda_{2})^{i}$$

Les constantes $C_{1}$ et $C_{2}$ sont déterminées par les conditions aux limites $u_{0}=0$ et $u_{N}=1$ . La solution exacte de l'équation aux différences (3.42) s'écrit:

$$u_{i}=\frac{\left(P_{e}^{h}+1\right)^{i}-1}{\left(P_{e}^{h}+1\right)^{N}-1}$$ (3.44)

Cette solution numérique converge vers la solution exacte aux noeuds du maillage (3.39), mais uniquement à l'ordre 1 (comme prévu par l'erreur de troncature). Un développement limité de l'écart $u(\frac{i}{N})-u_{i}$ montre que cet écart est d'ordre 1:

$$u(\frac{i}{N})-u_{i}=O(dx)$$

contrairement à la solution numérique du schéma centré (3.38) qui converge à l'ordre 2. Si on compare maintenant cette solution numérique $u_{i}$ à la solution exacte $u_{eq}$ de l'équation équivalente (3.43) fournie par l'analyse de troncature:

$$u_{eq}=\frac{1-e^{\frac{Pe}{1+\frac{P_{e}}{2}dx}x}}{1-e^{\frac{Pe}{1+\frac{P_{e}}{2}dx}}}$$

on montre que l'écart $u_{eq}(\frac{i}{N})-u_{i}$ est alors d'ordre 2:

$$u_{eq}(\frac{i}{N})-u_{i}=O(dx^{2})$$

Le nombre de Péclet effectif de la solution numérique est donc plus petit que le nombre de Péclet exacte et vaut $P_{eff}=\frac{Pe}{1+P_{e}dx/2}$ . Le schéma décentré a donc introduit de la diffusion numérique. Grâce à cette diffusion numérique, la solution $u_{i}$ reste toujours comprise entre les deux valeurs des conditions limites (0 et 1), sans oscillations.

3.4.6 Expérimentation numérique avec Matlab

programme Matlab 3.4.6: Résolution par différences finies décentrées de l'équation 3.35

% resolution de l'equation de convection diffusion
%clear;
% parametres
N=11; Peh=1/2
L=1.0; dx=1/(N-1); Pe=Peh/dx;
X=[0:dx:L];
% definition du schema D.F. operateur
a=[-(Peh+1) 2+Peh  -1]
% construction des matrices tridiag
e=ones(1,N);
A=[a(1)*e; a(2)*e; a(3)*e];
B=zeros(1,N);
% C.L dirichlet en 1 et N
A(:,1)=0; A(2,1)=1;  B(1)=0;
A(:,N)=0; A(2,N)=1;  B(N)=1;
% resolution
U=tridiag(A,B);
% trace solution finale
clf; figure(1);
Xe=[0:1/100:L];
Ue=(exp(Pe*Xe)-1)/(exp(Pe)-1);
plot(X,U,Xe,Ue);

Nous avons programmé le schéma décentré avec Matlab (programme 3.4.6). Sur la figure 3.21, on a comparée la solution numérique calculée avec $N=11$ points et la solution exacte pour différentes valeurs de $P_{e}^{h}$ .

En comparant avec la solution numérique obtenue avec le schéma centré (figure 3.20), on constate que pour $P_{e}^{h}\le2$ , la solution décentrée est moins précise que la solution centrée. Par contre pour $P_{e}^{h}>2$ , la solution décentrée reste monotone, mais en restant assez loin de la solution exacte.

Figure 3.21: Comparaison de la solution exacte et de la solution numérique avec le schéma décentré 3.42

[$P_e^h<2$ avec $N=11$ ]

[$P_e^h>2$ ]

Enfin sur la figure (3.22), on a comparé la solution numérique ( $N=11,$ $P_{e}^{h}=3$ ) à la solution exacte calculée avec le nombre de Péclet $P_{e}$ et le nombre de Péclet équivalent $P_{eff}$ . On constate que la solution numérique décentrée est beaucoup plus proche de la solution avec le nombre de Péclet éffectif $P_{eff}$ , qu'avec le nombre de Péclet exacte $P_{e}$ , ce qui confirme l'analyse précédente.

Figure 3.22: Comparaison de la solution numérique avec la solution exacte calculée avec le nombre de Péclet éffectif $P_{eff}=\frac {P_{e}}{1+P_{e}dx/2}$ et le nombre de Péclet exacte $P_{e}$