On veut calculer la répartition de température transversale dans une
flamme (figure 3.32). Pour cela on utilise un modèle simple,
en supposant que dans la direction transverse (notée
) on a des
échanges de chaleur uniquement par diffusion et par rayonnement.
Avec ces hypothèses, l'équation de conservation de l'énergie à l'état
stationnaire dans la direction
s'écrit:
Dans ce modèle, l'énergie
produite dans la flamme de largeur
par réaction chimique est diffusée par conduction et rayonné
vers l'air extérieur à température
. Pour le rayonnement,
nous avons adopté un modèle simple de rayonnement de corps noir proportionnel
à
avec une constante de rayonnement
. Les variations
de température entre la flamme et l'extérieur étant importantes (avec
un rapport de l'ordre de 3 à 4), le coefficient de conduction
dépend de la température à travers une loi en puissance
.
Enfin pour le terme source, nous choisirons un terme de réaction constant
dans la flamme et nul à l'extérieur.
En supposant que la flamme est symétrique, les conditions aux limites s'écrivent:
étant une distance grande par rapport à l'épaisseur
de la flamme.
En posant
, le problème modèle associé s'écrit
alors (en choisissant
et
):
avec
et
.
Le problème (3.67) étant non linéaire, il n'existe pas de
solution analytique simple. Cependant dans le cas de faibles variations
de
,
on peut écrire un problème linéaire équivalent (en posant
):
C'est le problème de diffusion linéaire avec source étudié au paragraphe c3chalsource.
Le programme Maple (3.7.2) calcule la solution générale de ce problème (3.68) en utilisant la fonction dsolve (ligne 7) et en déterminant les constantes d'intégration à l'aide des conditions aux limites (lignes 8 à 9).
> restart; # Diffusion dans une flamme: modele lineaire > T0:=1; > -diff(kappa*diff(T(x),x),x)+alpha*(T(x)-T0)=Q(x);eq:=%: > Q:=x->beta*Heaviside(delta-x); delta:=2/10; # Resolution et prise en compte des C.L. > dsolve(eq,T(x)); T1:=unapply(rhs(%),x); > D(T1)(0)=0; T1(1)=T0; solve({%,%%},{_C1,_C2});rel:=%: > Tex:=unapply(subs(rel,T1(x)),x); # Trace du resultat > subs(kappa=1/100,alpha=10,beta=30,T0=1,Tex(x));plot(%,x=0..1);
On a tracé sur la figure (3.33) cette solution analytique
pour
et différentes valeurs de
pour montrer
l'importance du rayonnement qui a tendance à raidir le front de température.
Pour résoudre numériquement le problème non linéaire (3.67),
on le discrétise tout d'abord par différences finies centrées avec
un maillage de pas
sous la forme:
A ces équations on ajoute les 2 conditions aux limites:
Pour résoudre numériquement ces équations, on transforme le problème en un problème équivalent de point fixe (i.e. ayant la même solution que 3.70):
A partir de ce système (3.71), on construit une suite de
valeurs
telle que:
Si la suite
converge, elle converge vers
un point fixe (i.e. une solution) de (3.71) et donc vers
la solution du problème non linéaire initial (3.70). La condition
de convergence de la suite de point fixe (3.72) est donnée
par le théorème classique du point fixe:
Une première approche consiste à considérer la solution du problème initial (3.67) comme solution stationnaire du problème instationnaire associé:
On discrétise cette équation avec du schéma différences finies explicite en temps:
Ce schéma explicite s'écrit s'écrit sous la forme d'un problème de point fixe:
avec
La matrice jacobienne
de
est la matrice tridiagonale
d'ordre
suivante:
avec
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
Par définition, cette matrice Jacobienne relie une petite variation
de la solution à l'étape
à la variation correspondante
à l'étape
:
En interprétant
comme une perturbation, la relation
précédente montre que la matrice
est la matrice
d'amplification du schéma différences finies explicite (3.74),
i.e. une perturbation
de la solution de (3.73)
vérifie
Or pour que le schéma (3.74) soit stable, il faut que les
valeurs propres de la matrice
soient en module inférieures
à 1. Pour notre schéma, la condition de convergence de la méthode
de point fixe est en faite une condition de stabilité sur le schéma
explicite associé.
Pour pouvoir faire l'analyse simplement, on simplifie les coefficients
de la matrice
en négligeant les termes en
devant ceux en
:
L'équation sur la perturbation
s'écrit alors:
Cette équation est analogue à l'équation de la perturbation du schéma
explicite pour le problème de convection linéaire étudié au paragraphe
(c3chalsource) avec les paramètres
et
.
L'étude de stabilité à été faîte au paragraphe 3.2.4 pour ce schéma explicite linéaire, et on a trouvé une condition de stabilité donnée par la relation (3.7).
Par analogie, la condition de stabilité pour le schéma explicite non linéaire (3.75) s'écrit donc:
Contrairement au cas linéaire, cette condition dépend de la solution
. En notant
, la valeur maximale
de la solution à l'étape
, on choisit le pas en temps
à
l'étape
tel que:
Le schéma explicite converge pour de très petites valeurs du
pas en temps, vérifiant la condition (3.76). Cette condition
devient très sévère lorsque le nombre de points de calcul
augmente
(i.e. lorsque
diminue).
Pour améliorer la convergence, on peut utiliser un schéma implicite
linéarisé, basé sur le schéma implicite linéaire (3.5) avec
.
Ce schéma s'écrit aussi sous la forme d'une itération de point fixe:
où
est la matrice tridiagonale d'ordre
suivante:
avec
La convergence de l'itération de point fixe (3.78) est lié
à la stabilité du schéma (3.77). La matrice Jacobienne
de l'itération de point fixe s'écrit:
Dans le cas linéaire, la matrice
On peut donc en conclure que le schéma implicite linéarisé converge pour des petits pas en temps. Cependant la limite de convergence est beaucoup plus grande que pour le schéma explicite, et on a un schéma plus efficace.
Pour rechercher les racines des équations non linéaires (3.70), on peut utiliser la méthode de Newton-Raphson, qui consiste à construite la suite itérative de point fixe suivante:
est la matrice jacobienne des fonctions
:
calculée à l'itération
. Cette relation s'écrit sous la forme
matricielle:
A chaque itération de Newton, il faut résoudre un système linéaire
d'ordre
.
Dans notre cas, la matrice Jacobienne
est tridiagonale
et s'écrit:
avec:
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
|
![]() |
![]() |
![]() |
Comme pour le schéma explicite, on peut simplifier la calcul de ces
coefficients en en négligeant les termes en
devant ceux en
, ce qui donne:
La méthode de Newton converge au voisinage de la racine
.
En effet la fonction de point fixe
associée s'écrit:
et sa matrice Jacobienne
vaut:
![]() |
![]() |
![]() |
|
![]() |
![]() |
Cette matrice Jacobienne est identiquement nulle en
car
. Donc sa norme vaut zéro à la racine,
et la méthode de Newton converge au voisinage de cette racine.
Cette méthode de point fixe converge au voisinage de la racine (à condition que la racine soit une racine simple), et possède une vitesse de convergence quadratique. C'est la méthode la plus efficace.
Les trois méthodes précédentes ont été implémentées sous Matlab, et comparées sur deux cas tests:
Le programme Matlab (3.7.7) calcule la solution du problème
non linéaire (3.67) à l'aide du schéma explicite (3.74).
Le schéma est écrit pour les noeuds internes à la ligne 24, et les
conditions aux limites sont imposées pour les noeuds
(ligne
26) et
(ligne 27). Pour tester la convergence vers la solution
du problème non linéaire discret (3.70), on calcul à chaque
itération en temps
la norme euclidienne du résidu divisée par
le nombre de points
(ligne 36):
Cette erreur correspond à une approximation de la norme intégrale
du carré du résidu
, qui tend vers zéro lorsque la suite converge.
La solution exacte est obtenue lorsque le résidu est nul, i.e. :
.
Numériquement, si ce résidu est inférieure à une valeur
fixé (
), on arrête les itérations (ligne 32) et
on considère que la solution numérique du problème non linéaire est
la solution calculée à cette itération.
% resolution probleme non lineaire de diffusion % schema explicite clear; N=51; L=1; gamma=0.9; dx=L/(N-1); X=[0:dx:L]; % parametre K0=0.01; sigma=1; beta=300; % nds internes I=2:N-1; % second membre delta=0.2; Q=(X<delta)*beta; % loi K K=inline('t.^2'); % C.I. Un=ones(1,N); dUn=zeros(1,N); % iterations nitmax=5000; epsilon=1.0e-05; for it=1:nitmax % calcul du dt Umax=max(Un); dt=gamma*2/(4*sigma*Umax^3+4*K0*K(Umax)/dx^2); % nds internes K12=0.5*(K(Un(1:N-1))+K(Un(2:N))); dUn1(I)=K0/dx^2*(K12(I-1).*(Un(I-1)-Un(I))+K12(I).*(Un(I+1)-Un(I)))-sigma*(Un(I).^4-1.0)+Q(I); % C.L en 0 dUn1(1)=K0/dx^2*(2*K12(1)*(Un(2)-Un(1)))-sigma*(Un(1).^4-1.0)+Q(1); dUn1(N)=0; % erreur Err=norm(dUn1)/sqrt(N); % fin Un=Un+dt*dUn1; if (Err<epsilon) break; end; end;
Le pas en temps est choisie proportionnelle à la limite de stabilité (3.76) (ligne 21).
Pour le cas 2, nous avons testé la convergence du schéma en fonction
de
et on a tracé le résultat sur la figure (3.35).
On constate que pour
, le schéma ne converge pas
ou diverge
. Pour
, la convergence est d'autant
plus grande que
est proche de 1. On a choisi une valeur
de
proche de 1 pour avoir la meilleur convergence
(ligne 4). Cela confirme l'analyse de stabilité du paragraphe 3.7.4.
Le programme Matlab (3.7.7) calcule la solution du problème non linéaire (3.67) à l'aide du schéma implicite linéarisé (3.77). La structure du programme est la même que précédemment. A chaque itération en temps, on calcule le second membre et la matrice tri-diagonale du système linéaire (3.78) (lignes 29 à 32), que l'on résout à l'aide de la fonction tridiag.
% resolution probleme non lineaire de diffusion % schema implicite linearise clear; N=51; L=1; gamma=10; dx=L/(N-1);X=[0:dx:L]; % parametre K0=0.01; sigma=1; beta=300; % nds internes I=2:N-1; % second membre delta=0.2; Q=(X<delta)*beta; % loi K K=inline('t.^2'); % C.I. Un=ones(1,N); dUn=zeros(1,N); Fn=zeros(1,N); % iterations nitmax=1000; epsilon=1.0e-5; for it=1:nitmax % calcul du dt Umax=max(Un); dt=gamma*2/(4*sigma*Umax^3+4*K0*K(Umax)/dx^2); % calcul du 2nd membre Fn(I)=Un(I)+dt*(Q(I)+sigma); % C.L. Fn(1)=Un(1)+dt*(Q(1)+sigma); Fn(N)=1; % calcul de la matrice 3Diag Kn12=0.5*(K(Un(1:N-1))+K(Un(2:N))); A=[1+dt*(K0/dx^2*(2*Kn12(1))+sigma*Un(1)^3),1+dt*(K0/dx^2*(Kn12(I-1)+Kn12(I))+sigma*Un(I).^3),1]; B=[-dt*K0/dx^2*(2*Kn12(1)), -dt*K0/dx^2*(Kn12(I)),0]; C=[0,-dt*K0/dx^2*(Kn12(I-1)),0]; J=[C;A;B]; % resolution Un1=tridiag(J,Fn); % erreur Err=norm((Un1-Un)/dt)/sqrt(N); % fin Un=Un1; if (Err<epsilon) break; end; end;
Le pas en temps est choisie proportionnelle à la limite de stabilité
(3.76) (ligne 21) suivant la relation (
).
Pour le cas 2, nous avons testé la convergence du schéma en fonction
de
et on a tracé le résultat sur la figure (3.35).
On constate que contrairement au schéma implicite linéaire, le schéma
n'est pas inconditionnellement stable, i.e. on ne peut pas choisir
un pas en temps
arbitrairement grand, i.e. pour
,
le schéma ne converge plus sans toutefois diverger. Ainsi pour
,
le résidu oscille. Cette oscillation du résidu est due au choix d'un
pas en temps fonction de la valeur maximum de la solution (
).
Lorsque le schéma se met à diverger, la valeur
augmente
et le pas en temps
diminue. Le pas en temps ayant diminué, le
résidu diminue. A l'itération suivante le pas en temps augmente donc,
ce qui fait croître à nouveau le résidu, et ainsi de suite. D'où la
série d'oscillations sur le résidu que l'on observe sur la figure
(3.35) pour
. On constate aussi que pour
,
la vitesse de convergence diminue lorsque
décroit. Pour
on note aussi que la vitesse de convergence est beaucoup
plus grande que dans le cas du schéma explicite.
Le schéma implicite linéarisé a donc une limite de stabilité
, mais aussi une vitesse de convergence beaucoup plus
grande que le schéma explicite.
Le programme Matlab (3.7.7) calcule la solution du problème
non linéaire (3.67) à l'aide d'une méthode de Newton (3.79).
La structure du programme est la même que précédemment. Pour le calcul
de la matrice Jacobienne, on a utilisé l'expression simplifiée (3.81).
On note que l'utilisation d'une matrice jacobienne approchée influe
uniquement sur la vitesse de convergence et non sur la solution du
problème non linéaire, puisque lorsque la suite de Newton
(3.80) converge elle vérifie l'équation non linéaire
.
% resolution probleme non lineaire de diffusion % schema Newton Ralphson clear; N=51; L=1; dx=L/(N-1); X=[0:dx:L]; % parametre K0=0.01; sigma=1; beta=300; % nds internes I=2:N-1; % second membre delta=0.2; Q=(X<delta)*beta; % loi K K=inline('t.^2'); % C.I. Un=ones(1,N); dUn=zeros(1,N); Fn=zeros(1,N); % iterations nitmax=1000; epsilon=1.0e-5; for it=1:nitmax % calcul du 2nd membre -Fi % nds internes Kn12=0.5*(K(Un(1:N-1))+K(Un(2:N))); Fn(I)=K0/dx^2*(Kn12(I-1).*(Un(I-1)-Un(I))+Kn12(I).*(Un(I+1)-Un(I)))-sigma*(Un(I).^4-1.0)+Q(I); % C.L en 0 Fn(1)=K0/dx^2*(Kn12(1)*(Un(2)-Un(1))+Kn12(1)*(Un(2)-Un(1)))-sigma*(Un(1).^4-1.0)+Q(1); Fn(N)=0; % calcul de la matrice Jacobienne A=[K0/dx^2*(2*Kn12(1))+4*sigma*Un(1)^3,K0/dx^2*(Kn12(I-1)+Kn12(I))+4*sigma*Un(I).^3,1]; B=[-K0/dx^2*(2*Kn12(1)), -K0/dx^2*(Kn12(I)),0]; C=[0,-K0/dx^2*(Kn12(I-1)),0]; J=[C;A;B]; % resolution dUn=tridiag(J,Fn); % erreur Err=norm(Fn)/sqrt(N); % fin Un=Un+dUn; if (Err<epsilon) break; end; end;
On a comparé la vitesse de convergence de la méthode de Newton avec celle du schéma implicite linéarisé et du schéma explicite pour les deux cas tests (figure 3.37).
Sur cette figure, on note que la méthode de Newton est la méthode
qui converge le plus rapidement, et le schéma explicite celle qui
converge le plus lentement. Ainsi pour le cas 2 (le plus difficile),
il faut uniquement
itérations de Newton pour atteindre un résidu
de
, alors qu'il en faut
pour le schéma implicite
linéarisé, soit
fois plus, et
pour le schéma explicite,
soit
fois plus.