Méthodes Numériques pour la Mécanique

POLYCOPE DE COURS+TD - M1 MEPP

Olivier BOTELLA
LEMTA UMR CNRS 7563
2 avenue de la Forêt de Haye - TSA 60604 - 54518 VANDOEUVRE
Email : [Link]@[Link]
 Partie II

Différences finies
3 Introduction naïve aux différences finies

3.1 Exemple : fléchissement d’une poutre

Fig. 3.1. Fléchissement d’une poutre.

Considérons le problème représenté sur la figure 3.1, et que l’on retrouve dans de nombreux bouquins
(e.g. [2, 10]) : une poutre de longueur l = 1, étirée selon son axe par une force P et simplement appuyée à ses
extrémités 0 et 1, est soumise à une charge transversale f (x) par unité de longueur dx. La théorie des poutres
nous dit que le déplacement u(x) au point x est solution du problème aux limites :

−u00 (x) + c(x) u(x) = f (x), dans Ω =]0, 1[, (3.1a)

u(0) = 0, u(1) = 0. (3.1b)

avec c(x) = P/EI(x), où E est le module d’Young de la poutre et I(x) le moment principal d’inertie de la
section de la poutre au point x.
Sauf dans de très rares situations, on ne connait pas de solution analytique à un problème aux limites.
Les méthodes d’approximation numériques, et en particulier la méthode des différences finies, permettent
d’approcher les valeurs de la solution d’aussi près qu’on veut.

3.2 Résolution par la méthode des différences finies

Afin de discrétiser le problème (3.1), on définit l’ensemble
1
Ωh = {xi = (i − 1)h, i = 1, . . . , N }, avec h = , (3.2)
N −1
comme le maillage de Ω dont le pas de maillage, ou taille de la maille, est h. Ce maillage est représenté sur la
figure 3.2 ; les points xi sont appelés les nœuds du maillage.

Fig. 3.2. Maillage utilisé pour le problème de la poutre. Le symbole ◦ désigne les nœuds intérieurs et • désigne les
nœuds de la frontière.

La valeur du©pas h est destinée ª à tendre vers 0. La méthode des différences finies est une façon d’obtenir une
approximation ui , i = 1, . . . , N de la solution u(xi ) aux nœuds xi . D’après les conditions aux limites (3.1b),
on a u1 = uN = 0, et il reste à chercher le vecteur
 
u2
 
U =  ...  ∈ RN −2 , (3.3)
uN −1
24 3 Introduction naïve aux différences finies

tel que les ui soit une “bonne” approximation de la solution u(xi ) aux points intérieurs du maillage. De plus, la
qualité de l’approximation aux différences finies sera d’autant meilleure que h → 0, i.e. que le nombre de points
N du maillage est grand.
Pour cela, on va établir une approximation par différences finies de l’opérateur L = −d2 /dx2 + c(x) par
application de la formule de Taylor :
h2 00 h3 000 h4 (4) h5 (5)
u(xi + h) = u(xi ) + hu0 (xi ) + u (xi ) + u (xi ) + u (xi ) + u (xi + θi+ h), (3.4a)
2! 3! 4! 5!
h2 h3 000 h4 (4) h5 (5)
u(xi − h) = u(xi ) − hu0 (xi ) + u00 (xi ) − u (xi ) + u (xi ) − u (xi + θi− h), (3.4b)
2! 3! 4! 5!
avec θi+ , θi− ∈]0, 1[. En combinant ces deux expressions pour éliminer u0 (xi ), on obtient l’expression :
u(xi + h) − 2 u(xi ) + u(xi − h)
u00 (xi ) = 2
+ h2 R(xi ) , (3.5a)
| h{z } | {z }
formule aux différences finies Reste ou résidu

tel que

−1 (4) h
R(xi ) = u + (u(5) (xi + θi+ h) − u(5) (xi − θi− h)). (3.5b)
12 i 120
On rappelle les définitions suivantes :
Définition 3.1 (Grand “Oh” et petit “oh”).
On dit qu’une fonction f est telle que f = O(g) s’il existe K ≥ 0 tel que
|f |
≤K lorsque h → 0.
|g|
D’autre part, f = o(g) si
|f |
→0 lorsque h → 0.
|g|
On utilisera essentiellement les grands “Oh” dans la suite du cours.
D’après ces définitions, R(xi ) = −u(4) /12 + O(h), donc le résidu est une fonction bornée. On pourrait aussi
écrire que
u(xi + h) − 2 u(xi ) + u(xi − h)
u00 (xi ) = + O(h2 ). (3.6)
h2
Il ne reste plus qu’à reporter la formule (3.5a) dans l’équation différentielle (3.1a). En tenant compte des
conditions aux limites (3.1b), soit u1 = uN = 0, on obtient le système linéaire

 2u2 − u3

 + c2 u2 = f2 + h2 R2 ,

 h2


−ui−1 + 2ui − ui+1
+ ci ui = f i + h2 R i , pour i = 3, . . . , N − 2, (3.7)

 h2



 −uN −2 + 2uN −1 + cN −1 uN −1 = fN −1 + h2 RN −1 ,

h2
dans lequel on a posé ci = c(xi ), fi = f (xi ), et Ri = R(xi ). Ce système s’écrit sous la forme matricielle

Ah U = F + τh , (3.8a)

où Ah est la matrice tridiagonale définie par

 
2 + c2 h2 −1
 
 −1 2 + c3 h2 −1 
 
 
 · · · 
 
1  −1 2 + ci h2 −1

Ah = 2  , (3.8b)
h 


 · · · 
 
 −1 2 + cN −2 h2 −1 
 
 
−1 2 + cN −1 h2
 3.2 Résolution par la méthode des différences finies 25

et en posant
   
f2 R2
   
F =  ...  , τh = h2  ...  . (3.8c)
fN −1 RN −1
Le vecteur τh , qui tend vers 0 lorsque h → 0, est naturellement conduit à être négligé : on aura donc remplacé
la résolution du problème aux limites linéaire (3.1) par le système algébrique linéaire
Ah U = F. (3.9)
Ce système correspond à la discrétisation du problème (3.1) par la méthode des différences finies. Il est souvent
appelé schéma aux différences finies.
Le vecteur τh est appelé résidu ou erreur de troncature du schéma. En pratique, il est toujours négligé dans
l’expression du schéma. Néanmoins, il donne des renseignements sur la qualité de l’approximation obtenue : ici
τh = O(h2 ), on dit alors que le schéma est précis au second ordre. Ces notions seront énoncées de manière plus
formelle dans le paragraphe 5.4.
Exercice 3.2 (Trigonométrie pénible mais résultat important). On considère la matrice tridiagonale de
taille N × N :  
2 −1
 
 −1 2 −1 
 
K= . . .
.. .. .. . (3.10)

 
 −1 2 −1 
−1 2
1. On pose h = 1/(N + 1). Montrez que1 les valeurs propres de K sont données par :
λk = 2 − 2 cos kπh, pour k = 1, . . . , N, (3.12a)

dont le vecteur propre associé est

T
xk = (sin jkπh, j = 1, . . . , N ) . (3.12b)
2. En remarquant que :
cos N πh = − cos πh,

déduire de la question précédente que le nombre de conditionnement de K s’écrit

1 + cos πh
χ(K) = ,
1 − cos πh
et que, lorsque h → 0,
4
χ(K) ' .
π 2 h2
3. Que pouvez-vous déduire de ces questions sur la performance des méthodes du gradient appliquées à la
résolution du système linéaire issu de l’approximation par différences finies du problème (3.1) avec c(x) ≡ 0
?
4. Question subsidiaire : plus généralement, on considère la matrice de taille N × N :
 
b c
a b c 
 
 · · · 
 
 a b c 
 
 · · · 
 
 a b c 
a b
1
On rappelle à ce sujet les formules
sin(a ± b) = sin a cos b ± cos a sin b,
cos(a ± b) = cos a cos b ∓ sin a sin b,
sin(a − b) + sin(a + b) = 2 sin a cos b.
26 3 Introduction naïve aux différences finies

Montrez que, pour k = 1, 2, . . . , N , la k ième valeur propre est :

¡ ¢1/2 kπ
λk = b − 2 a c cos λk , αk = ,
N +1
de vecteur propre associè :
³ ³ a ´j−1 ´
2
xk = sin jαk , j = 1, . . . , N.
c
4 Notions sur les équations aux dérivées partielles (EDP)

4.1 Les équations de bilan de la mécanique

On rappelle que Ω désigne le domaine physique (ou de calcul), qui est un sous-domaine de R2 ou R3 . On
désignera sa frontière par ∂Ω, Γ , ou S.
Le but de cette partie est d’établir la loi de bilan d’une propriété φ quelconque d’une particule de matière
fluide ou solide, par exemple φ = T sa température. On supposera que φ est définie par unité de masse, et ρ
désigne la masse volumique. La loi de bilan pour φ est traditionnellement exprimée pour un volume V ⊂ Ω
constitué des mêmes particules matérielles, et qui évolue dans Ω avec la vitesse v. On aura besoin de la relation
fondamentale suivante, énoncée sous la forme d’un théorème :
Théorème 4.1 (Théorème du transport de Reynolds).
Z Z µ ¶
d ∂ρφ
∀ V ⊂ Ω, ρφ dV = + ∇ · ρφv dV. (4.1)
dt V V ∂t
Dans le volume matériel V , la loi de conservation exprime que φ peut uniquement varier sous les effets d’échange
entre les particules matérielles et la présence de sources (ou puits) externes. La loi de bilan pour la quantité φ
peut alors s’exprimer comme :

Taux de variation = Echanges + Sources,

qui s’écrit sous forme mathématique :

Z Z Z
d
ρφ dV = t · n dS + q dV , (4.2)
dt V
| {z } | S {z } | V {z }
I1 I2 I3
tel que
• I1 représente la variation de la grandeur physique au cours du temps.
• I2 désigne le flux au travers de la surface S, dues aux forces de pression ou aux contraintes agissant sur le
volume V , qui sont représentées par le vecteur t.
• I3 représente la production ou dissipation à l’intérieur de V , dues à des forces extérieures q (e.g. gravité
q ≡ ρg, forces centrifuges ou électromagnétiques, . . . ) agissant sur V .
En utilisant le théorème du transport 4.1, on obtient la forme intégrale d’une loi de bilan :
Z Z Z
∂ρφ
+ (ρφv − t) · n dS = q dV. (4.3)
V ∂t S V

De plus, il est possible de transformer l’intégrale de surface en intégrale de volume en utilisant le théorème de
Gauss :
Théorème 4.2 (Théorème de la divergence de Gauss).
Pour tout volume V dont la surface S est fermée, nous avons pour tout vecteur u :
Z Z
∇ · u dV = u · n dS, (4.4)
V S

où n désigne la normale extérieure unitaire de S.

On obtient alors l’expression : Z · ¸ Z
∂ρφ
+ ∇ · (ρφv − t) dV = q dV, (4.5)
V ∂t V

d’ou on déduit la forme locale1 d’une loi de bilan, valable en tout point x de V :
1
Cette forme est aussi qualifiée de “conservative”, pour mettre en avant le lien avec la forme globale (4.3) qui est
exprimée par le théorème de la divergence.
28 4 Notions sur les équations aux dérivées partielles (EDP)

∂ρφ
+ ∇ · (ρφv − t) = q, (4.6)
∂t
qui prends alors la forme d’une équation aux dérivées partielles (EDP).
Les formes intégrale (4.3) ou locale (4.6) représentent les formes génériques des équations de bilan de la
mécanique.
Exemple 4.3 (Equation de conservation de la masse). L’équation de continuité en mécanique des fluides s’obtient
en posant φ ≡ 1, t ≡ 0, et q ≡ 0 dans (4.6), soit :
∂ρ
+ ∇ · (ρv) = 0, (4.7)
∂t
et dans le cas où la masse volumique est constante, on retrouve la contrainte d’incompressibilité :

∇ · v = 0.

Exemple 4.4 (Equation de Navier-Stokes). Pour la conservation de la quantité de mouvement d’un fluide New-
tonien, on prend φ ≡ v dans (4.6) et, par exemple, q ≡ 0. Lorsque t prend la forme du tenseur de Cauchy :

t = −pI + τ , (4.8)

où I est l’identité et τ le tenseur des contraintes visqueuses défini en (11.5), on obtient l’équation de Navier-
Stokes
∂ρv
+ ∇ · ρvv = −∇p + ∇ · τ . (4.9)
∂t
Exemple 4.5 (Equation de température). Pour des phénomènes de diffusion de la chaleur, le tenseur t suit la loi
de Fourier
t = κ∇T,
où T est la température et κ la conductivité thermique. En prenant q ≡ 0 et φ ≡ cp T , tel que cp est la chaleur
spécifique à pression constante, on obtient l’équation
µ ¶
∂T
ρcp + ∇ · (vT ) = κ∇2 T, (4.10)
∂t
si l’on suppose que les divers coefficients prennent des valeurs constantes.

4.2 Une équation générique : l’équation d’advection-diffusion

Le principal objectif de la section précédente est de montrer que les équations fondamentales de la mécanique,
qu’elles concernent la masse, l’énergie, la quantité de mouvement, etc . . . , ont la même structure. Il est alors
possible de développer des méthodes numériques générales pour résoudre (4.3) ou (4.6). Précisons que la forme
locale (4.6) est plus particulièrement utilisée dans la méthode des différences finies. La forme intégrale (4.3) est
à l’origine de la méthode des volumes finis et, dans une certaine mesure, des éléments finis.
Afin de simplifier notre présentation des méthodes numériques, on considérera dans un premier temps
l’équation d’advection-diffusion2 multidimensionnelle :
∂ρφ
+ ∇ · (ρφv) = ∇ · (ν∇φ) + q, (4.11)
∂t
qui est un modèle linéaire de la loi de bilan locale (4.6) pour la loi de Fick t = ν∇φ, où ν est un coefficient de
diffusion de signe positif strictement. Le terme de gauche représente le transport de φ par advection le long de
l’écoulement, et le premier terme à droite représente le transport par diffusion.
Exercice 4.6. Determinez la forme globale équivalente de (4.11).
Afin de simplifier encore plus le problème, nous allons considérer la version unidimensionelle linéaire de cette
équation, en posant ρ ≡ 1, φ ≡ u(x, t) et v ≡ a, telle que a est une constante d’advection de signe quelconque,
et ν > 0 est supposé constant. L’equation d’advection prend alors la forme :

∂u ∂u ∂2u
+ a − ν = q. (4.12)
∂t ∂x}
| {z ∂x2}
| {z
Advection Diffusion
2
Elle est aussi connue sous le non d’équation de convection-diffusion. A ne pas confondre avec l’équation de réaction-
diffusion utilisée en combustion.
 4.2 Une équation générique : l’équation d’advection-diffusion 29

Une grande partie de la dynamique des solutions de problèmes de mécanique est présente dans cette simple
équation linéaire. Par exemple, cherchons une solution de l’equation (4.12), sans terme de force (q ≡ 0), sous la
forme du mode de Fourier
b(t) e i kx ,
u(x, t) = u i2 = −1. (4.13)
L’équation devient
db
u
u + νk 2 u
+ i akb b = 0,
dt
soit ¡ ¢
2
u b0 e− aik + νk t ,
b(t) = u
où u
b0 est une condition initiale donnée. On obtient finalement la solution
2
b0 e|−νk
u(x, t) = u t ik (x − at) ,
{z } e| {z } (4.14)
I II

telle que :
• I est un terme de diffusion ou dissipation, strictement décroissant en temps.
• II est un terme de propagation.

Fig. 4.1. Solution de l’équation d’advection-diffusion au temps initial (traits pleins) et au temps t > 0 (traits pointillés).

Le graphe de cette solution est représenté sur la figure 4.1. Deux cas particuliers sont à noter :
Remarque 4.7. Dans le cas a ≡ 0, Eq. (4.12) devient

∂u ∂2u
− ν 2 = 0, (4.15)
∂t ∂x
et la solution s’écrit
2
b0 e−νk t e ikx .
u(x, t) = u (4.16)
La caractéristique principale de cette solution est son amortissement exponentiel dans le temps, dû à la présence
de la viscosité ν > 0 :
lim u(x, t) = 0 quand t → ∞.
L’équation (4.15) est connue sous le nom de l’équation de la chaleur. Comme il sera vu dans la section suivante,
ce problème représente le prototype des équations paraboliques.
Remarque 4.8. Dans le cas ν = 0 et a > 0, Eq. (4.12) devient
∂u ∂u
+a = 0. (4.17)
∂t ∂x
Elle prend le nom d’équation d’advection, prototype des équations hyperboliques, et la solution s’écrit

b0 e ik (x − at) .
u(x, t) = u (4.18)

Cette solution prend la forme d’une onde qui se translate à la vitessse constante a, sans déformation ou amor-
tissement : phénomène de propagation.
30 4 Notions sur les équations aux dérivées partielles (EDP)

4.3 EDP : le point de vue mathématique

La résolution numérique des équations aux dérivées partielles nécessite, au préalable, un minimum de connais-
sance sur leurs propriétés analytiques. Dans la suite du cours, les EDP seront souvent représentées sous la forme
générique
Au = f,
où A désigne un opérateur aux derivées partielles quelconque, et f est un terme de force donné. La solution de
ce problème prends la forme
u(t, x) : [0, T ] × Ω → R.
On fait souvent la distinction entre la variable t qui est associée au temps et la variable spatiale x qui appartient
au domaine physique Ω, sous-ensemble de R2 ou R3 . La frontière de Ω est souvent notée Γ ou ∂Ω.

4.3.1 Classification des EDP

Les EDP sont mathématiquement classées en plusieurs groupes, ou classes. Les problèmes relevant de chacune de
ces classes possèdent des caractéristiques physiques et mathématiques communes. L’intérêt de cette classification
est que l’étude de problèmes types relevant de chacune de ces classes permet alors de développer des méthodes
numériques universelles.
Plus précisement, les EDP linéaires du second ordre sont regroupées en trois classes distinctes : parabolique,
elliptique, et hyperbolique. Les problèmes physiques qui présentent une évolution au cours du temps sont gou-
vernés par des équations paraboliques ou hyperboliques, dont les exemples types sont :
• L’équation de la chaleur,
∂u
− ν∇2 u = f, ν > 0, (4.19)
∂t
de solution u(t, x) : [0, T ] × Ω → R est représentative des problèmes paraboliques.
• L’équation des ondes,
∂2u ∂2u
2
− c2 2 = f, (4.20)
∂t ∂x
et l’équation d’advection,
∂u ∂u
−c = f, (4.21)
∂t ∂x
de solution
u(t, x) : [0, T ] × R → R,
sont des exemples importants de problèmes hyperboliques.
Les problèmes physiques stationnaires sont généralement décrits par des équations elliptiques. L’équation de
Laplace,
−∇2 u = f avec u(x) : Ω → R, (4.22)
est le prototype de cette classe de problèmes.

Remarque 4.9. L’équation d’advection-diffusion (4.12) est parabolique pour ν > 0 ; si ν = 0 elle est hyperbolique.
Pour plus d’info, cf. Chap. 9.

4.3.2 Conditions initiales et aux limites, problèmes bien posés

Afin de définir un problème aux limites bien posé (existence et unicité de la solution), il faut imposer des
conditions initiales et des conditions de bords à la solution.

Conditions aux limites

Les conditions que l’on impose sur la frontière ∂Ω sont de plusieurs types. Un premier type consiste à imposer
la valeur de u sur la paroi. Ce type de conditions aux limites est appelé conditions de Dirichlet :

u(t, · )|∂Ω = g(t, · ), ∀ t ∈ [0, T ], (4.23)

où g : [0, T ] × ∂Ω → R est une donnée du problème. Par exemple, si u désigne la vitesse d’un fluide, alors (4.23)
avec g ≡ 0 exprime la condition d’adhérence du fluide à la paroi.
Un autre type important consiste à imposer la composante normale du gradient de u à la frontière. Il s’agit
de conditions de Neumann :
 4.3 EDP : le point de vue mathématique 31
¯
∂u ¯
(t, ·, )¯¯ ≡ (∇u · n)|∂Ω = g(t, ·, ), ∀ t ∈ [0, T ], (4.24)
∂n ∂Ω

où n désigne la normale extérieure, et g est une fonction connue. Par exemple, si u désigne la température d’un
matériau, alors (4.24) avec g = 0 revient à imposer une condition de flux nul à la paroi.
Citons aussi les conditions de Robin, utilisées notamment en thermique :
¯
∂u ¯
(t, ·, )¯¯ + λ u(t, ·, )|∂Ω = g(t, ·, ), ∀ t ∈ [0, T ], (4.25)
∂n ∂Ω

où λ et g sont connues. Enfin, précisons qu’il est possible de combiner ces différents types de conditions aux
limites. Si ∂Ω se décompose en deux parties disjointes3 ∂Ω1 et ∂Ω2 , on peut alors imposer, par exemple, des
conditions de Dirichlet sur ∂Ω1 et des conditions de Neumann sur ∂Ω2 . Il s’agit alors de conditions mixtes de
Dirichlet-Neumann.
Remarque 4.10. Un peu de vocabulaire : on qualifie d’homogènes des conditions aux limites où la donnée g est
identiquement nulle.

Conditions initiales

Il faut aussi imposer des conditions initiales à un problème d’évolution. Le nombre de conditions initiales
est généralement égal à l’ordre de différentiation en temps de l’EDP. Pour l’équation de la chaleur (4.19) et
l’équation d’advection (4.21), une seule condition suffit :

u(t = 0, x) = u0 (x) ∀x ∈ Ω, (4.26)

où u0 (x) est une donnée du problème. Pour l’équation des ondes (4.20), il faut formuler une condition supplé-
mentaire sur ∂u/∂t à l’instant initial.

Problème bien posé

Pour que la solution d’un problème aux limites possède un sens physique, il est nécessaire que le problème soit
mathématiquement bien posé, c.à.d.
• la solution doit exister,
• la solution doit être unique,
• la solution doit être stable au sens de Hadamard : elle doit dépendre continument des données4 .
On ne peut arbitrairement imposer n’importe quel type de conditions aux limites. L’intuition physique, suivie
d’une analyse mathématique, doivent servir de guide pour choisir les conditions adéquates, conduisant à un
problème bien posé. Ces conditions dépendent essentiellement de la classe de l’EDP considérée.
Exemple 4.11. L’équation de Laplace (4.22) est bien posée avec :
1. les conditions de Dirichlet (4.23), ou
2. les conditions de Neumann (4.24) si les données du problème vérifient la condition de compatibilité :
Z Z
g dS + f dV = 0,
S V

qui s’obtient en appliquant le théorème de la divergence 4.2 à l’équation (4.22). Cette condition exprime
un principe physique de conservation : le flux au travers de la frontière s’équilibre avec la production à
l’intérieur de Ω. Dans ce cas, la solution est connue à une constante (additive) près.
Dans la suite du cours, nous admettrons que tous les problèmes aux limites que nous considérons sont bien
posés.

3
C’est-à-dire : ∂Ω = ∂Ω1 ∪ ∂Ω2 , ∂Ω1 ∩ ∂Ω2 = ∅.
4
C’est-à-dire que si on perturbe quelque peu les données du problème, l’écart entre la solution originale et la solution
perturbée doit rester commensurable aux perturbations, et ne peut s’amplifier de manière dramatique. Un système
chaotique est-il, a priori, stable au sens de Hadamard ?
5 Concepts de base des différences finies

5.1 Introduction
L’approximation par différences finies d’un problème aux limites comprend les étapes suivantes :
1. Générer un maillage du domaine physique Ω, par exemple Ωh = {(xi , yj )} où (xi , yj ) est un nœud en lequel
une approximation de la solution sera calculée.
2. Discrétiser le problème sur la grille calcul Ωh . Cela est effectué lors des étapes suivantes :
a) Substituer les dérivées de l’EDP par un schéma de différences finies. L’idée derrière cette approximation
vient directement de la définition d’une dérivée, p. ex. :

du u(xi + ∆x) − u(xi )

(xi ) = lim .
dx ∆x→0 ∆x
b) Implémenter les conditions aux limites du problème.
L’EDP, que nous supposons linéaire, est alors transformée en un système algébrique d’équations linéaires1 .
3. S’assurer de la validité du schéma utilisé. Cela suppose d’effectuer une analyse des erreurs du schéma
numérique.
4. Résoudre le système discret par une méthode directe (pivot de Gauss, . . . ), itérative (Jacobi, Gauss-Seidel,
Gradient conjugué, . . . ), ou multigrille (cf. [8]).
5. Implémenter le programme sur un ordinateur et le débugger !

5.2 Grilles de calcul pour les différences finies

5.2.1 Premières notions sur les maillages

Les méthodes de différences finies sont particulièrement adaptées aux géométries simples de R2 (ou R3 ), telles
que le domaine Ω peut s’écrire dans un système de coordonnées (x, y)

Ω = [ax , bx ] × [ay , by ].

Par exemple, si Ω est un rectangle, alors (x, y) sont les coordonnées cartésiennes, pour un domaine annulaire
alors (x, y) = (r, θ), les coordonnées polaires. Ce type de maillages est appelé maillage cartésien. Une discussion
plus approfondie sur les différents types de maillage utilisés en pratique sera effectuée dans le chapitre 10.
En deux dimensions (2D), le maillage ou grille de calcul Ωh est défini par

Ωh = {(xi , yj ), i = 1, . . . Nx ; j = 1, . . . , Ny } ⊂ Ω. (5.1)

Un exemple est représenté sur la figure 5.1. Le pas de maillage suivant la direction x est noté ∆xi , défini par

∆xi = xi − xi−1 .

Le maillage est uniforme dans la direction x si :

∀ i = 2, . . . , Nx , ∆xi = h = 1/(Nx − 1). (5.2)

De même, dans la direction y, le pas de maillage est noté ∆yj .

1
Une EDP non-linéaire est transformée en un système algébrique d’équations non-linéaires, à résoudre par méthodes
itératives de Newton et ses variantes, cf. [10, 9].
34 5 Concepts de base des différences finies

Fig. 5.1. Un exemple de maillage cartésien bidimensionnel. Le symbole ◦ désigne les nœuds intérieurs et • désigne les
nœuds de la frontière. Noter que le pas de maillage est non-uniforme dans les deux directions.

5.2.2 Notions sur la discrétisation temporelle

Lorsqu’on veut calculer des problèmes aux limites instationnaires, une dimension supplémentaire doit être prise
en considération : le temps. De même que le domaine spatial a été discrétisé par un maillage, le domaine
temporel, e.g. [0,T], doit lui aussi être représenté par une “grille” temporelle :
{tn = n∆t, n = 0, 1, . . . } , (5.3)
où ∆t représente le pas de temps.
La majeur différence entre le temps et l’espace est que la solution au temps (n + 1)∆t est uniquement
influencée par la solution aux instants précédents : l’information ne peut voyager dans le passé ! C’est pour cela
que la plupart des schémas de différences finies utilisent des méthodes classiques pour les problèmes aux valeurs
initiales2 des équations aux dérivées ordinaires (EDO), de la forme :
½ 0
y (t) = f (t, y(t)),
(5.4)
y(0) = y 0 .
Par la suite, nous utiliserons essentiellement des méthodes à un pas, telles la méthode d’Euler progressive,
y n+1 = y n + ∆tf n , (5.5a)

la méthode d’Euler rétrograde,

y n+1 = y n + ∆tf n+1 , (5.5b)

et la méthode du trapèze ou de Crank-Nicolson,

∆t ¡ n ¢
y n+1 = y n + f + f n+1 . (5.5c)
2
Les applications de ces méthodes aux différences finies sont présentées dans les paragraphes 5.4.3 et 6.3.2. Des
perspectives sur l’utilisation de schémas plus sophistiqués, comme les méthodes de Runge-Kutta, sont présentées
dans le paragraphe 6.2. Pour plus de renseignements sur les schémas pour EDO, consulter [9].
Certaines approximations par différences finies, comme les schémas ADI ou de Lax-Wendroff, n’ont pas de
lien avec les méthodes pour EDO. Ces schémas sont présentés dans les paragaphes 8.4 et 9.6 respectivement.

5.2.3 Notation compacte de la solution

Pour représenter la solution sur la grille de calcul, on utilise la notation compacte de u :

uni,j ≡ u(xi , yj , tn ), (5.6a)

et de ses dérivées, e.g. la dérivée première :

µ ¶n
∂u ∂u
≡ (xi , yj , tn ). (5.6b)
∂x i,j ∂x
2
aussi appelés problèmes de Cauchy.
 5.3 Construction de formules aux différences finies 35

5.3 Construction de formules aux différences finies

Le concept d’une formule de différences finies provient directement de la définition d’une dérivée, par exemple :
µ ¶
∂u u(xi + ∆x) − u(xi )
= lim , (5.7)
∂x i ∆x→0 ∆x

et leur construction repose sur l’utilisation de développements de Taylor du type

X hk
u(x + h) = u(k) (x). (5.8)
k!
k≥0

Dans le cas d’une fonction u(x, t) sur maillage uniforme de pas ∆x, nous pouvons écrire les développements en
espace
µ ¶n µ ¶n µ ¶n
∂u ∆x2 ∂ 2 u ∆x3 ∂ 3 u
uni+1 ≡ u(xi + ∆x, tn ) = uni + ∆x + + + O(∆x4 ), (5.9a)
∂x i 2 ∂x2 i 6 ∂x3 i
µ ¶n µ ¶n µ ¶n
n n ∂u ∆x2 ∂ 2 u ∆x3 ∂ 3 u
ui−1 ≡ u(xi − ∆x, tn ) = ui − ∆x + − + O(∆x4 ), (5.9b)
∂x i 2 ∂x2 i 6 ∂x3 i

mais aussi en temps

µ ¶n µ ¶n µ ¶n
∂u ∆t2 ∂2u ∆t3 ∂3u
un+1
i ≡ u(xi , tn + ∆t) = uni + ∆t + + + O(∆t4 ). (5.9c)
∂t i 2 ∂t2 i 6 ∂t3 i

Il est bien sûr possible de combiner des développements en espace et en temps dans une même formule. Ces
développements vont servir à construire des formules aux différences pour les opérateurs différentiels usuels.
Pour un même opérateur, il est possible de construire un grand nombre de schémas de différences finies. Ces
formules diffèrent essentiellement par les nœuds qui sont utilisés, et leur erreur de troncature, i.e. le premier
terme négligé dans le développement de Taylor. Les formules les utilisées pour les opérateurs spatiaux usuels
sont décrites dans les paragraphes suivants.

5.3.1 Approximation de la dérivée première

Pour l’opérateur L = d/dx, on tire de l’Eq. (5.9a) :

µ ¶
ui+1 − ui ∂u
= + τi , (5.10a)
∆x ∂x i

avec τi l’erreur d’approximation ou erreur de troncature, déjà rencontrée dans le chapitre 3, qui s’écrit :
µ ¶
∆x ∂2u
τi = + O(∆x2 ). (5.10b)
2 ∂x2 i

D’une manière similaire, l’équation (5.9b) donne

µ ¶ µ ¶
ui − ui−1 ∂u ∆x ∂ 2 u
= − + O(∆x2 ). (5.11)
∆x ∂x i 2 ∂x2 i

Pour ces deux formules aux différences, le terme dominant de l’erreur de troncature est proportionnel à ∆x :
on dit que ces formules sont au premier ordre. L’erreur de troncature mesure la précison de l’approximation
et détermine le taux de décroissance de l’erreur lorsque l’espacement des points ∆x diminue. Une présentation
plus formelle des notions de troncature et de convergence est effectuée dans le paragraphe 5.4.
Il est possible d’obtenir des approximations de la derivée première d’ordre supérieur à 1. En effet, la combi-
naison des Eqs. (5.9a) et (5.9b) donne
µ ¶ µ ¶
ui+1 − ui−1 ∂u ∆x2 ∂ 3 u
= + + O(∆x3 ). (5.12)
2∆x ∂x i 6 ∂x3 i

Cette approximation est au second-ordre, et elle est formellement plus précise que (5.10) ou (5.11).
Les formules de différences finies sont maintenant obtenues en négligeant l’erreur de troncature. On obtient
alors :
Définition 5.1 (Approximations usuelles de la dérivée première).
36 5 Concepts de base des différences finies

Fig. 5.2. Approximations usuelles de la dérivée première au nœud xi .

• Schéma décentré amont (“upwind”, ou “backward”), en O(∆x) :

ui − ui−1
δx− ui ≡ , (5.13)
∆x
• Schéma décentré aval (“forward”), en O(∆x) :
ui+1 − ui
δx+ ui ≡ , (5.14)
∆x
• Schéma centré (“central”), en O(∆x2 ) :
ui+1 − ui−1
δx0 ui ≡ . (5.15)
2∆x
Une interprétation géométrique de ces formules est présentée sur la figure 5.2. Il est évident d’après cette figure
que certaines formules sont meilleures que d’autres. Les propriétés de chacune des ces approximations seront
mises en évidence dans le chapitre 9.

5.3.2 Approximation de la dérivée seconde

L’approximation usuelle de la dérivée seconde est la formule centrée à trois points, déjà présentée dans le
paragraphe 3.1 : µ 2 ¶ µ ¶
ui+1 − 2 ui + ui−1 ∂ u ∆x2 ∂ 4 u
= + + O(∆x3 ). (5.16)
∆x2 ∂x2 i 12 ∂x4 i
Définition 5.2 (Approximation de la dérivée seconde). L’approximation usuelle est la formule centrée à
trois points :
0 ui+1 − 2 ui + ui−1
δxx ui ≡ . (5.17)
∆x2
dont l’erreur de troncature est en O(∆x2 ).

Remarque 5.3. Il est possible de construire des schémas pour les dérivées secondes (ou d’ordre plus élevé) à
partir de formules pour les dérivées d’ordre inférieur. Par exemple,

0 1 +
δxx ui = δx+ (δx− ui ) = δx− (δx+ ui ) = (δ − δx− )ui ,
∆x x

par contre

ui+2 − 2 ui + ui−2
δx0 δx0 ui = ¡ ¢2
0
6= δxx ui .
2∆x
 5.4 Analyse des schémas aux différences 37

5.3.3 Technique générale

Il existe des manières plus systématiques pour obtenir des formules aux différences. Par exemple, pour obtenir
une approximation la dérivée première de u(x) au point xi qui utilise les valeurs aux nœuds i − 1, i et i + 1, on
écrit µ ¶
∂u
a ui+1 + b ui + c ui−1 = + τi , (5.19)
∂x i
où a, b et c sont les coefficients à déterminer et τi est le terme de troncature indiquant la précision de la formule.
La substitution de développements de Taylor dans (5.19) conduit à :
³ µ ¶ µ ¶ µ ¶
¡ ¢ 1 ´ ∂u ¡ ¢ ∆x2 ∂ 2 u ¡ ¢ ∆x3 ∂ 3 u
τi = a + b + c ui + a − c − ∆x + a+c + a − c + O(∆x4 ), (5.20)
∆x ∂x i 2 ∂x2 i 6 ∂x3 i

en ordonnant le second-membre en puissances croissantes de ∆x. Pour obtenir la formule la plus précise, il faut
annuler autant de termes que possible dans l’erreur de troncature. Puisque nous disposons de trois coefficients,
il est possible de mettre les trois premiers termes de (5.20) à zéro, ce qui conduit à :
1
a + b + c = 0, a−c= , a + c = 0.
∆x
Ce système linéaire possède l’unique solution :
1 1
a= , b = 0, c=− ,
2∆x 2∆x
et la substitution de ces valeurs dans (5.20) conduit à
µ ¶
∆x2 ∂ 3 u
τi = + O(∆x4 ).
6 ∂x3 i

On aura bien sûr reconnue l’approximation centreé δx0 (Eq. (5.12)).

Cett technique permet de construire des formules de différences finies dont l’ordre de précision est arbitraire.
Elle est très utilisée pour construire des formules d’ordre élevé3 , mais aussi pour construire des formules sur
maillages non-uniformes.
Exercice 5.4. Utilisez la technique générale pour retrouver les formules et erreurs de troncature des sections
5.3.1 et 5.3.2.

Exercice 5.5. Déterminez la formule décentrée du second ordre qui exprime (∂u/∂x)i en fonction de ui , ui+1 ,
et ui+2 . Montrez que l’erreur de troncature s’écrit :
µ ¶
∆x2 ∂ 3 u
τi = − .
3 ∂x3 i

Exercice 5.6. Montrez que sur un maillage non-uniforme, la technique générale permet de construire la formule
centrée suivante pour la dérivée première :
µ ¶n µ ¶n
∂u (∆xi )2 n (∆xi+1 )2 − (∆xi )2 n (∆xi+1 )2 n ∆xi+1 ∆xi ∂ 3 u
= ui+1 + ui − ui−1 − ,
∂x i D D D 6 ∂x3 i

avec D = ∆xi+1 ∆xi (∆xi+1 + ∆xi ) et ∆xi défini en (5.2). Retrouvez l’expression de la dérivée première centrée
dans le cas particulier d’un maillage uniforme.

Exercice 5.7. Montrez que la formule donnée en (7.21) est bien du quatrième ordre.

5.4 Analyse des schémas aux différences

Considérons que, pour résoudre le problème aux limites générique

A u(t, x) = f (t, x), (5.21)

les paragraphes précédents ont permis de construire un schéma aux différences que nous écrivons en toute
généralité sous sa forme nodale :
Ah uni = fin . (5.22)
3
On qualifie traditionnellement d’“ordre élevé” toute méthode dont l’ordre de précision est strictement supérieur à 2.
38 5 Concepts de base des différences finies

Ce schéma peut s’écrire sous la forme du système linéaire (cf. Chap. 3) :

Ah U n = F n , (5.23)

avec U n = (uni , i = 1, . . . , N ) et F n = (fin , i = 1, . . . , N ), et Ah la représentation matricielle du schéma Ah .

La solution discrète uni est une approximation de la solution exacte u(tn , xi ) aux nœuds du maillage, i.e. :

uni ∼
= u(tn , xi ),

et on a bien sûr pour le second membre fni = f (tn , xi ). La question centrale de l’analyse du schéma (5.22) est
la suivante : comment se convaincre qu’il nous donnera une “bonne” solution uni , qui convergera vers la solution
exacte lorsque les pas de discrétisation ∆t, ∆x tendent vers zéro, si on ne connait pas la solution exacte ? Nous
verrons dans ce paragraphe que le théorème de Lax permet d’assurer la convergence du schéma vers la bonne
solution dès que :
• le schéma est “correct” ou consistant, i.e. il représente la bonne EDP, et
• il est stable, i.e. une petite perturbation, ou erreur, dans les données initiales ne va pas s’amplifier de façon
incontrolée. Ceci est surtout pertinent pour les EDP instationnaires.
Mais avant tout, ce paragraphe va introduire une certain nombre de concepts fondamentaux tels la consis-
tance et la stabilité, qui permettront d’évaluer la qualité d’un schéma numérique avant de le programmer sur
un ordinateur, et d’éviter de mauvaises surprises quand on se rend compte qu’il est inefficace pour le problème
qui nous intéresse !

5.4.1 Consistance d’un schéma numérique

Puisque Ah représente une approximation tronquée de A, la solution discrète uni ne peut être strictement égale
à la solution exacte u(tn , xi ). Comme dans l’exemple de la poutre (Eq. 3.8), on définit alors

Ah u(tn , xi ) = fin + τin , (5.24)

où τin est l’erreur de troncature en chaque nœud (tn , xi ) du maillage. Puisque fin = A u(tn , xi ) d’après (5.21),
on obtient :
Ah u(tn , xi ) = A u(tn , xi ) + τin .
On définit alors :
Définition 5.8 (Erreur de troncature d’un schéma). L’erreur de troncature locale τin , est la différence
entre le schéma numérique et l’équation à résoudre, soit :
³ ´
τin = Ah u(tn , xi ) − fin − A u(tn , xi , ) − f (tn , xi ) . (5.25)
| {z } | {z }
Schéma numérique Equation en (tn , xi )

Cette erreur est définie uniquement aux nœuds de discrétisation (tn , xi ).

L’erreur de troncature va nous servir à évaluer la capacité d’un schéma à représenter correctement l’équation à
résoudre. En effet, il est naturel de s’attendre à ce que pour un schéma “correct”, ou consistant, l’erreur τin soit
petite. La consistance d’un schéma, qui est la propriété de représenter “correctement” une EDP, est définie de
la manière suivante :
Définition 5.9 (Consistance d’un schéma).
Un schéma est dit consistant si l’erreur de troncature tends vers zéro lorsque tous les pas de discrétisation
tendent vers zéro, i.e.
τin −→ 0 pour ∆t, ∆x −→ 0. (5.26)
De plus, on dira qu’un schéma est consistant d’ordre p en temps et d’ordre q en espace si l’erreur de troncature
vérifie
τin = O(∆tp ) + O(∆xq ). (5.27)
 5.4 Analyse des schémas aux différences 39

5.4.2 Analyse des schémas stationnaires

Reprenons l’exemple du fléchissement de poutre du chapitre 3. La discrétisation de l’équation (3.1a) avec c(x) = 0
par le schéma centré à trois points s’écrit :
0
−δxx ui = fi . (5.28)
L’erreur de troncature du schéma est :
³ ´ ³ d2 u ´
0
τi = −δxx u(xi ) − fi − − 2 (xi ) − f (xi ) (5.29)
| {z } | dx {z }
Schéma numérique Equation en xi
2 4
∆x ∂ u
=− (xi ) + O(∆x3 ), (5.30)
12 ∂x4
en utilisant (5.16). L’erreur de troncature est donc de la forme τi = O(∆x2 ), ce qui montre que le schéma centré
est consistant au second ordre.
Cette propriété de consistance va nous servir à montrer que le schéma est convergent, c.à.d. :

|uni − u(tn , xi )| −→ 0 lorsque ∆x −→ 0,

en tout point xi du maillage. Pour cela, remarquons que :

Ah ui − Ah u(xi ) = fi − Ah u(xi ),
= Ah u(xi ) − τi − Ah u(xi ) d’après (5.24),
= −τi .

On en déduit l’inégalité vectorielle :

(−1)
||E||2 ≤ |||Ah |||2 ||τ ||2 , (5.31)
¡ ¢
où E = ui − u(xi ), i = 1, . . . , N désigne l’erreur aux nœuds du maillage. Il suffit maintenant de montrer que
(−1)
la norme de Ah reste bornée pour ∆x petit, soit :
(−1)
∃C > 0 tel que |||Ah |||2 ≤ C lorsque ∆x −→ 0. (5.32)

pour que (5.31) devienne pour ∆x petit :

||E||2 ≤ C ||τ ||2 . (5.33)
Si le schéma est consistant, alors ||τ || −→ 0 et la solution numérique convergera vers la solution exacte. La
propriété (5.32) est la stabilité d’un schéma numérique ; on aura ainsi demontré pour cet exemple simple que :

consistance + stabilité =⇒ convergence,

qui est l’essence du théorème de Lax, qui est un des théorèmes fondamentaux de l’analyse des schémas
numériques.
Vérifier la consistance d’un schéma est la partie aisée de l’analyse de convergence. La stabilité est la partie
difficile. Pour cet exemple simple, montrons-la dans le cas des conditions de Dirichlet (3.1b). D’après les résultats
de l’exercice 3.2, on a :
(−1) 1 ∆x2
|||Ah |||2 = = .
ρ(Ah ) min λk
k=1,...,N

où les λk sont donnés par (3.12a). On obtient alors

(−1) ∆x2 ∼ 1
|||Ah |||2 = = 2,
1 − cos πh π
lorsque ∆x → 0. L’estimation d’erreur (5.31) devient donc :

||τ ||
||E||2 ≤ = O(∆x2 ). (5.34)
π2
Le schéma est donc convergent au second ordre pour des conditions de Dirichlet. Le cas des conditions de
Neumann sera abordé dans la section 7.2.
40 5 Concepts de base des différences finies

5.4.3 Analyse des schémas d’évolution

Intéressons-nous maintenant aux problèmes instationnaires, en considérant l’équation (parabolique) de la chaleur

∂u ∂2u
−ν = f (t, x). ν ≥ 0. (5.35)
∂t ∂x2
La première étape de l’approximation par différences finies de cette EDP consiste à discrétiser uniquement la
0
partie spatiale de l’EDP. L’utilisation du schéma centré δxx conduit à

dui 0
(t) − ν δxx ui (t) = fi (t). (5.36)
dt
Cette étape intermédiaire, qui est habituellement appelée semi-discrétisation, vient de transformer l’EDP
(5.35) en un système d’équations différentielles linéaires. Il est alors possible d’utiliser les schémas classiques
d’intégration des équations différentielles ordinaires.

Un exemple de schémas d’évolution : les schémas d’Euler

Dans la suite, nous allons appliquer les méthodes d’Euler (5.5a) et (5.5b), qui sont certainement parmi les plus
simples que l’on puisse considérer. Ces schémas peuvent aussi être utilisés pour l’équation d’advection-diffusion
générale (4.11).
Exemple 5.10 (Schéma d’Euler progressif ).

Lorsqu’on applique le schéma temporel (5.5a) à la sémi-discretisation (5.36), on obtient simplement :

un+1 − uni 0
i
− ν δxx uni = fin . (5.37)
∆t
Pour établir l’erreur de troncature de ce schéma de façon pratique, il faut porter le développement de Taylor
(5.9c) et la formule (5.16) dans (5.37) pour obtenir
µ ¶n µ 2 ¶n
uin+1 − uni 0 n ∂u ∂ u
− ν δxx ui = −ν + τin ,
∆t ∂t i ∂x2 i
avec
µ ¶n µ ¶n
∆t ∂2u ∆x2 ∂4u
τin = −ν + O(∆t2 , ∆x3 ),
2 ∂t2 i 12 ∂x4 i
2
= O(∆t, ∆x ).

Le schéma d’Euler progressif (5.37) est donc consistant au premier ordre en temps, et au second ordre en espace.
Afin d’utiliser ce schéma dans un programme informatique, il est intéressant de le mettre sous sa forme
“itérative” µ ¶
n+1 ν∆t n 2ν∆t ν∆t n
ui = u + 1− uni + u + ∆t fin , (5.38)
∆x2 i+1 ∆x2 ∆x2 i−1
soit, matriciellement,

U n+1 = B U n + G, (5.39a)

avec

B = I + ν∆t Hh , G = ∆t F n . (5.39b)

De plus, I désigne la matrice identité et Hh est définie par

 
1 −2 1
1  .. .. .. 
Hh =  . . . . (5.40)
∆x2
1 −2 −1

Le calcul de U n+1 , solution du problème au temps tn+1 , est obtenu en fonction de U n . Ce calcul par récurrence,
ou avancement en temps, de la solution aux instants successifs t0 , t1 , . . . , tn ne nécessite pas l’inversion de
systèmes linéaires. On dit que le schéma d’Euler progressif est explicite en temps.
 5.4 Analyse des schémas aux différences 41

Remarque 5.11. La matrice tridiagonale Hh est la discrétisation de l’opérateur d2 /dx2 par le schéma centré du
second ordre. On a déjà rencontré cette matrice dans l’exemple de la poutre (Eq. (3.8), avec ci = 0 et si on
“oubliait” les conditions aux limites).

Définition 5.12 (Schémas multiniveaux).

Plus généralement, un schéma aux différences finies pour une équation d’évolution qui s’écrit sous la forme
itérative :
B1 U n+1 = B2 U n + · · · + Bl U n+2−l + G, (5.41)
est un schéma à l niveaux de temps. Si B1 = I, le schéma est dit explicite. Sinon, le schéma est implicite, et le
calcul de U n+1 nécessite l’inversion d’un système linéaire.

Exemple 5.13 (Schéma d’Euler rétrograde).

Le schéma temporel (5.5b) appliqué à la semi-discrétisation (5.36) conduit à :

un+1 − uni
i 0
− ν δxx un+1
i = fin+1 . (5.42)
∆t
L’étude de sa consistance est laissée en exercice (comparer avec le schéma (5.37)). Il s’écrit matriciellement

(I − ν∆t Hh ) U n+1 = U n + ∆t F n+1 . (5.43)

Il est donc de la forme (5.41) avec

B1 = I − ν∆t Hh , B2 = I, G = ∆t F n+1 .

C’est donc un schéma implicite à deux niveaux de temps.

Les schémas d’Euler progressif et rétrograde sont tous deux des schémas à deux niveaux de temps dont l’ordre
de troncature est O(∆t, ∆x2 ). Le premier, de type explicite, est moins coûteux en temps de calcul puisqu’il ne
nécessite pas la résolution de systèmes linéaires. Il paraît donc plus avantageux sur le plan algorithmique. Il y
a néanmoins un prix à payer pour pouvoir l’utiliser : la stabilité conditionnelle de la solution calculée.

Stabilité d’un schéma d’évolution

Pour un problème d’évolution, la stabilité est la propriété du schéma à fournir une solution numérique qui reste
bornée, sur un maillage donné, quand t → ∞ (∆t fixé) ou quand ∆t → 0 (à t fixé).
Tous les schémas d’évolution, même les schémas multiniveaux de la forme (5.41) (cf. TD), peuvent se réécrire
sous la forme itérative :
U n+1 = B U n + G, (5.44)
où B est la matrice d’itération du schéma. Par exemple, la forme itérative du schéma d’Euler progressif (5.37)
est donnée par (5.39).
Lorsqu’on enlève le terme de force (G ≡ 0), le schéma (5.44) prend la forme équivalente :

U n+1 = B(n+1) U 0 ,

dont on déduit les inégalités

kU n+1 k ≤ |||B (n+1) ||| kU 0 k ≤ |||B|||(n+1) kU 0 k. (5.45)

Pour que la solution reste bornée à tout temps et que les erreurs initiales ne s’amplifient pas, on énonce la
définition :
Définition 5.14 (Stabilité d’un schéma). Le schéma (5.44) est stable si pour ∆t et ∆x suffisamment petits
sa matrice d’itération vérifie :
|||B||| ≤ 1. (5.46)
Le schéma est inconditionnellement stable si (5.46) est vraie pour toute valeur de ∆t. Le schéma est condition-
nellement stable si (5.46) est vraie pour certaines valeurs de ∆t seulement. Dans ce dernier cas, une condition
de stabilité, analogue à celle qu’on trouve pour les EDO, restreint les valeurs que ∆t peut prendre.
Dans la pratique, nous utiliserons l’étude de stabilité la plus usuelle, qui est celle de Von Neumann. La présen-
tation de cette méthode, et son application concrète aux schémas d’Euler du paragraphe 5.4.3, sera effectuée
dans le chapitre 6 consacré aux équations paraboliques.
42 5 Concepts de base des différences finies

5.4.4 Convergence d’un schéma aux différences

Dans le paragraphe 5.4.2, nous avons observé que pour un exemple stationnaire simple, on avait

consistance + stabilité d’un schéma =⇒ convergence.

Ces notions vont être maintenant formalisées et étendues au cas d’un problème d’évolution.

Le point de vue formel

La convergence d’un schéma est la caractéristique principale de la qualité de l’approximation discrète uni par
rapport à la solution exacte u(tn , xi ). Plus précisement :
Définition 5.15 (Convergence).
Soit Uen = (u(tn , xi ), i = 1, . . . , N ) le vecteur représentant la solution exacte aux nœuds du maillage. Un schéma
est dit convergent si
∀ n, kU n − Uen k −→ 0 pour ∆t, ∆x → 0.
Un schéma est dit convergent d’ordre p en temps et d’ordre q en espace si

∀ n, kU n − Uen k = O(∆tp ) + O(∆xq ).

La convergence d’un schéma est, en général, très difficile à demontrer mathématiquement. Néanmoins, on dispose
du théorème suivant :
Théorème 5.16 (Lax).
Pour un problème aux limites bien posé et un schéma aux différences qui possède la propriété de consistance, la
stabilité est la seule condition nécessaire et suffisante pour assurer sa convergence.

Preuve. De la même manière qu’en (5.24), la solution exacte satisfait l’équation

Uen+1 = B Uen + Gn + ∆t τ n , (5.47)

où τ n désigne l’erreur de troncature du schéma (5.44). Soustraire (5.47) à (5.44) donne l’équation vérifiée par
l’erreur E n = U n − Uen :
E n+1 = B E n − ∆t τ n .
Par induction, on obtient :
n
X
E n = B (n) E 0 − ∆t B (n−k) τ k−1 ,
k=1

dont on déduit,
n
X
||E n || ≤ |||B|||(n) ||E 0 || + ∆t |||B|||(n−k) ||τ k−1 ||,
k=1
n o
(n) 0
≤ |||B||| ||E || + tn max |||B|||(n−k) ||τ k−1 || ,
k=1,...,n

avec tn = n∆t. Si la methode est consistante (i.e. ||τ k−1 || → 0 pour k = 1, . . . , n) et si elle est stable (|||B|||(n) →
0 pour n → ∞), alors la convergence est assuré.

L’importance de ce théorème est grande, car il est bien plus aisé de vérifier la consistance et la stabilité d’un
schéma que la convergence de son approximation. Néanmoins, ce théorème a ses limitations. Les EDP de
la mécanique sont très souvent non-linéaires et possèdent des conditions aux limites complexes, et dans ces
conditions il est difficile d’appliquer le théorème de Lax rigoureusement. On doit l’interpréter comme si les
propriétés de stabilité et de convergence sont des modes de sélection des schémas que l’on construit.

Un exemple mathématique

Pour l’équation (3.1) avec c(x) = 0, discrétisée par le schéma centré (3.9), voici le type d’estimation d’erreur
que l’on peut trouver dans la littérature mathématique (cf. [11]) :
Proposition 5.17. Si on suppose que la solution u(x) est dans C 4 (]0, 1[), on dispose alors de la majoration
d’erreur suivante :
∆x2
ku(xi ) − ui kl∞ ≤ max |u(4) (x)|. (5.48)
96 x∈[0,1]
 5.4 Analyse des schémas aux différences 43

Un résultat similaire a été obtenu dans la section 5.4.2.

La proposition 5.48 établit la convergence au second ordre de la solution sous l’hypothèse que la solution
soit suffisament régulière4 . L’estimation d’erreur (5.48) introduit la norme discrète k · kl∞ , définie sur la grille
de calcul Ωh = {xi , i = 1, . . . , N } comme :

kφkl∞ = max |φ(xi )|. (5.49a)

xi ∈Ωh

Une autre norme utilisée en pratique est la norme discrète l2 , définie par :
à N
!1/2
1 X
kφkl2 = |φ(xi )|2 . (5.49b)
N i=1

On remarquera que ces normes sont des approximations ad hoc des normes de Sobolev continues :

kφkL∞ (Ω) = max |φ(x)|, (5.50a)

x∈Ω
µZ ¶1/2
kφkL2 (Ω) = |φ(x)|2 dx . (5.50b)
Ω

Un exemple concret

Comme illustration de ce paragraphe, on considère la résolution du problème d’Helmholtz dans Ω =] − 1, 1[

−νu00 (x) + u(x) = 1, ν = 10−2 , (5.51a)

u(−1) = 1, u(1) = 0, (5.51b)

dont la solution est C ∞ , avec les méthodes numériques suivantes :

1. Différences finies centrées du second ordre,
2. Différences finies compactes du sixième ordre (cf. e.g. [7, 9]),
3. Méthode de collocation-Tchebychev [6].

0
10

−5
10

−10
10

−15
10
10 20 30 50
2
Fig. 5.3. Erreur l en fonction du nombre de points N (coordonnées logarithmiques) pour le problème (5.51) résolu
par différentes méthodes numériques; 2 : différences finies centrées; 3 : différences finies compactes; ° : collocation-
Tchebychev.

La figure 5.3 montre que, pour les schémas aux différences finies 1. et 2., une erreur du second et sixième ordre
est asymptotiquement obtenue en accord avec les estimations d’erreur théoriques. La méthode de collocation-
Tchebychev démontre pour ce problème très régulier une convergence exponentielle : ce type de convergence est
typique des méthodes spectrales (e.g. : séries de Fourier).

4
Ici u ∈ C 4 (]0, 1[), i.e. la dérivée quatrième de u(x) existe et est continue dans ]0,1[.
6 Les problèmes paraboliques

6.1 Généralités
Les problèmes d’évolution temporelle qui sont de nature dissipative sont gouvernés par des équations aux dérivées
partielles paraboliques.
Le prototype de ce type de problèmes est l’équation de la chaleur pour la température T (t, x)

∂T
ρcp − κ∇ · κ∇T = q, ∀ t ∈ [0, tmax ], ∀ x ∈ Ω ⊂ R2 , (6.1)
∂t
où cp désigne la chaleur spécifique à pression constante, κ > 0 la conductivité thermique et q une source de
chaleur. On adjoint à cette équation une condition initiale de la forme

T (0, x) = T0 (x), ∀ x ∈ Ω, (6.2)

et on considère généralement deux types de conditions sur la paroi ∂Ω.

• Condition de Dirichlet :
∀ t, T (t, x) = Tw (t, x), ∀ x ∈ ∂Ω, (6.3)
où Tw est une température prescrite à la paroi.
• Condition de Neumann :
∂T
∀ t, −h (t, x) = qw , ∀ x ∈ ∂Ω, (6.4)
∂n
où h est le coefficient de transfert thermique et qw le flux de chaleur à travers la paroi. Cette condition aux
limites modélise la quantité de chaleur qui est échangée aux frontières du domaine. Lorsque qw = 0, cette
condition aux limites est appelée adiabatique, et énonce qu’aucun flux de chaleur ne traverse la paroi.
Sur le plan mathématique, ce problème aux limites est bien posé pour tout κ > 0.

Fig. 6.1. Distribution instantanée de la température dans une barre de métal.

Exemple 6.1. Considérons le problème représenté sur la figure 6.1, qui concerne la distribution instantanée de
la chaleur dans une barre de métal qui a été portée à la température T0 sous l’effet d’une source extérieure,
qui a ensuite été désactivée au temps t = 0. La température au temps initial est donc T (x, 0) = T0 (x) et, pour
t > 0, les extrémités x = 0 et x = L de la barre sont maintenues à la température constante Tw . Ce problème
est gouverné par l’EDP
∂T ∂2T
− κ 2 = 0. (6.5)
∂t ∂x
46 6 Les problèmes paraboliques

L’évolution temporelle de la solution est représentée sur la figure 6.1 (droite). Comme il a été observé dans le
paragraphe 4.2, on observe une décroissance exponentielle de la condition initiale jusqu’à obtenir une distribution
uniforme T (x) = Tw dans la barre lorsque t → ∞. Cette solution stationnaire est, en fait, solution de l’EDP

d2 T
−κ = 0, (6.6)
dx2
qui est de nature elliptique. Cette dernière équation est obtenue en posant ∂T /∂t = 0 dans (6.5).
Cette exemple illustre le lien entre les EDP paraboliques et elliptiques. Il est donc naturel d’utiliser les mêmes
opérateurs de différences finies pour la partie spatiale de ces problèmes, en particulier les schéma centrés δx0 et δxx
0

du paragraphe 5.3. Le présent chapitre porte essentiellement sur l’approximation temporelle dans les équations
paraboliques. Les conditions aux limites s’implémentent de la même manière que pour les équations elliptiques,
et sont traitées dans le chapitre 7. Les extensions multidimensionnelles seront traitées dans le chapitre 8.

6.2 Semi-discrétisation et lien avec les équations différentielles

Considérons, pour simplifier, l’équation de la chaleur 1D
∂u ∂2u
− ν 2 = f (x, t), ν > 0. (6.7)
∂t ∂x
Si on discrétise la partie spatiale de l’opérateur uniquement par un schéma aux différences finies quelconque
Ah , on obtient matriciellement
dU
(t) − νAh U (t) = F (t), (6.8)
dt
avec U (t) = (ui (t), i = 1, . . . , N ).
On vient de réduire l’EDP (6.7) à un système d’équations aux dérivées ordinaires pour le vecteur U , qui
ne dépend que du temps. Il est alors possible d’utiliser les schémas d’intégration classique pour EDO (Euler,
Runge-Kutta, . . . cf. [10, 9]), pour construire un schéma de discrétisation complet. Dans la section 5.4.3, le
principe de semi-discrétisation a été appliquée à la construction des schémas d’Euler progressif et rétrograde
pour l’équation de la chaleur.
Cette méthode de construction, très générale, est indépendante du schéma de discrétisation spatial A, et est
utilisable pour d’autres équations d’évolution, comme par exemple les équations hyperboliques. Elle applicable
pour d’autre discrétisations que les différences finies (e.g. volumes finis, élements finis, . . . ). Une fois le schéma
construit, il reste néanmoins essentiel de s’assurer de sa stabilité.

6.3 L’étude de stabilité de Von Neumann

D’après la section 5.4.3, la stabilité d’un schéma numérique écrit sous sa forme itérative :

U n+1 = B(n+1) U 0 ,

est la propriété à donner une solution U n+1 qui reste bornée lorsque n → ∞. La définition 5.14 enonce comme
condition nécessaire de stabilité de la solution numérique :

|||B||| ≤ 1. (6.9)

L’étude de stabilite d’un schéma revient à determiner les valeurs de ∆t pour que (6.9) soit vérifiée et que le
calcul n’explose pas. Ce critêre depend de la norme matricielle utilisée. Le plus souvent, on utilise la norme
||| · |||2 qui pour une matrice symétrique s’écrit (cf. Sect. 1.3.1) :

|||B|||2 = ρ (B) , (6.10)

où ρ ≥ 0 est le rayon spectral, défini par :

ρ (B) = max { |λi |, i = 1, . . . , N } , (6.11)

tel que les λi sont les valeurs propres de B.

Cette étude peut néanmoins se révéler fastidieuse pour un maillage non-uniforme et des conditions aux
limites compliquées. Dans la pratique, l’analyse de stabilité est effectuée dans le cas d’un maillage uniforme, en
“oubliant” les conditions aux limites : en toute rigueur, cela équivaut à considérer des conditions de périodicité.
Une analyse de Fourier est alors possible et les conditions de stabilité du schéma sont facilement obtenues : c’est
la méthode de Von Neumann.
 6.3 L’étude de stabilité de Von Neumann 47

6.3.1 Présentation formelle

Pour un problème non forcé (f ≡ 0), nous allons considérer sans perte de généralité le schéma à deux niveaux
sous forme itérative :
B1 U n+1 = B2 U n . (6.12)
Cherchons une solution à cette équation sous la forme du mode de Fourier

bn e i kxi ,
uni = u (6.13)

avec i2 = −1 et xi = i∆x. L’équation (6.12) devient donc

bn+1 = G2 u
G1 u bn , (6.14a)

qu’on réécrit sous la forme

bn+1 = G u
u bn , (6.14b)

où G = (G1 )−1 G2 est appelée matrice d’amplification. La définition suivante énonce simplement que u
bn+1 reste
bornée dans la limite n → ∞ lorsque les valeurs propres de G sont de module inférieur à 1 :
Définition 6.2 (Stabilité au sens de Von Neumann).
Le schéma (6.12) est stable si
ρ(G) ≤ 1, ∀k tel que − π ≤ k∆x ≤ π, (6.15)
où ρ(G) est le rayon spectral de la matrice d’amplification.

Remarque 6.3. Dans le cas particulier d’un schéma 1D possédant 2 niveaux de temps, la matrice G se réduit au
nombre g = g2 /g1 , appelé facteur d’amplification, qui peut prendre des valeurs complexes. La condition de Von
Neumann (6.15) se réduit alors à

|g| ≤ 1, ∀k tel que − π ≤ k∆x ≤ π. (6.16)

Remarque 6.4. Si le schéma à analyser possède plus de 2 niveaux de temps, il est toujours possible de se ramener
à une expression du type (6.12) en introduisant des variables auxiliaires (cf. TD).

6.3.2 Application aux schémas d’Euler

Dans cette partie, on va utiliser la méthode de Von Neumann pour étudier la stabilité des schémas d’Euler du
paragraphe 5.4.3.
Exemple 6.5 (Suite de l’exemple 5.10). On a vu précédemment que le schéma d’Euler explicite (5.37) pour
l’équation de la chaleur (avec f ≡ 0) peut s’écrire sous la forme itérative (5.38), soit :

un+1
i = runi+1 + (1 − 2r)uni + runi−1 , (6.17)

avec
ν∆t
r= . (6.18)
∆x2
En utilisant (6.13), l’Eq. (6.17) s’écrit
bn+1 = g u
u bn , (6.19)
avec le facteur d’amplification
g = (1 − 2r) + 2r cos k∆x. (6.20)
Il est souvent pratique d’utiliser la variable auxiliaire
k∆x 1 − cos k∆x
ξ = sin2 = , (6.21a)
2 2

qui vérifie ξ ≥ 0. Notez qu’une autre relation utile pour la suite est

sin2 k∆x = 4ξ(1 − ξ). (6.21b)

Le facteur d’amplification g devient alors

g = 1 − 4rξ. (6.22)
La condition de Von Neuman (6.16) s’écrit alors
48 6 Les problèmes paraboliques
2
|g| ≤ 1 (⇐⇒ |g| ≤ 1) pour ξ ∈ [0, 1]. (6.23)

Apres quelques calculs1 , cette condition s’écrit

rξ(2rξ − 1) ≤ 0, ∀ξ ∈ [0, 1], (6.24a)

qui devient, puisque ξ ≥ 0 et r ≥ 0

2rξ ≤ 1, ∀ξ ∈ [0, 1]. (6.24b)

Puisque le membre de gauche de l’inégalité est une fonction convexe de ξ, la condition (6.24b) est vraie dès que
l’inégalité est vérifiée en ξ = 0 et ξ = 1, soit

1
r≤ . (6.24c)
2
On en déduit finalement la condition de stabilité du schéma d’Euler progressif :
1
∆t ≤ ∆x2 . (6.25)
2ν
Ce schéma est donc conditionellement stable.

Exemple 6.6 (Suite de l’exemple 5.13).

Le schéma d’Euler implicite (5.42) sans terme force peut s’écrire sous la forme itérative :

−r un+1 n+1
i+1 + (1 + 2r) ui − run+1 n
i−1 = ui , (6.26)

Son facteur d’amplification est

1
g= . (6.27)
1 + 4rξ
Comme le module de g est toujours inférieur à l’unité, ce schéma est stable sans condition.

6.3.3 Bilan des courses

• Schéma explicite
− Avantage : Facile à mettre en œuvre.
− Désavantage : Limitation de la forme ∆t ≤ f (∆x). En pratique, on peut être amené à utiliser un pas de
temps très petit pour maintenir la stabilité du calcul, ce qui entraîne un très long temps de calcul2 pour
aller de t = 0 a t = T .
• Schema implicite
− Avantage : Très bonne stabilité : autorise des ∆t beaucoup plus grand qu’un schéma explicite.
− Désavantage : Entraîne la reśolution d’un système linéaire (e.g. Eq. (5.43)) à chaque niveaux de temps. En
conséquence, il faut utiliser des algorithmes de résolution de systèmes linéaires très performants (le plus
souvent des algorithmes itératifs) pour minimiser le coût CPU par cycle de temps.

1
Conseil : élever (6.22) au carré.
2
Aussi appelé temps CPU.
7 Les problèmes elliptiques

7.1 Généralités
Les problèmes de mécanique qui décrivent des phénomènes physiques à l’état d’équilibre, où le temps n’intervient
pas, sont souvent gouvernés par des EDP elliptiques. Parmi ces problèmes, citons la distribution de température
dans une barre de métal (cf. Fig. 7.1), la distribution des contraintes dans un solide soumis à une déformation,
ainsi que de nombreux problèmes d’écoulements de fluides. Le prototype des problèmes elliptiques est l’équation
de Laplace, qui s’écrit en 2D
∂2u ∂2u
+ 2 = 0,
∂x2 ∂y
et qui décrit aussi bien l’écoulement irrotationel d’un fluide incompressible, que le transfert de chaleur station-
naire.
Exemple 7.1 (Les équations de l’élasticité linéaire).

∇ · [σ] + ρf = 0, (7.1a)
[σ] = λ tr [²] + 2µ [²], (7.1b)

où [σ] est le tenseur des contraintes et [²] le tenseur des déformations qui a pour composantes :
µ ¶
1 ∂ui ∂uj
²ij = + , (7.1c)
2 ∂xj ∂xi

forment un système elliptique de 3 équations pour le déplacement u. En 1D, on obtient :

µ ¶
∂ ∂u
(λ + 2µ) = −ρf. (7.2)
∂x ∂x

Les équations de l’élasticité linéaire sont des problèmes elliptiques du second ordre puisque des dérivées secondes
au plus interviennent dans les équations.

Exemple 7.2. Un petit peu de vocabulaire : considérons l’équation aux dérivées partielles générique

−∇2 u + λu = f, avec λ ≥ 0. (7.3)

• Si λ ≡ 0 et f ≡ 0, elle est appelée équation de Laplace,

• Si λ ≡ 0 et f 6≡ 0 : c’est l’équation de Poisson,
• Si λ > 0 : c’est l’équation de Helmholtz.

Fig. 7.1. Distribution de température dans une barre de métal à l’équilibre.

50 7 Les problèmes elliptiques

Un exemple très simple de phénomène d’équilibre est le problème de conduction de chaleur stationnaire
dans une barre de métal représenté sur la figure 7.1, dont les extrémités x = 0 et x = L sont portées à des
températures constantes T0 et TL , telles que TL > T0 . Ce problème est gouverné par l’équation de Laplace 1D
∂2T
κ = 0.
∂x2
La solution de ce problème est, bien sûr, la droite représentée sur la figure 7.1.
Une caractéristique importante des équations elliptiques est qu’une perturbation de la solution en un point
du domaine, par exemple un changement de la température dû à l’apparition d’une source de chaleur ponctuelle,
entraîne un changement de la solution dans le domaine tout entier. C’est particulièrement évident si on perturbe
les conditions aux limites du problème de la figure 7.1.
L’impact de ce phénomène sur l’imposition des C. L. dans un schéma numérique sera illustré dans la section
7.2. Une autre conséquence importante est que les schémas numériques doivent permettre à la solution en
un nœud d’être influencée par ses voisins, de manière que l’information se propage dans toutes les directions.
C’est la raison pour laquelles les formules centrées (5.15), (5.17) et (8.3) sont particulièrement utilisées pour les
problèmes elliptiques.

7.2 Implémentation des conditions aux limites, et leur influence sur la précision
globale du schéma
Cette partie est consacrée à l’implémentation pratique des conditions aux limites dans un schéma aux différences.
Elle donne aussi un aperçu des problèmes que peut poser le traitement des frontières, et présente des notions
sur l’utilisation de schémas décentrés spéciaux près des bords. On montre aussi la chute de la précision globale
de la solution qui peut être observée si on n’est pas rigoureux dans l’implémentation des schémas de bords.
Considérons le problème de Helmholtz unidimensionnel
λu(x) − u00 (x) = f (x) dans Ω =]0, 1[, (7.4a)

avec les conditions aux limites de Neumann-Dirichlet non homogènes

u0 (0) = g− et u(1) = g+ . (7.4b)

Introduisons le maillage
Ωh = {xi = (i − 1) h, 1 ≤ i ≤ N } avec h = 1/(N − 1). (7.5)
On a donc N inconnues, les ui . Afin d’écrire un système discret complet, il faut disposer de N équations. En
règle générale :
• (N − 2) équations sont obtenues à partir de l’EDP (7.4a) aux points intérieurs i = 2, . . . , N − 1.
• 2 équations sont obtenues à partir des conditions aux limites (7.4b).
• On combine ces équations de manière à éliminer les valeurs au bord u1 et uN : on obtient alors un système
linéaire de taille N − 2 × N − 2 pour les inconnues u2 , . . . , uN −1 aux points intérieurs.
Au point courant xi , i = 2, . . . , N − 1, l’approximation par différences finies de (7.4a) s’écrit
0
λui − δxx ui = fi , (7.6)
0
où δxx est l’opérateur centré défini par (5.17). Cette formule est en O(h2 ), et fait intervenir toutes les inconnues
u1 , . . . , uN . Afin d’éliminer u1 et uN du système, il faut combiner cette équation avec les conditions aux limites.
Discrétisation de la condition de Dirichlet

Commençons par le plus facile : en i = N − 1, Eq. (7.6) s’écrit

uN −2 − 2 uN −1 + uN
λuN −1 − = fN −1 , (7.7)
h2
et la condition aux limites (7.4b) en x = 1 donne
uN = g+ . (7.8)
On peut donc éliminer uN de (7.7) qui devient
uN −2 − 2 uN −1 g+
λuN −1 − 2
= fN −1 + 2 . (7.9)
h h
Il est aisé de vérifier que cette formule est en O(h2 ). L’imposition de conditions de Dirichlet ne modifie donc
pas l’erreur de troncature locale en i = N − 1. Comme on l’a vu dans la section 5.4.2 pour le cas λ = 0, le
schéma démontrerait une précision globale en O(h2 ) si des conditions de Dirichlet étaient imposées à chacun
des deux bords.
 7.2 Implémentation des conditions aux limites, et leur influence sur la précision globale du schéma 51

Discrétisation de la condition de Neumann

La présence d’une condition de Neumann en x = 0 présente quelques difficultés. En i = 2, Eq. (7.6) s’écrit
u3 − 2 u2 u1
λu2 − 2
= f2 + 2 . (7.10)
h h
Maintenant, pour connaitre u1 , il va falloir discrétiser la condition de Neumann (7.4b) en x = 0 par différences
finies. On ne peut utiliser de formule centrée en x = 0 puisque cela ferait intervenir un nœud qui sort du domaine
Ω (cf. Fig 7.2). On va donc utiliser des formules décentrées.

Fig. 7.2. Maillage près de la frontière x = 0.

Exemple 7.3 (Condition de Neumann en O(h)).

Une première méthode utilise la formule

u2 − u1
u0 (0) = , (7.11)
h
qui possède une erreur de troncature en O(h), et d’où on tire

u1 = u2 − h g− . (7.12)

L’équation (7.10) devient alors

u3 − u2 g−
λu2 − = f2 − . (7.13)
h2 h
Le schéma complet s’écrit donc :
• Pour u2 , . . . , uN −1 : résoudre Eq. (7.6) pour i = 3, . . . , N − 2, Eq. (7.9) et Eq. (7.13).
• Ensuite, u1 et uN sont donnés respectivement par Eqs. (7.12) et (7.8).
On a ainsi éliminé u1 du système avec une formule de bord du premier ordre. Que devient l’erreur de troncature
du schéma, en particulier au point i = 2 ? D’après la définition 5.8, cette erreur s’écrit

τ2 = (Eq. (7.13)) − (λu2 − u002 − f2 ) , (7.14a)

−u3 + u2 + h g− + h2 u002
= . (7.14b)
h2
En utilisant les développements limités au nœud i = 2 :

h2 00 h3 (3)
u3 = u2 + hu02 + u + u2 + O(h4 ), (7.15a)
2 2 6
2
h (3)
g− ≡ u01 = u02 − hu002 + u2 + O(h3 ), (7.15b)
2
l’erreur de troncature devient
u002
τ2 = − + O(h), (7.16)
2
52 7 Les problèmes elliptiques

ce qui montre que le schéma n’est pas consistant en i = 2 !

En conséquence, même si le schéma démontre une consistance du second ordre aux autres points du mail-
lage, on observe une perte de la précision globale du schéma. Il s’agit là d’une manifestation d’un phénomène
caractéristique des équations elliptiques qui a été énoncé dans le paragraphe 7.1, c.à.d. une perturbation de la
solution en un point entraîne un changement de la solution dans la totalité du domaine. On observe alors que
la solution est seulement précise en O(h) dans la norme l2 , définie par (5.49b). Pour cet exemple simple, ce
dernier resultat peut être démontré mathématiquement, cf. [11] pour le cas λ = 0. Nous nous contenterons de
le vérifier numériquement en TD.

Exemple 7.4 (Condition de Neumann en O(h2 )).

La perte de précision de la méthode est une conséquence de l’utilisation de la formule au premier ordre (7.11).
Pour y remédier, on peut utiliser la formule suivante,
−u3 + 4 u2 − 3 u1
u0 (0) = , (7.17)
2h
qui est une formule décentrée au second ordre, qui ne fait intervenir que des points de Ωh (cf. Exercice 5.5). La
discrétisation de la condition de Neumann (7.4b) donne
1
u1 = (4u2 − u3 − 2hg− ) . (7.18)
3
Au point i = 2, la formule courante (7.6) devient
2 2
λu2 − 2
(u3 − u2 ) = f2 − g− . (7.19)
3h 3h
Le schéma complet s’écrit maintenant :
• Pour u2 , . . . , uN −1 : résoudre Eq. (7.6) pour i = 3, . . . , N − 2, Eq. (7.9) et Eq. (7.19).
• Ensuite, u1 et uN sont respectivement donnés par Eqs. (7.18) et (7.8).
Examinons, comme pour l’exemple précédent, la consistance en i = 2 :
2 ³ ´
τ2 = 2
−u3 + u2 + hg− + 23 h2 u002 , (7.20a)
3h
2 (3)
= u2 h + O(h2 ), (7.20b)
9
en utilisant les développements limités (7.15). Le schéma est maintenant consistant au premier ordre en i = 2.
Les calculs numériques montrent que la solution est globalement au second ordre dans la norme l2 .

Sur cet exemple simple, on voit le problème posé par l’implémentation des conditions aux limites, qui peuvent
nécessiter l’utilisation de schémas décentrés près du bord. Pour le schéma centré au second ordre (7.6), seule
la condition de Neumann a posé problème. Néanmoins, des difficultés analogues peuvent se poser si on veut
utiliser des schémas aux différences plus précis, par exemple :
µ 2 ¶
∂ u −ui+2 + 16 ui+1 − 30 ui + 16 ui−1 − ui−2
= , (7.21)
∂x2 i 12 h2

qui est précis en O(h4 ), mais fait intervenir 5 voisins ! L’utilisation de cette formule nécessiterait de construire
des schémas spéciaux près des bords. Ces schémas de bord influent sur la précision globale du schéma, mais
peuvent aussi influer sur la stabilité de la solution pour les problemes d’evolution.

7.3 Résolution de l’équation d’Helmholtz unidimensionnelle par l’algorithme de

Thomas
Le schéma complet de l’exemple 7.3 s’écrit
³ 1´ 1 g−
λ + 2 u2 − 2 u3 = f2 − , (7.22a)
h h h
1 ³ 2´ 1
− 2 ui−1 + λ + 2 ui − 2 ui+1 = fi , pour i = 3, . . . , N − 2, (7.22b)
h h h
1 ³ 2´ g+
− 2 uN −2 + λ + 2 uN −1 = fN −1 + 2 , (7.22c)
h h h
 7.3 Résolution de l’équation d’Helmholtz unidimensionnelle par l’algorithme de Thomas 53

après avoir éliminé les valeurs au bord :

u1 = u2 − hg− , (7.23a)
uN = g+ . (7.23b)

Le système (7.22) s’écrit en toute généralité sous la forme

ai ui−1 + bi ui + ci ui+1 = di , pour i = 2, . . . , N − 1, (7.24)

soit, matriciellement
    
b2 c2 u2 d2
 a3 b3 c3   
u3  d3
    
 · · ·   
·  ·
    
 ai bi ci  =
ui .di (7.25)
    
 · · ·   
·  ·
    
 aN −2 bN −2 cN −2   uN −2   dN −2 
aN −1 bN −1 uN −1 dN −1

Le système de l’exemple 7.4 peut aussi se mettre sous cette forme.

Ce système tridiagonal peut se résoudre de manière très efficace (et très simple à programmer) par
l’algorithme de Thomas, qui est un cas particulier de l’algorithme de factorisation LU pour les systèmes tridi-
agonaux. Cette méthode directe s’effectue en deux étapes successives :
• Premiere étape (forward sweep) : On effectue une élimination de Gauss, dans le sens des i croissants, pour
éliminer les ai et mettre les bi égaux à 1 : le système (7.25) devient
    0 
1 c02 u2 d2
 0 1 c03   u3   d03 
    
 · · ·  ·   · 
     
 0 1 c0i   ui  =  d0  (7.26)
    i 
 · · ·  ·   · 
     
 0 1 c0N −2   uN −2   d0N −2 
0 1 uN −1 d0N −1
avec
c2 d2
c02 = , d02 = , (7.27a)
b2 b2

et, pour i = 3 % N − 1 croissants :

ci di − ai d0i−1
c0i = , d0i = , (7.27b)
pivot(i) pivot(i)

avec

pivot(i) = bi − ai c0i−1 . (7.27c)

• Seconde étape (backward sweep) : après la première étape, la dernière ligne du système s’écrit

uN −1 = d0N −1 . (7.28a)

On peut alors remonter le système dans le sens des i décroissants pour calculer la solution ui de proche en
proche. Cette récurrence inverse est donnée par la formule

ui = d0i − ui+1 c0i , (7.28b)

pour i = N − 2 & 2 décroissants.

Cet algorithme, particulièrement peu coûteux en temps de calcul, est énormement utilisé en CFD, car il est à
la base de bon nombre d’algorithmes de résolution d’EDP multidimensionnelles sur maillages cartésiens.
8 Notions sur les schémas numériques multidimensionnels

8.1 Le schéma à 5 points pour l’opérateur d’Helmholtz

On considère la résolution de l’équation de Helmholtz
µ 2 ¶
∂ u ∂2u
λu − + 2 = f (x, y), λ ≥ 0, (8.1)
∂x2 ∂y

dans le domaine bidimensionnel Ω =]0, 1[2 par la méthode des différences finies. Dans ce but, on introduit le
maillage

Ωh = { (xi , yj ), 1 ≤ i ≤ Nx , 1 ≤ j ≤ Ny } , (8.2a)

dont les nœuds sont définis par

xi = (i − 1)∆x , yj = (j − 1)∆y , (8.2b)

avec les pas de discrétisation

1 1
∆x = , ∆y = . (8.2c)
Nx − 1 Ny − 1

Au point courant (i, j), on définit le schéma à 5 points pour (8.1) comme

Fig. 8.1. Le schéma à 5 points pour un maillage régulier.

¡ 0 0
¢
λ uij − δxx uij + δyy uij = fij , (8.3)

0
où δxx 0
et δyy sont les opérateurs de différences finies 1D définis en (5.17), pour ∂/∂x2 et ∂/∂y 2 respectivement.
Le schéma s’écrit aussi
µ ¶
ui+1,j − 2 ui,j + ui−1,j ui,j+1 − 2 ui,j + ui,j−1
λ uij − + = fij . (8.4)
∆x2 ∆y 2

Il est appelé schéma à 5 points car il fait intervenir le nœud courant (i, j) et ses 4 voisins le long des lignes du
maillage (cf. Fig. 8.1). Ce schéma possède une erreur de troncature en O(∆x2 , ∆y 2 ).
56 8 Notions sur les schémas numériques multidimensionnels

Fig. 8.2. Détail du maillage près de la paroi x = 0.

Les conditions aux limites se traitent de manière analogue au cas 1D. Considérons, par exemple, la condition
de Neumann ¯
∂u ¯¯
= g − (y), (8.5)
∂x ¯x=0
0
sur la paroi gauche du domaine Ω (cf. Fig. 8.2). Au nœud (2, j), il faut utiliser un schéma special pour δxx u2,j .
En utilisant la méthode décrite dans l’exemple 7.4, le schéma à 5 points s’écrit au point (2, j) :
2 u3,j − 2 u2,j 2
λ u2,j − 0
− δyy u2,j = f2,j − g− . (8.6)
3∆x2 3∆x j
Pour un nœud courant (i, j), il est commode d’introduire les notations géographiques du schéma à 5 points,
représentées sur la figure 8.3, avec P = (i, j), N = (i, j + 1) (North), etc . . . . Le schéma (8.4) s’écrit alors pour
un nœud suffisament loin des bords :

aS uS + aW uW + aP uP + aE uE + aN uN = fP , (8.7)

avec
1 1 2 2
aS = − , aW = − , aP = + + λ, (8.8a)
∆y 2 ∆x2 ∆x2 ∆y 2
1 1
aE = − , aN = − 2 , fP = fij . (8.8b)
∆x2 ∆y
Cette notation sera utile pour décrire la résolution itérative du schéma à 5 points dans le paragraphe 8.5.

Fig. 8.3. Notations géographiques en trois dimensions.

8.2 Les schémas implicites pour les équations paraboliques : lien avec l’équation
de Helmholtz
On considère, par exemple, l’équation de la chaleur 2D :
 8.3 Perspectives sur la résolution numérique du schéma à 5 points 57

∂u
− ν∇2 u = f, (8.9)
∂t
résolue par le schéma d’Euler implicite (ou rétrograde, cf. Chap. 6)

un+1
ij − unij ¡ 0 n+1 0 n+1
¢ n+1
− ν δxx uij + δyy uij = fij . (8.10)
∆t
En regroupant au second membre les termes déja connus à l’instant t = (n + 1)∆t, on est alors amené à résoudre
le système linéaire µ ¶
1 n+1 ¡ 0 n+1 0 n+1
¢ 1 n+1
unij
uij − δxx uij + δyy uij = fij + , (8.11)
ν∆t ν ∆t
qui est une équation de Helmholtz (discrète) de type (8.1) avec λ = 1/(ν∆t). Ainsi, à chaque niveau de temps,
l’utilisation du schéma d’Euler implicite conduit à la résolution d’une équation elliptique de type Helmholtz.
Remarque 8.1. Une discrétisation de (8.9) par le schéma de Crank-Nicolson (cf. TD), ou plus généralement tout
schéma implicite, conduit à la résolution d’un problème de Helmholtz analogue à (8.11).

8.3 Perspectives sur la résolution numérique du schéma à 5 points

Afin de résoudre le système (8.3), il est intéressant de remarquer qu’il peut s’écrire comme un système linéaire
lorsqu’on a arrangé les inconnues aux nœuds,

{ uij ; i = 1, . . . , Nx , j = 1, . . . , Ny } (8.12)

sous la forme du vecteur

u = {Ul , l = 1, . . . , Nx Ny } , (8.13a)

en posant

Ul = uij , avec l = (j − 1)Nx + i. (8.13b)

On fait de même pour le second membre f = {fl }. Cet arrangement des inconnues est connu sous le nom d’ordre
lexicographique, et il correspond à faire un balayage du domaine ligne par ligne dans le sens des i croissants
(cf. Fig. 8.4). Le schéma (8.3) peut alors s’écrire sous la forme matricielle

Fig. 8.4. Ordre lexicographique pour les nœuds intérieurs d’un maillage cartésien.

A u = f, (8.14)

où A est une matrice de taille Nx Ny × Nx Ny très creuse, puisqu’elle ne comporte que 5 éléments non nuls par
lignes. En reprenant les notations géographiques de la section 8.2, le système (8.14) est représenté sur la figure
8.5.
La résolution efficace de systèmes linéaires qui sont issus, tels (8.14), de la discrétisation des EDP est
un point crucial de l’efficacité d’une méthode numérique. En effet, la résolution de ces systèmes est une des
58 8 Notions sur les schémas numériques multidimensionnels

Fig. 8.5. Système linéaire correspondant au schéma à 5 points (8.3) ordonné lexicographiquement. Les 5 diagonales
hachurées représentent les élements non nuls de la matrice.

parties de l’algorithme les plus coûteuses en temps de calcul : par exemple, on observe que dans une simulation
d’écoulements de fluides incompressibles, la résolution de système linéaires occupe au mieux 60% du temps de
calcul, et qu’une résolution mal adaptée ou peu optimisée peut faire multiplier ce temps de calcul par 10 !
On peut diviser les méthodes de résolutions de systèmes linéaires issus de la mécanique en trois catégories :
A) Les méthodes directes : elles comprennent les méthodes bien connues de Gauss et de factorisation LU, dont
l’algorithme de Thomas (cf. section 7.3) est un cas particulier adapté aux matrices tridiagonales.
B) Les méthodes itératives, qu’elles soient de base comme les algorithmes de Jacobi, Gauss-Seidel, S.O.R., etc
. . . , qui seront presentees dans le paragraphe 8.5), ou plus sophistiquées, comme celles issues de la théorie
de l’optimisation telle la méthode du gradient conjugué du chapitre 2.
C) On rappelle que le système (8.14) est de taille Nx Ny ×Nx Ny . Son équivalent en 3D serait de taille Nx Ny Nz ×
Nx Ny Nz . Pour des problèmes de grandes taille (i.e. grilles tres fines, i.e. Nx Ny Nz très grands), on ne peut
stocker la totalité de la matrice A dans la mémoire de l’ordinateur. On doit alors se ramener à ne garder
uniquement que les éléments non-nuls de A et leur emplacement dans cette matrice : il s’agit de méthodes
de stockage de systèmes creux. Ce type d’adressage indirect complique la programmation et dégrade les
performances des supercalculateurs.
Dans de nombreux cas, il est intéressant d’exploiter les caractéristiques de l’EDP (i.e. séparabilite des
opérateurs) et du maillage (i.e. maillage cartésien ou structuré) pour se ramener à ne résoudre qu’une
série de problèmes unidimensionnels. Ces problèmes 1D sont typiquement résolus par l’algorithme direct de
Thomas. Parmi ces méthodes de résolution, on cite :
• Les méthodes du type “solveur de Poisson rapide” (hors programme)
• Les méthodes de “splitting” du type ADI, qui seront présentées dans le paragraphe 8.4.
Ces dernières méthodes sont bien adaptées aux EDP en géométries simples, telle l’équation de la chaleur
(8.9) dans un domaine rectangulaire.

8.4 La méthode ADI pour les équations paraboliques

8.4.1 Présentation de la méthode

La méthode ADI (“Alterning Direction Implicit”) et ses variantes sont des méthodes très répandues pour résoudre
les équations de type advection-diffusion en géométrie simple. Dans cette partie, la méthode ADI sera décrite
pour l’équation parabolique de la chaleur (8.9) dans Ω =]0, 1[2 . Elle peut aussi être interprétée comme une
méthode de résolution de l’équation de Poisson
µ 2 ¶
∂ u ∂2u
−ν + 2 = f,
∂x2 ∂y

lorsque t −→ ∞ : lorsque la solution est stationnaire, le processus d’évolution temporel peut être assimilé à un
processus itératif.
 8.4 La méthode ADI pour les équations paraboliques 59

Dans le paragraphe 8.2, on a vu que la résolution de (8.9) par le schéma d’Euler implicite entrainait la
résolution du problème d’Helmholtz 2D :

un+1
ij − unij ¡ 0 n+1 0 n+1
¢ n+1
− ν δxx uij + δyy uij = fij . (8.15)
∆t
pour passer du temps tn = n∆t au temps tn+1 = (n + 1)∆t. On est alors amené à résoudre un système linéaire
2D de la forme (8.14), car tous les termes de l’opérateur de Laplace ∇2 u = ∂ 2 u/∂x2 + ∂ 2 u/∂y 2 sont pris de
manière implicite.
En contraste, la méthode ADI permet de découpler en deux étapes l’opération tn −→ tn+1 de manière que,
à chaque étape, seuls les termes associés à une direction spatiale particulière sont traités implicitement. De cette
façon, on se ramène à calculer uniquement des problèmes de Helmholtz 1D, qui sont usuellement résolus par
l’algorithme de Thomas, cf. paragraphe 7.3. Ces étapes successives sont schématisés par :
A) Première étape, tn −→ tn+ 12 = (n + 21 )∆t : schéma implicite en x, et explicite en y,

n+ 12 ³ ´
uij − unij 0 n+ 1 0 n n
− ν δxx uij 2 + δyy uij = fij . (8.16)
∆t/2

B) Deuxième étape, tn+ 12 −→ tn+1 : schéma explicite en x, et implicite en y,

n+ 12 ³ ´
un+1
ij − uij 0 n+ 1 0 n+1 n+1
− ν δxx uij 2 + δyy uij = fij . (8.17)
∆t/2
1 1
Examinons l’ordre de consistance du schéma ADI : en formant 2 Eq. (8.16) + 2 Eq. (8.17), on obtient
µ ¶
un+1
ij − unij 0 n+ 1 1 0 ¡ n+1 ¢ 1 ¡ n+1 ¢
− ν δxx uij 2 + δyy uij + unij = n
fij + fij . (8.18)
∆t 2 2

En utilisant les développements suivants au temps (n + 12 )∆t :

µ ¶n+ 12 µ ¶2 µ ¶n+ 12 µ ¶3 µ ¶n+ 21
n n+ 12 ∆t ∂u 1 ∆t ∂2u 1 ∆t ∂3u
u =u − + − + O(∆t4 ),
2 ∂t 2! 2 ∂t2 3! 2 ∂t3
µ ¶n+ 12 µ ¶2 µ ¶n+ 12 µ ¶3 µ ¶n+ 21
n+1 n+ 12 ∆t ∂u 1 ∆t ∂2u 1 ∆t ∂3u
u =u + + + + O(∆t4 ),
2 ∂t 2! 2 ∂t2 3! 2 ∂t3
on obtient
µ ¶n+ 12 µ ¶n+ 21
un+1 − un ∂u ∆t2 ∂3u
= + + O(∆t3 ),
∆t ∂t 24 ∂t3

et pour φ = uij , fij ,

µ ¶n+ 12
1 ¡ n+1 ¢ 1 ∆t2 ∂2u
φ + φn = φn+ 2 + + O(∆t3 ).
2 8 ∂t2

Maintenant, en sachant que le schéma à 5 points est précis en O(∆x2 , ∆y 2 ), l’équation de troncature (8.18)
sՎcrit : ̵
µ ¶n+ 21 ¶n+ 12 µ 2 ¶n+ 12 !
∂u ∂2u ∂ u n+ 1
−ν + = fij 2 + O(∆t2 , ∆x2 , ∆y 2 ). (8.20)
∂t ij ∂x2 ij ∂y 2 ij
Le schéma ADI est consistant à l’ordre 2 en espace et en temps. Une étude de stabilité montrerait que ce schéma
est inconditionnellement stable.

8.4.2 L’algorithme de la méthode ADI

Pour la première étape, l’équation (8.16) peut s’écrire

n+ 12 0 n+ 12
uij − ν∆t
2 δxx uij = unij + ν∆t
2
0 n
δyy uij + ∆t
2
n
fij . (8.21)

Cela correspond à résoudre, pour chaque ligne J du maillage, J ∈ {j / j = 2, . . . , Ny − 1} , cf. Fig. 8.6, le
problème de Helmholtz 1D :
60 8 Notions sur les schémas numériques multidimensionnels
n+ 12 n+ 12 n
uJ − ν∆t
2 H uJ = unJ + ν∆t
2 [ uyy
J ] +
∆t
2 f nJ , (8.22a)

où la matrice H est donnée par (5.40), et

   
u1,J f1,J
   
uJ =  ...  , f J =  ...  . (8.22b)
uNx ,J fNx ,J

n
Le vecteur [ uyy yy n
J ] = (uJ (i), i = 1, . . . , Nx ) représente la dérivée seconde de u par rapport à y, qui est calculée
colonne par colonne, i.e.

ui,J+1 − 2 ui,J + ui,J−1

uyy 0
J (i) ≡ δyy ui,J = . (8.22c)
∆y 2

Enfin, on n’oublie pas d’associer à (8.22a) les conditions de Dirichlet sur les bords gauche (x = 0) et droit
(x = 1) du domaine :

n+ 1 n+ 1
u1,J 2 = gJW , uNx ,J
2
= gJE . (8.22d)

De même, pour la seconde étape, Eq. (8.17) s’écrit

n+ 12 n+ 12
un+1
ij − ν∆t
2
0 n+1
δyy uij = uij + ν∆t
2
0
δxx uij + ∆t
2
n+1
fij . (8.23)

Cela correspond à résoudre Nx − 2 problèmes d’Helmhlotz 1D sur chaque colonne du maillage, avec conditions
de Dirichlet sur les parois supérieure et inférieure du domaine.

Fig. 8.6. La méthode ADI.

8.5 Les méthodes itératives de base pour le schéma à 5 points

On va s’intéresser à la résolution du système linéaire

A u = f, (8.24)

correspondant à la discrétisation de l’équation de Helmholtz 2D par le schéma à 5 points (cf. section 8.1). On
va décrire les méthodes stationaires de Jacobi, Gauss-Seidel et S.O.R., qui utilisent la décomposition suivante
de A :
A = D − L − U, (8.25)
 8.5 Les méthodes itératives de base pour le schéma à 5 points 61

Fig. 8.7. Décomposition de la matrice A en trois sous-matrices disjointes −Nl , M, et −Nu . Le principe des méthodes
itératives de ce paragraphe est de n’inverser que la sous-partie M de A : M doit donc être une matrice “facile” à inverser.
Dans la déconposition LU que nous utilisons, M est la diagonale principale de A, qui est triviale à inverser.

où D, L1 et U 2 sont les sous-matrices de A contenant respectivement les termes diagonaux, supérieurs et in-
férieurs par rapport à la diagonale principale (cf. Fig 8.7). Dans le cas du schéma à 5 points, cette décomposition
appliquée à l’équation (8.7) s’écrit :

aS uS + aW uW + aP uP + aE uE + aN uN = fP . (8.26)
| {z } | {z } | {z }
−Lu Du −Uu

La plus simple (simpliste ?) des méthodes itératives est la méthode de Jacobi, qui s’écrit

Dum+1 = f + (L + U ) um , (8.27)

où l’indice m se réfère au nombre d’itérations de la méthode. Dans le cas d’un maillage cartésien, une itération
m −→ m + 1 correspond à un balayage (sweep) du maillage dans l’ordre lexicographique (cf. section 8.3).
Lorsqu’elle est appliquée à l’équation (8.26), la méthode de Jacobi s’écrit
1 £ ¤
um+1
P = fP − (aS um m m m
S + aW uW ) − (aE uE + aN uN ) . (8.28)
aP
Cette méthode est aujourd’hui très peu utilisée car elle se révèle très coûteuse en temps CPU : elle se révèle
être moins efficace qu’une méthode directe.
Une amélioration de la méthode de Jacobi consiste à utiliser les valeurs de um+1 dans le second membre dès
qu’elles sont disponibles : c’est la méthode de Gauss-Seidel, qui s’écrit :

Dum+1 = f + Lum+1 + Uum , (8.29)

soit, pour le schéma à 5 points,

1 £ ¡ ¢ ¤
um+1
P = fP − aS um+1
S + aW um+1
W − (aE um m
E + aN uN ) , (8.30)
aP

car, lorsqu’à l’itéré m + 1 on arrive au point P, les valeurs um+1

S et um+1
W sont déja connues, mais pas encore
m m
uE et uN . La méthode de Gauss-Seidel assure en général une convergence deux fois plus rapide que la méthode
de Jacobi.
La méthode S.O.R.3 apporte une amélioration significative de la méthode de Gauss-Seidel en évaluant um+1
comme la moyenne pondérée de um et (um+1 )GS , cette dernière quantité étant le résultat de l’itération de
Gauss-Seidel (8.30). La méthode S.O.R. s’écrit
¡ m+1 ¢
um+1
P = (1 − ω) umP + ω uP GS
, (8.31a)

soit
ω£ ¡ ¢ ¤
um+1
P = (1 − ω) um
P + fP − aS um+1
S + aW um+1
W − (aE um m
E + aN uN ) . (8.31b)
aP

dire un mot sur ω .

1
L : Lower (inférieur).
2
U : Upper (supérieur).
3
Successive Over-Relaxation.
9 Les problèmes hyperboliques

9.1 Généralités
Reprenons le modèle des équations de bilan de la mécanique du paragraphe 4.2, c.à.d. l’équation d’advection-
diffusion :
∂u ∂u ∂2u
+ a − ν = 0. (9.1)
∂t | {z∂x } ∂x2 }
| {z
a6=0, advection ν≥0, diffusion

On rappelle que dans le cas ν 6= 0, cette équation est de type parabolique, et pour ν = 0 elle prends la forme
de l’équation d’advection, ou de transport :
∂u ∂u
+a = 0. (9.2)
∂t ∂x
qui est de type hyperbolique. Comme il a été vu dans le paragraphe 4.2, une solution de Fourier de cette équation
s’écrit
b0 e i k (x−at) ,
u(x, t) = u
ce qui correspond à une propagation de la condition initiale à la vitesse a sans déformation ou dissipation
visqueuse.

Remarque 9.1. Revenons au cas ν 6= 0, en introduisant le nombre de Péclet :

|a|
Pe = , (9.3)
ν
qui exprime le rapport entre les forces d’advection et de diffusion, et qui est à rapprocher du nombre de Reynolds
Re des équations de Navier-Stokes. Lorsque Pe À 1 (ou |a| À ν), c’est-à-dire lorsque la dissipation est très
faible, des effets hyperboliques apparaissent dans la solution, meme si l’équation (9.1) est, stricto senso, de type
parabolique. Ce mélange de comportement parabolique et hyperbolique est typique de l’équation d’advection-
diffusion et, plus généralement, des équations de la mécanique des fluides. C’est la raison pour laquelle les
méthodes du chapitre 6 pour les problèmes paraboliques deviennent moins performantes dans le cas Pe À 1
(ou Re À 1), et que certaines des méthodes pour les équations hyperboliques que nous verrons ci-après, sont
utilisées en CFD pour résoudre les équations de Navier-Stokes.

Exemple 9.2. Les équations d’Euler instationnaires pour la dynamique des gaz (non visqueux). En aéronautique,
on fait souvent l’hyp que le fluide est non-visqueux ... Pour calculer les écoulements autour d’ailes d’avion, on
applique généralement une variante des schémas qui seront discutés dans ce chapitre à la résolution des équations
d’Euler, qui sont non-linéaires.

Exemple 9.3 (Vibration de torsion d’une poutre de section circulaire). Les équations hyperboliques sont om-

Fig. 9.1. Torsion d’une poutre de section S.

64 9 Les problèmes hyperboliques

niprésentes dans les problèmes de vibrations mécaniques. Par exemple, si on considère le problème de torsion
représenté sur la figure 9.1, l’équation vérifiée par l’angle de rotation ϕ de la section S est
∂2ϕ ∂2ϕ µ
2
− c2 = 0, avec c2 = , (9.4)
∂t ∂x2 ρ
où µ est le coefficient de Lamé et ρ la p
masse volumique de la poutre. Cette équation, appelée équation des ondes,
est de nature hyperbolique, et c = ± µ/ρ est la célérité des ondes de torsion.
Plus généralement, les équations de la mécanique qui modélisent des phénomènes de conservation sont de nature
hyperbolique, et non-linéaires. Leur solution est caractérisée par des paquets d’ondes qui se propagent sans perte
d’amplitude ni de phase. La résolution numérique de ces problèmes est délicate car il est important d’utiliser un
schéma numérique dont les erreurs altèrent le moins possible l’amplitude de ces ondes (phénomène de dissipation
numérique), et leur célérité (phénomène de dispertion numérique).
Afin de se familiariser à la résolution de problèmes hyperboliques, nous allons examiner une variété de
schémas explicites pour l’équation de transport (9.2) tels qu’ils ont été historiquement proposés, et examiner
leurs performances. Ces considérations plus ou moins théoriques seront approfondies (et éclaircies !) en TP et
TD.
Les notions de stabilité (analyse de Von Neumann, paragraphe 6.3) et de consistance (paragraphe 5.4.1)
sont, bien sûr, toujours valides pour les équations hyperboliques.

9.2 Solution générale et lignes caractéristiques

On considère donc, dans la suite de ce chapitre, l’équation de transport
∂u ∂u
+a = 0. (9.5)
∂t ∂x
avec a > 0, et la condition initiale à t = 0
u(x, 0) = u0 (x). (9.6)
Il est aisé de vérifier qu’une solution générale de (9.5) s’écrit sous la forme :
u(x, t) = f (x − at) = f (ξ), (9.7)
avec ξ = x − at, et f une fonction pour l’instant quelconque. Une conséquence importante est que la solution
u(x, t) reste constante sur la droite d’équation x − at = Cte , appellee ligne caractéristique dans la théorie des
EDP.

Fig. 9.2. Lignes caractéristiques pour l’équation d’advection (9.5).

Il est possible de déterminer complètement la solution en utilisant la condition initiale (9.6). Soit ξ? > 0 et
plaçons nous sur la ligne caractéristique d’équation
x − at = ξ? ,
représentée sur la figure 9.2. Le point courant (x? , t? ) tel que x? − at? = ξ? est sur cette ligne de même que, au
temps initial t = 0, le point (ξ? = x? − at? , 0). La solution vérifie donc
u(x? , t? ) = u(x? − at? , 0) ≡ u0 (x? − at? ),
ce qui signifie que l’information se propage à partir de l’état initial u0 le long des caractérisitiques, et le problème
(9.5)-(9.6) admet finalement pour solution
u(x, t) = u0 (x − at). (9.8)
 9.4 Le schéma “upwind” (décentré) et la condition CFL 65

9.3 Le schéma explicite centré est-il stable ?

Le schéma le plus simple que l’on pourrait utiliser pour résoudre l’équation de transport (9.5), avec a de signe
quelconque, est le schéma explicite centré défini par

un+1 − uni
i
+ a δx0 uni = 0, (9.9)
∆t
où l’opérateur δx0 est défini en (5.15). Il est facile de vérifier que l’erreur de troncature de ce schéma est en
O(∆t, ∆x2 ).
Exercice 9.4. Montrez que le facteur d’amplification g du schéma (9.9) s’exprime comme

g = 1 − i C sin k∆x,

= où C est le nombre de Courant C, défini par

∆t
C=a . (9.10)
∆x
Le schéma explicite est il stable ?

9.4 Le schéma “upwind” (décentré) et la condition CFL

Un premier remède à l’instabilité du schéma (9.9) consiste à décentrer l’approximation de la partie spatiale de
l’équation de transport suivant le sens de propagation de l’information (i.e. signe de a). Ce principe est commun
à tous les schémas dits “upwind” pour les problèmes hyperboliques.
• Cas a > 0 : schéma décentré amont (“upwind”)

un+1 − uni
i
+ a δx− uni = 0, (9.11a)
∆t

avec δx− uni défini en (5.13). Ce schéma s’écrit aussi

¡ ¢
un+1
i = uni − C uni − uni−1 , (9.11b)

avec C le nombre de Courant (9.10).

• Cas a < 0 : schéma décentré aval (“downwind”)

un+1 − uni
i
+ a δx+ uni = 0, (9.12a)
∆t

avec δx+ uni défini en (5.14). Ce schéma s’écrit aussi

¡ ¢
un+1
i = uni − C uni+1 − uni , (9.12b)

Ces deux schémas sont du premier ordre en espace et en temps. Une étude de la stabilité du schéma (9.11) pour
le cas a > 0 montre que le facteur d’amplification s’écrit

g = 1 − C (1 − cos k∆x) − i C sin k∆x. (9.13)

La condition de Von Neumann |g|2 ≤ 1 est vérifiée lorsque

C(1 − C)(1 − cos k∆x) ≥ 0.

Comme
a ∆t
C= ≥0 et 1 − cos k∆x ≥ 0,
∆x
la condition de stabilité devient
1 − C ≥ 0,
soit
∆t
C=a ≤ 1. (9.14)
∆x
66 9 Les problèmes hyperboliques

Remarque 9.5. Si a était négatif, la condition de Von Neumann ne serait pas vérifiée, et le schéma (9.11) serait
inconditionnellement instable.

Remarque 9.6. Pour le cas a < 0, la condition de stabilité du schéma (9.12) est :
∆t
−a ≤ 1. (9.15)
∆x
En conclusion, les schémas décentrés (9.11) et (9.12) sont stables sous la condition CFL (Courant-Friedrich-
Levy, 1928) :
∆t
|C| = |a| ≤ 1. (9.16)
∆x
Cette condition de stabilité est omniprésente pour les schémas hyperboliques explicites. On la retrouve aussi
énormément dans les méthodes numériques pour les écoulements de fluides, où les termes de convection (hyper-
boliques) sont très souvent pris en compte explicitement, cf. le logiciel FLUENT.

9.5 Schémas décentrés et dissipation numérique

Dans le paragraphe 9.3, il a été observé que le schéma explicite centré est inconditionnellement instable. S’il était
utilisé dans un calcul numérique, on observerait au cours de l’intégration en temps des oscillations parasites de
la solution qui iraient en s’amplifiant, jusqu’à atteindre des amplitudes “infinies” pour l’ordinateur (phénomène
d’“overflow”). Ces oscillations sont amplifiées car il n’y a pas de termes paraboliques (dissipatifs) pour les amortir.
Cette instabilité est due à l’approximation centrée du terme a ∂u/∂x, et elle est éliminée lorsqu’on utilise les
schémas upwind (9.11) ou (9.12) et que la condition CFL est vérifiée.
Dans le cas a > 0, par exemple, on peut interpréter l’effet de décentrage amont en écrivant δx− ui comme :
ui − ui−1 ui+1 − ui−1 ui+1 − 2 ui + ui−1
δx− ui ≡ = − , (9.17a)
∆x 2∆x 2∆x
∆x 0
= δx0 ui − δ ui . (9.17b)
2 xx
Le schéma upwind (9.11) s’écrit alors :

un+1 − uni
i
+ a δx0 uni − ² δxx
0
un = 0, (9.18a)
∆t | {z i}
terme dissipatif

avec

∆x
²=a > 0. (9.18b)
2
Le schéma upwind s’interprète alors comme le schéma explicite centré (9.9) auquel on ajoute (artificiellement)
un terme dissipatif qui va stabiliser la solution numérique, tout en préservant la consistance du schéma. Une
des conséquences de la présence de ce terme dissipatif est l’amortissement au cours du temps de la solution
numérique, qui est un effet non observé par la solution exacte. Le terme ² > 0 est appelé viscosité artificielle
(ou numérique) par analogie à la mécanique des fluides (terme dissipatif : ν∇2 u).
Remarque 9.7. La notion de dissipation numérique est valable pour d’autres schémas (cf. TP, TD). On l’observe
lors de la résolution numérique de toutes les EDP de la mécanique, mais elle est essentiellement prononcée pour
les EDP à haut caractère hyperbolique, par exemple l’équation d’advection diffusion (9.1) à haut nombre de
Péclet.

Remarque 9.8. Cette notion de dissipation numérique sera approfondie en TD. On introduira aussi la notion de
dispersion numérique.

9.6 Un schéma explicite du second ordre : le schéma de Lax-Wendroff

D’après la section précédente, une technique de stabilisation du schéma explicite centré (9.9) consiste à lui
rajouter un terme de dissipation numérique. Or, les essais numériques montrent que la dissipation artificielle
des schémas décentrés (9.11) et (9.12), qui sont du premier ordre, peut diminuer de façon excessive la précision
de la solution calculée (cf. TP, TD). Lax & Wendroff (1960) ont proposé une autre stabilisation du schéma
explicite centré, construit de la manière suivante :
 9.7 Commentaires et perspectives 67

On effectue premièrement un développement en temps de la solution, soit

∂u ∆t2 ∂ 2 u
u(xi , tn + ∆t) = u(xi , tn ) + ∆t (xi , tn ) + (xi , tn ) + O(∆t3 ). (9.19a)
∂t 2 ∂t2

D’après (9.5), on a

∂u ∂u ∂2u ∂2u
(x, t) = −a (x, t), 2
(x, t) = a2 2 (x, t). (9.19b)
∂t ∂x ∂t ∂x

On reporte ces expressions dans (9.19a) pour obtenir

∂u ∆t2 ∂ 2 u
u(xi , tn + ∆t) = u(xi , tn ) − a∆t (xi , tn ) + a2 (xi , tn ) + O(∆t3 ). (9.19c)
∂x 2 ∂x2

En tronquant ce développement au second-ordre et en approchant les dérivées spatiales par différences finies
centrées, on obtient :

∆t2 0 n
un+1 = uni − a∆t δx0 uni + a2 δ u . (9.19d)
i
2 xx i
Le schéma de Lax-Wendroff s’écrit alors, pour a > 0 ou a < 0,

un+1 − uni ∆t 0 n
i
+ a δx0 uni − a2 δxx ui = 0. (9.20)
∆t | {z2}
²>0 : viscosité numérique

Son facteur d’amplification a pour expression :

∆x
g = 1 − 2 C 2 sin2 k − i C sin k∆x. (9.21)
2
On déduit facilement que, quelque soit le signe de a, le schéma est stable sous la condition CFL (9.16). Il est
aisé de vérifier que l’erreur de troncature de ce schéma est en O(∆t2 , ∆x2 ).
Historiquement, lorsqu’il fut appliqué aux équations d’Euler, il fut le premier schéma qui a permis le calcul
précis d’écoulements autour d’ailes d’avion.

9.7 Commentaires et perspectives

Il existe, bien sûr, d’autres schémas explicites au second ordre que le schéma de Lax-Wendroff (ex : le schéma
saute-mouton, etc . . . ). Pour des solutions régulières, les essais numériques montrent que le schéma de Lax-
Wendroff donne les résultats les plus satisfaisants.
Néanmoins, dans certaines circonstances, les problèmes hyperboliques ont la particularité de posséder des
solution discontinues (appelées chocs ou ondes de choc). Pour les équations linéaires comme celles considérées
dans ce chapitre, les chocs sont uniquement présents lorsque la condition initiale est discontinue : ils se propagent
alors le long des caractéristiques.
Par contre, les équations non-linéaires peuvent présenter des chocs même si les données (i.e. conditions
initiales et aux limites) sont régulières. Par exemple, en mécanique des fluides compressibles, dans les écoulements
régis par les équations d’Euler, des ondes de discontinuité sont génerées par un avion en vol supersonique.
En presence de chocs, les schémas présentés dans ce chapitre montrent leurs faiblesses : en particulier, le
schéma de Lax-Wendroff exhibe de nombreuses oscillations au voisinage de la discontinuité....
 1.5
1

1
0.5
0.5

0 0

–0.5 –0.5

–1
–1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1
(a) Schéma Upwind (a) Schéma Upwind

1.5
1

1
0.5
0.5

0 0

–0.5 –0.5

–1
–1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1
(b) Schéma de Lax-Wendroff (b) Schéma de Lax-Wendroff

1.5
1

1
0.5
0.5

0 0

–0.5 –0.5

–1
–1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1
(c) Schéma MinMod (c) Schéma MinMod

Figure 1: Figure 2:

4 5
 1. Motivations & Objectifs (1/2)
Z Motivation de ce travail : Solveur précis & efficace pour géométries complexes (écoulement
visqueux incompressible, méthode volumes/différences fini(e)s)
Z (Un des) enjeux majeurs : Génération de grille pour géométries complexes
Grilles conformes : qui épousent la frontière des géométries complexes

Maillage curviligne Maillage non-structuré

V V
V
Maturité des solveurs structurés Génération de maillage très automatisée
Contrôle aisé de la densité de maillage V
X
Utilisation sur de nombreux types de
Génération de grille «chronophage» géométries
en terme de main d’œuvre X Absence de données structurées
⇒ coût CPU élevé

Y. CHENY & O. BOTELLA (LEMTA) La méthode LS-STAG CFM Grenoble 2008 2 / 19

1. Motivations & Objectifs (2/2)

Grilles non-conformes : méthodes aux frontières immergées (IB)

Méthodes cartésiennes : Frontières non-alignées avec la grille

Cellules solides ∈ Ωs ⇒ Pas de discrétisation

Cellules fluides ∈ Ωf ⇒ Discrétisation cartésienne habituelle

Cellules mixtes (ou « Cut-cells » ) ⇒ IB prise en compte dans la discrétisation

Point central des méthodes IB : Discrétisation dans les cut-cells

Y. CHENY & O. BOTELLA (LEMTA) La méthode LS-STAG CFM Grenoble 2008 3 / 19
 Partie III

Les éléments finis

10 Introduction à la méthode des éléments finis

10.1 Introduction et motivation

10.1.1 La méthode des résidus pondérés

Les met́hodes de différences finies pour résoudre une E.D.P consistent à approcher les opérateurs différentiels
de l’équation avec une formule aux différences et à chercher une approximation de la solution aux nœuds du
maillage. En tout rigueur, la solution n’est définie qu’aux points de ce maillages.
Au contraire, les méthodes de résidus pondérés supposent que la solution u(t, x) d’une EDP peut être
représentée sous la forme analytique :
N
X
u(t, x) ∼
= uh (t, x) ≡ aj (t)ϕj (x), (10.1)
j=1

où les aj (t) sont des coefficients à determiner et les ϕ(x) sont des fonctions analytiques connues. Les ϕ(x) sont
appelées fonctions de base.
Considérons l’approximation de l’EDP générique :
∂u
− L(u) − f = 0, (10.2)
∂t
ou L est un opérateur différentiel spatial1 et f un second-membre donné. Pour déterminer les N coeffcients
aj (t), on a besoin de construire N équations à partir de l’EDP (10.2). Pour cela, nous injectons l’approximation
de la solution (10.1) dans (10.2) et obtenons :

∂uh
− L(uh ) − f ≡ R(uh ), (10.3)
∂t
où R est le résidu de l’équation, que l’on va entreprendre d’annuler ou, du moins, rendre aussi petit que possible.
Pour cela, nous allons exprimer que l’intégrale du residu sur tout le domaine, pondéré par des fonctions test
ϕ(x),
e soit nulle : Z
R(uh ) ϕ
ei dV = 0, pour i = 1, . . . , N, (10.4)
Ω

ce qui nous fournit les N équations permettant de déterminer les coefficients aj (t), j = 1, . . . , N de la solution
approchée (10.1).
L’équation (10.4) est appelée formulation faible (ou variationnelle) de l’EDP (10.2). Suivant la forme des
fonctions test ϕ,
e on obtient différents types de méthodes, notamment :
A) Méthodes des sous-domaines (volumes finis)

Considérons une partition du domaine Ω en N sous-domaines (ou mailles) Ωi , i = 1, . . . N . Si les fonctions

test sont choisies comme : ½
1 dans Ωi ,
ϕ
ei = (10.5)
0 en dehors,
alors Eq. (10.4) s’écrit :
Z ³ ´
∂uh
− L(uh ) − f dV = 0, pour i = 1, . . . , N, (10.6)
Ωi ∂t

et exprime la conservation des quantités transportées dans une maille Ωi . C’est le point de départ de la
méthodes des volumes finis, où Ωi est souvent appelé volume de contrôle.

1
Par exemple, L = −∇2 .
72 10 Introduction à la méthode des éléments finis

B) Méthode Galerkin (Eléments finis)

Si les fonctions test sont choisies égales aux fonctions de base, i.e.

ϕ
ej = ϕj ,

on obtient alors la méthode de Galerkin. Cette méthode exprime que le résidu soit orthogonal à chaque
fonction test ϕj dans le produit scalaire L2 (Ω), c.à.d. orthogonal à l’espace

Vh = {ϕj , j = 1, . . . , N } ,

formé par les fonctions de base. Lorsque N → ∞ et que

Vh → V,

où V représente l’espace fonctionnel des solutions de l’EDP (10.2), alors la méthode de Galerkin exprime
que
uh → u,
la solution exacte de l’EDP.

C) Il existe d’autres méthodes de résidus pondérés, en particulier la méthode de collocation (hors programme).
Les approximations par sous-domaines et de Galerkin sont à l’origine des méthodes de volumes finis et des
éléments finis respectivement. Leur principal avantage sur les différences finies et de pouvoir s’appliquer à des
géométries arbitrairement complexes (cf. paragraphe suivant).
Lorsqu’on considère un maillage de type cartésien comme ceux que nous avons vu pour les différences finies
(i.e. Fig. 5.1), la méthode de volumes finis est très semblable aux différences finies.
Par contre, la méthode des éléments finis repose sur des techniques très spécifiques qui seront présentées
dans la suite de ce chapitre, et notamment :
• La construction de la formulation variationnelle (10.4) par méthode de Galerkin, permettant de fixer les
conditions aux limites et de déterminer l’espace fonctionnel V où chercher sa solution.
• La construction de l’espace d’éléments finis Vh tel que Vh ∼
= V , qui permet de représenter la solution discrète
(10.1) dans sa base d’éléments finis.
• La construction du problème discret à partir de la formulation variationnelle.

10.1.2 Les limitations d’un maillage de différences finies

On considère, par exemple, l’écoulement autour d’une aile d’avion représenté sur la figure 10.1.

Fig. 10.1. Ecoulement supersonique autour d’un profil d’aile d’avion. Dans cette approximation on considère une fluide
non visqueux compressible, modélisé par les équations d’Euler.

Maillage structuré (Différences finies)

Pour ce type de géométries, on considère un maillage structuré représente sur la figure 10.2. Il est possible de
trouver un changement de variable
(ξ, η) = f (x, y), (10.7)
entre le domaine physique (x, y) et le domaine de calcul (ξ, η) tel que, dans ce dernier repère, le maillage soit
de la forme représenté sur la figure 10.3.
Il n’est néanmoins pas très aisé de contrôler la distribution des nœuds du maillage, par exemple si on veut
mettre beaucoup de points dans la couche limite ou le sillage de l’aile. Ces maillages structurés sont donc limités
à des géométries relativement simples, où il est possible de trouver, analytiquement ou numériquent, une relation
du type (10.7).
 10.1 Introduction et motivation 73

Fig. 10.2. Maillage curviligne utilisé pour l’écoulement autour d’une aile d’avion.

Fig. 10.3. Représentation du maillage curviligne dans le domaine de calcul (ξ, η).

Fig. 10.4. Maillage non-structuré (triangulation) pour l’écoulement autour d’une aile d’avion.

Les maillage non-structurés

Ce type de maillage est très bien adapté aux approximations par élèments finis ou volumes finis. Dans
l’exemple de la figure 10.4, par l’utilisation d’éléments triangulaires, il est possible d’augmenter la densité du
maillage dans la couche limite de l’aile, sans mettre beaucoup d’éléments dans les zones éloignées du profil ou
rien d’important n’apparait.
74 10 Introduction à la méthode des éléments finis

Ce type de maillages est réalisé de manière (semi)-automatique par un logiciel de maillage appellé mailleur,
par exemple GAMBIT est le mailleur utilisé par le logiciel FLUENT. Avec quelques semaines de travail, vous
pouvez arriver à mailler le type de géométries représenté sur la figure 10.5.

Fig. 10.5. Maillage non-structuré (héxaèdral) pour une géométrie complexe tridimensionnelle (skud).

10.2 Les bases d’éléments finis

10.2.1 L’élément fini IP1 en une dimension

Un espace d’éléments finis Vh est une base de polynômes par morceaux permettant de représenter des fonctions
continues sur le domaine Ω. On prendra dans cet exemple Ω =]0, 1[.
Nous allons définir au préalable ce qu’est un polynôme par morceaux. Considérons d’abord le maillage de
Ω, tel que :
© ª
Ωh = xi ∈ Ω, i = 1, . . . , N , (10.8a)

avec

0 = x1 < x2 < · · · < xN −1 < xN = 1. (10.8b)

Un polynôme par morceaux p(x) sur Ωh est tel que :

• p est une polynôme sur [xi , xi+1 ],
• p est continu aux nœuds xi .
L’espace Vh1 des éléments finis IP1 est l’ensemble des polynômes par morceaux définis par
© ª
Vh1 = p, tel que p ∈ C 0 (0, 1), et p|[xi ,xi+1 ] ∈ IP1 pour 1 ≤ i ≤ N − 1. (10.9)

On a utilisé les définitions suivantes :

• IP1 est l’espace des polynômes de degré inférieur ou égal à 1 (i.e. , de la forme ax + b).
• C 0 (0, 1) est l’espace des fonctions continues sur ]0, 1[.
Un exemple de fonction de Vh1 est représenté sur la figure 10.6. L’espace Vh1 est un espace vectoriel de
dimension N . Il est alors très commode de représenter ses éléments dans la base canonique de Lagrange :
© ª
ϕj (x) ∈ Vh1 , j = 1, . . . , N , (10.10)

définie par :
 10.2 Les bases d’éléments finis 75

Fig. 10.6. Polynôme par morceaux appartenant à Vh1 .







 0 si x ≤ xj−1 ,





 x − xj−1


 pour x ∈ [xj−1 , xj ],
xj − xj−1
ϕj (x) = (10.11)


 xj+1 − x


 pour x ∈ [xj , xj+1 ],
 xj+1 − xj







 0 si x ≥ xj+1 .

Fig. 10.7. Représentation de la base de Lagrange de Vh1 .

Cette base est représentée sur la figure 10.7. Clairement, une fonction de base vérifie

φj (xi ) = δij , i = 1, . . . , N,

et un élément p de Vh1 s’écrit dans cette base comme

N
X
p(x) = p(xj )ϕj (x). (10.12)
j=1

c’est-à-dire que les coefficients du polynôme par morceaux sont égaux à la valeur de p aux nœuds du maillage.
Puisque ϕj (x) est non-nulle dans l’intervalle [xj−1 , xj+1 ] uniquement, elle s’exprime tres facilement dans les
éléments Aj−1 = [xj−1 , xj ] et Aj = [xj , xj+1 ] comme :

p(x) = p(xj−1 )ϕj−1 (x) + p(xj )ϕj (x), dans Aj−1 ,

et
76 10 Introduction à la méthode des éléments finis

p(x) = p(xj )ϕj (x) + p(xj+1 )ϕj+1 (x), dans Aj .

C’est dû au caractère local de la représentation dans la base canonique que l’approximation des opérateurs par
éléments finis est creuse (plus particulièrement, en 1D, les matrices sont bandes).

10.2.2 Approximation par éléments finis et espaces de Sobolev

Il est possible de représenter une fonction continue quelconque v(x) par son approximation polynômiale vh (x) ∈
Vh1 , telle que
N
X
vh (x) = v(xj )ϕj (x). (10.13)
j=1

On note parfois l’interpolant vh sous la forme :

vh = Πh1 v.

Fig. 10.8. Interpolation par éléments finis IP1 d’une fonction C 0 .

Un exemple d’interpolation de fonction continue est représenté sur la figure 10.8. Il est important de mesurer
la qualité de cette approximation polynômiale (i.e. , l’ordre de convergence de cette approximation). Ces mesures
sont effectuées dans les espaces de Sobolev L2 (Ω) et H 1 (Ω).

Rappels sur les espaces de Sobolev

Soit f une fonction de Ω ⊂ Rn à valeur dans R. L’espace vectoriel L2 (Ω) est l’ensemble des fonctions de carré
intégrables, i.e. Z
n o
2
L (Ω) = f tel que |f |2 dV < +∞ ,
Ω

qui admet pour produit scalaire Z

(f, g) = f g dV,
Ω
et pour norme
kf kL2 (Ω) = (f, f )1/2 .
L’espace vectoriel H 1 (Ω) est l’ensemble des fonctions admettant des dérivées de carré intégrable
n ¡ ¢n o
H 1 (Ω) = f tel que f ∈ L2 (Ω) et ∇f ∈ L2 (Ω) .

L’espace H 1 (Ω) ne possède pas de produit scalaire, mais possède une norme, définie par :
³ ´1/2
||f ||H 1 (Ω) = ||f ||2L2 (Ω) + ||∇f ||2(L2 (Ω))n .

Les fonctions de H 1 (Ω) sont par définition plus régulières que celle de L2 (Ω) , car les dérivées d’une fonction
de L2 peuvent être des distributions (i.e. , des éléments de H −1 (Ω)) ! On a en particulier l’inclusion
 10.3 Formulation variationnelle d’un problème aux limites 77

H 1 (Ω) ⊂ L2 (Ω).

Une conséquence de cette plus grande régularité est qu’on peut, contrairement aux fonctions L2 , définir la
valeur au bord ∂Ω (i.e. , la trace) des fonctions de H 1 (Ω). Un sous-espace de H 1 (Ω) très utile pour imposer
les conditions aux limites de Dirichlet dans la méthode des éléments finis est :
n o
H01 (Ω) = f ∈ H 1 (Ω) tel que f |∂Ω = 0 . (10.14)

Précision de l’interpolation IP1

Les fonctions de base ϕj (x) sont dérivables une fois seulement et on a ϕ0j ∈ L2 (Ω). Donc ϕj appartient à H 1 (Ω)
et on a l’inclusion :
Vh1 ⊂ H 1 (Ω).
On dispose en plus du résultat de convergence suivant :
Théorème 10.1. Soit v ∈ L2 (Ω), avec Ω =]0, 1[, telle que v 0 et v 00 sont aussi dans L2 (Ω). Alors on dispose
des estimations d’erreur suivantes :

||v − vh ||L2 (Ω) ≤ h2 ||v 00 ||L2 (Ω) ,

||v − vh ||H 1 (Ω) ≤ h ||v 00 ||L2 (Ω) ,

avec h = maxi |xi+1 − xi |, c.à.d. h est la taille maximale d’un élément du maillage.
En language clair, ce théorème exprime que pour tout v dans H 1 , l’approximation polynomiale vh converge en
O(h2 ) vers v et en O(h) pour ses dérivées. Ce thèorème indique que Vh1 est une bonne approximation polynômiale
de H 1 , et qu’on peut raisonnablement représenter la solution u ∈ H 1 d’une EDP par sa représentation de type
éléments finis uh . On qualifie couramment Vh1 comme une approximation interne de H 1 (Ω).

Remarque 10.2. Les éléments finis IP1 donnent une précision sur la solution du même ordre que les différences
finies centrées. Des rapprochements entre ces deux méthodes seront énoncées en TD.

10.2.3 D’autres éléments finis

On peut généraliser le théorème 10.1 en énoncant un resultat classique de l’approximation polynomiale, qui dit
qu’une approximation de polynômes par morceaux de degré k est précise en O(hk+1 ).
Ainsi, pour augmenter la précision des méthodes d’éléments finis, on utilise généralement des polynômes par
morceaux de degré plus élevé, par exemple l’élément IP2 qui engendre l’espace :
© ª
Vh2 = p, tel que p ∈ C 0 (0, 1), et p|[xi ,xi+1 ] ∈ IP2 pour 1 ≤ i ≤ N − 1 ,

est précis en O(h3 ).

On peut aussi utiliser des polynômes par morceaux moins réguliers. Par exemple, une approximation interne
de L2 (Ω) est l’espace Vh0
© ª
Vh0 = p, tel que p ∈ C 0 (0, 1), et p|[xi ,xi+1 ] ∈ IP0 pour 1 ≤ i ≤ N − 1 ,

qui représente les fonctions constantes par morceaux. Les éléments de la base canonique de Vh0 sont les fonctions
en escalier. Cet espace est un espace vectoriel de dimension N − 1.

10.3 Formulation variationnelle d’un problème aux limites

Nous allons considérer comme équation modèle le problème (elliptique) de Helmholtz :

λu − u00 = f, dans Ω =]0, 1[, (10.15)

avec les conditions de Dirichlet :

u(0) = u(1) = 0, (10.16)
ou les conditions de Neumann
u0 (0) = u0 (1) = 0. (10.17)
Dans la suite, nous allons considérer le problème (10.15) avec λ = 0, i.e. l’équation de Poisson. Le cas important
λ 6= 0 sera vu en TD.
78 10 Introduction à la méthode des éléments finis

10.3.1 Cas de conditions de Dirichlet

Nous considérons ici le problème (10.15),(10.16) avec λ = 0. L’équation (10.15) devient alors

−u00 = f, dans Ω =]0, 1[, (10.18)

avec u ∈ V , un espace fonctionnel qui est pour l’instant indéterminé. L’espace de travail V n’est pas une simple
coquetterie mathématique. C’est lui qui va indiquer quel approximation interne Vh utiliser2 , mais il va aussi
servir à imposer les conditions aux limites du problème.
En multipliant (10.18) par une fonction test v ∈ V et en intégrant sur tout le domaine, on obtient :
Z Z
00
− u v dx = f v dx, ∀ v ∈ V. (10.19)
Ω Ω

Le but de la formulation variationnelle est d’approcher la solution de (10.19) par des éléments de Vh1 . Or ces
éléments ne sont dérivables qu’une fois. Nous allons donc réécrire (10.19) en ne faisant intervenir que des dérivées
premières.
En utilisant la formule d’intégration par parties :
Z b Z b
£ ¤b
f 0 g dx = f g a − f g 0 dx,
a a

l’équation (10.19) devient : Z Z

£ ¤1
u0 v 0
+ u0 v 0 dx = f v dx, ∀ v ∈ V.
Ω Ω

En imposant que les fonctions tests vérifient les conditions aux limites (10.16), on obtient :
Z Z
u0 v 0 dx = f v dx, ∀ v ∈ V. (10.20)
Ω Ω

On peut maintenant déterminer l’espace V dans lequel vivent u et v. Pour que (10.20) ait un sens, on doit avoir
au moins
u0 , v 0 ∈ L2 (Ω), v ∈ L2 (Ω),
et, de plus, u et v doivent vérifier les C. L. homogènes (10.16). On en déduit aisement que :

V = H01 (Ω),

où H01 (Ω) est le sous-espace de H 1 (Ω) défini par (10.14). Pour que (10.20) ait un sens, on doit aussi avoir la
condition sur le second-membre :
f ∈ L2 (Ω).
En conclusion, si f ∈ L2 (Ω), le problème (10.18),(10.17) admet pour formulation variationnelle :





 Trouver u ∈ V = H01 (Ω) tel que
(10.21)




 a(u, v) = L(v), ∀ v ∈ V = H01 (Ω),

avec a( · , · ) la forme bilinéaire définie par :

a : V × V −→ R
Z
(u, v) 7−→ u0 v 0 dx,
Ω

et la forme linéaire

L : V −→ R
Z
v 7−→ f v dx.
Ω
2
pour briser le suspense, on utilisera Vh = Vh1 l’espace des éléments finis IP1 .
 10.3 Formulation variationnelle d’un problème aux limites 79

10.3.2 Quelques propriétés importantes des formulations variationnelles

• On rappelle que a( · , · ) est une forme bilinéaire si et seulement si pour tout v ∈ V , les applications u 7→ a(u, v)
et u 7→ a(v, u) sont des formes linéaires.
• Il y a (généralement) équivalence entre la formulation forte et faible d’un problème : on peut montrer que, au
sens des distributions :
u solution du problème fort (10.18),(10.17) ⇐⇒ u solution du problème faible (10.21).
• La forme bilinéaire a( · , · ) est symétrique :
∀ u, v ∈ V, a(u, v) = a(v, u).
• La forme bilinéaire a( · , · ) est positive :
∀ v ∈ V, a(v, v) ≥ 0.
Ces deux dernières propriétés remarquables permettent de démontrer que la formulation faible (10.21) est
équivalente à un problème d’optimisation. Soit J la fonctionnelle définie par :
J : V −→ R
1
u 7−→ J(u) = a(u, u) − (f, u),
2
et considérons le problème d’optimisation

 Trouver u ∈ V tel que :
 J(u) = inf J(v).
v∈V

On peut montrer que ce problème d’optimisation est équivalent à la formulation faible (10.21)3 . En particulier,
les algorithmes issus de l’optimisation (par exemple, la méthode du gradient conjugué) sont bien adaptés aux
problèmes issus des éléments finis.

10.3.3 Cas des conditions de Neumann

D.I.Y.

3
Indication : former l’équation d’Euler du problème à minimiser (cf. section 1.4).
80 10 Introduction à la méthode des éléments finis

10.4 Discrétisation par éléments finis IP1 du problème variationnel

On va maintenant considérer la discrétisation du problème variationnel (10.21) par les éléments finis IP1 . Ce
problème est bien posé dans l’espace fonctionnel H01 (Ω) de dimension infinie. Pour obtenir une solution discrète
uh de (10.21), nous allons utiliser un espace polynômial Vh tel que :
Vh ∼
= V = H01 (Ω).
Cet espace sera construit à partir de l’espace d’éléments finis Vh1 défini dans la section 10.2 :
n o
1
V0h = v ∈ Vh1 , v(0 = v(1) = 0 . (10.22)

10.4.1 Le problème variationnel discret

L’approximation par éléments finis de (10.21) s’écrit alors :

 Trouver uh ∈ V0h 1
tel que
(10.23)
 a(u , v ) = L(v ), 1
∀ vh ∈ V0h .
h h h

D’après les définitions de a( ·, , ·, ) et de L, ce problème s’écrit aussi :

Z Z
1
Trouver uh ∈ V0h tel que : u0h vh0 dx = f vh dx, 1
∀ vh ∈ V0h .
Ω Ω
1
Afin de représenter (10.23) sous forme matricielle, nous avons besoin de déterminer la base canonique de V0h ,
qui s’écrit : © ª
ϕj (x) ∈ Vh1 , j = 2, . . . , N − 1 , (10.24)
et qui correspond à la base de Vh1 (Eq. (10.10)) à laquelle on a retiré les fonctions ϕ0 (x) et ϕN (x) de manière à
1
ce que les éléments de V0h prennent une valeur nulle aux bords.
En écrivant la solution discrète dans la base canonique
N
X −1
uh (x) = uj ϕj (x), (10.25)
j=2

1
et en choisissant les fonctions test vh comme les éléments de la base canonique de V0h , i.e.
vh (x) = ϕi (x), i = 2, . . . , N − 1, (10.26)
alors le problème (10.23) s’écrit :
 ¡ ¢

 Trouver u2 , . . . , uN −1 ∈ RN −2 tel que

NX−1
(10.27)

 a(ϕi , ϕj ) uj = L(ϕi ), ∀ i = 2, . . . , N − 1,

j=2

en utilisant la bilinéarité de a( · , · ). Ce problème prend la forme matricielle

KU = F, (10.28)
avec ¡ ¢T
U = u2 , . . . , uN −1 ,
le vecteur solution (les C.L. ont servi pour déterminer u1 = uN = 0). Le second membre F est souvent appelé
vecteur charge, et il a pour expression :
¡ ¢ ³Z Z ´T
F = L(ϕj ) j=2,...,N −2 = f ϕ2 dx , . . . , f ϕN −1 dx .
Ω Ω

La matrice K de taille N − 2 × N − 2 est appelée matrice de rigidité. Elle correspond à la discrétisation par
éléments finis de l’opérateur −d2 /dx2 avec conditions aux limites de Dirichlet ; ses composantes sont :
Z
Kij = a(ϕi , ϕj ) = ϕ0i ϕ0j dx. (10.29)
Ω

On vient de montrer que le problème variationnel discret (10.23) peut se formuler en un système linéaire. Puisque
la matrice de rigidité K est symétrique, définie positive (c’est une des conséquences des propriétés énoncées dans
la section 10.3.2), ce problème possède une unique solution et la méthode du gradient conjugué est très bien
adapté à sa résolution.
 10.4 Discrétisation par éléments finis IP1 du problème variationnel 81

10.4.2 Explicitation du système discret pour un maillage uniforme

La détermination des composantes de la matrice de rigidité K et du vecteur charge F nécessitent l’évaluation
des intégrales où apparaissent les fonctions de base ϕj et le second membre f .
Pour cet exemple 1D, il est facile de calculer de manière exacte les éléments de K, en remarquant que, d’après
la figure 10.7, pour tout i = 2, . . . , N − 1,
Z
Kij = ϕ0i ϕ0j dx 6= 0 pour j = i − 1, i, i + 1 uniquement.
Ω

La matrice K est donc tridiagonale symétrique. Dans le cas d’un maillage uniforme, puisque les ϕj sont identiques
à une translation près, il suffit de calculer Ki,i−1 et Kii , puis de poser par arguments de symétrie et de translation
Ki,i+1 = Ki,i−1 ,
pour déterminer K complètement.
En posant
b−a 1
h= = ,
N −1 xj − xj−1
les fonctions de base s’écrivent, d’après (10.11) :




 x − xj−1
 pour x ∈ [xj−1 , xj ],
ϕj (x) = h



 xj+1 − x
 pour x ∈ [xj , xj+1 ],
h

et leur dérivée




 1
 pour x ∈ [xj−1 , xj ],
0
ϕj (x) = h

 1


− pour x ∈ [xj , xj+1 ].
h
Pour j = i − 1, on a Z
Ki,i−1 = ϕ0i ϕ0i−1 dx,
Ω
qui s’écrit, en utilisant le fait que le produit ϕ0i ϕ0i−1 a support dans [xi−1 , xi ] uniquement :
Z xi Z xi
1 1 1
Ki,i−1 = ϕ0i ϕ0i−1 dx = × dx = .
xi−1 xi−1 h h h
De même, Z Z Z
xi xi+1
2 1 1 2
Kii = (ϕ0i ) dx = dx + dx = .
Ω xi−1 h2 xi h 2 h
On en déduit que pour des conditions de Dirichlet sur maillage uniforme, la matrice de rigidité K, de taille
N − 2 × N − 2, prend la forme :  
2 −1
 
 −1 2 −1 
1  
K=   . . . . . . ,
 (10.30)
h . . . 
 −1 2 −1 
−1 2
qu’on a déjà rencontré en différences finies (section 3.1).
Remarque 10.3 (Conditions de Neumann). Dans le cas des conditions de Neumann (section 10.3.3), la matrice
de rigidité est de taille N × N , et elle s’écrit :
 
1 −1
 
 −1 2 −1 
1 
K=  .. .. .. . (10.31)
h  . . . 

 −1 2 −1 
−1 1
82 10 Introduction à la méthode des éléments finis

Fig. 10.9. Intégration d’une fonction f par la formule du trapèze.

Le calcul du vecteur charge F fait intervenir le terme analytique f (x), ce qui rend son intégration exacte plus
compliquée que précédemment. Pour cela, dans les codes d’éléments finis, on utilise couramment des formules
de quadrature numérique (trapèze, Simpson, Gauss, . . . , etc). L’une des plus utilisée est la formule du trapèze.
Comme représenté sur la figure 10.9, l’intégrale de f sur [a, b] est écrite comme la somme de son aire sur chaque
élément du maillage, i.e.
Z b N
X −1 Z xi+1
f (x) dx = f (x) dx,
a k=1 xi

et chaque aire élémentaire est évaluée par la formule de quadrature du trapèze :

Z xi+1
xi+1 − xi ¡ ¢
f (x) dx ∼
= Ixxii+1 (f ) = f (xi+1 + f (xi ) .
xi 2

D’après Quarteroni [9], l’erreur de quadrature de cette approximation s’écrit

¯Z ¯
¯ b ¯ b−a
¯ ¯
E≡¯ f (x) dx − Iab (f ) ¯ = h2 f 00 (ξ), avec ξ ∈ [a, b],
¯ a ¯ 12

où, comme précédemment, h = maxi |xi+1 − xi |. Cette formule est donc précise au second-ordre pour f (x)
quelconque. De plus, elle est exacte pour les polynômes par morceaux de degré 1. Elle peut donc être utilisée
pour calculer la matrice de rigidité K de façon exacte, puisque ϕ0i ϕ0j est un polynôme par morceaux de degré
0. Si une formule de quadrature est introduite dans une méthode d’éléments finis pour calculer des termes de
façon approchée seulement, on parle alors de crime variationnel (cf. TD).
15 La méthode des volumes finis

La méthode des volumes finies est maintenant l’une des plus répandues pour simuler numériquement les écoule-
ments de fluides. Sa principale caractéristique est d’utiliser la formulation intégrale (ou “conservative”) des
équations de bilan de la mécanique (cf. Section 4.1). C’est une différence conceptuelle majeure avec la méthode
des différences finies, qui utilise la forme différentielle de ces équations).
Dans la méthode des volumes finis, le domaine physique est subdivisé en un nombre fini de volumes de
contrôle (CVs), et les lois de conservation du problème physique sont exprimées à l’intérieur de chaque CV : on
obtient alors une équation algébrique par volume de contrôle. Ce principe de discrétisation peut s’appliquer à
des CVs de forme arbitraire, et donc à des géométries complexes. C’est une seconde différence majeure avec les
différences finies, qui ne sont généralement appliquables qu’à des géométries rectangulaires.
Malgré les différences évoquées ci-dessus, ces deux discrétisations ont de nombreux points communs. Pour
des applications simples (maillages rectangulaires uniformes, coefficients constants, etc. . . ) les discrétisation par
volumes finis centrés en sommets (cf. Sect. 15.2.2) et par différences finies sont rigoureusement identiques. Des
divergences apparaissent lorsque, par exemple, les coefficients de l’équation à résoudre sont non constants ou
discontinus : il est alors recommandé d’utiliser une formulation volume finis.
Les notions de troncature, stabilité, dissipation artificielle, etc . . . , d’un schéma numérique restent valables
pour les volumes finis : les connaissances de base en différences finies se révéleront donc précieuses pour la suite
des événements.
Puisque les équations de bilan de la mécanique ont une forme semblable à l’équation de convection-diffusion
(cf. Chap. 4), la semi-discrétisation (i.e. discrétisation spatiale) par volumes finis sera présentée pour cette
dernière équation. L’intégration en temps sera abordé ensuite : elle est basée, comme pour les différences finies,
sur des schémas pour équations différentielles ordinaires (EDOs). Cette partie reprendra de manière unifiée la
plupart des schémas temporels présentés dans la partie II de ces notes.

15.1 Principes de la discrétisation : cas de l’équation de convection-diffusion

Le principe de la discrétisation par volumes finis est présentée pour l’équation de convection-diffusion (4.11)
avec ρ = 1, soit :
∂φ
+ ∇ · (φv) − ∇ · (ν∇φ) = q. (15.1)
∂t
Dans cette équation, le coefficient de diffusion ν ≥ 0, le terme source q et la vitesse de convection v = (u, v)
sont donnés. Cette équation sera résolue dans le domaine Ω ⊂ R2 , avec des conditions de Dirichlet

φ(t, x) = f (x), ∀x ∈ ΓD , (15.2a)

et/ou de Neumann

∂φ
(t, x) = g(x), ∀x ∈ ΓN . (15.2b)
∂n
Le vecteur n désigne (comme toujours) la normale extérieure, et les frontières ΓD et ΓN sont des parties disjointes
(ΓD ∩ ΓN = ∅) de la paroi Γ = ∂Ω du domaine. La condition initiale est :

φ(t = 0, x) = φ0 (x). (15.3)

La discrétisation de ce problème repose sur la forme conservative de (15.1), soit :

Z Z ³ ´ Z
∂φ
dV + F c − F d · n dΓ = q dV, (15.4)
Ω ∂t Γ Ω

où les vecteurs F c et F d désignent respectivement les flux convectif et diffusif, soit :

F c ≡ φv, F d ≡ ν∇φ. (15.5)

106 15 La méthode des volumes finis

Le point de départ de la discrétisation par volumes finis consiste à subdiviser Ω en un nombre fini de volumes
de contrôle (ou cellules) notés Ωl , tels que
N
Ω = ∪ Ωl , (15.6a)
l=1

et qui sont supposés disjoints :

Ωl ∩ Ωm = ∅ , l 6= m. (15.6b)

Cette subdivision est appelée le maillage de Ω, et on note Γl la frontière de chaque cellule élémentaire Ωl . La
forme conservative (15.4) s’écrit alors dans la cellule Ωl comme :
Z Z ³ ´ Z
d c d
φ dV + F − F · n dΓ = q dV. (15.7)
dt Ωl Γl Ωl

par application du théorème de Reynolds 4.1 pour un volume fixe en temps. © Afin d’obtenirª un système
d’équations algébriques, il est nécessaire de définir les inconnues discrètes φl , l = 1, . . . , N . Ces inconnues
sont couramment définies comme la valeur moyenne de φ dans la cellule Ωl , soit
Z
1
φl ≡ φ dV, (15.8)
Vl Ωl

où Vl désigne le volume de la cellule Ωl . De même, on définit le second-membre de (15.7) comme

Z
1
ql ≡ q dV. (15.9)
Vl Ω l

On montrera dans les paragraphes suivants qu’il est possible d’exprimer les intégrales de flux convectif ou diffusif
de (15.7) en fonction des inconnues φl : soit, formellement
Z
F · n dΓ ∼
= Vl Fh (φl ), (15.10)
Γl

où Fh représente une approximation discrète du flux F .

L’approximation par volumes finis de l’équation de convection-diffusion s’écrit alors :
d
φl + Fhc (φl ) − Fhd (φl ) = ql , l = 1, . . . , N, (15.11)
dt
ce qui correspond à une équation différentielle algébrique. Ce système peut alors être intégré en temps avec tout
schéma pour EDO convenable.
Les paragraphes suivants vont décrire pas à pas les étapes permettant d’arriver au système discret (15.11),
qui comprennent :
1. La définition du maillage de volumes finis,
2. L’approximation des flux convectifs et diffusifs au travers des parois d’une cellule,
3. L’approximation des flux au travers des parois du domaine : i.e. l’imposition des conditions aux limites.

15.2 Les maillages de volumes finis

Comme on l’a mentionné dans l’introduction, les volumes de contrôle peuvent avoir une forme arbitraire, et
faire donc partie d’un maillage non-structuré. Cependant, dans la suite du cours, on ne considérera que des
maillages quadrangulaires.

15.2.1 Maillages centrés en cellules (“cell centered”)

Afin de discrétiser le domaine Ω, on choisit en premier lieu une subdivision de Ω par cellules rectangulaires
¤ £ ¤ £
Ωij = xi− 12 , xi+ 12 × yj− 12 , yj+ 21 ,

pour i = 1, . . . , Nx , j = 1, . . . , Ny , qui respecte les propriétés (15.6).

On pose ∆xi = xi+ 12 − xi− 12 et ∆yj = yj+ 12 − yj− 12 : la cellules est Ωij est donc de taille Vij = ∆xi × ∆yj .
Ce maillage est de type structuré, formé de quadrangles de taille non-uniforme.
 15.2 Les maillages de volumes finis 107

Fig. 15.1. Maillage non-uniforme centré en cellules. Les nœuds ◦ sont placés au centre de chaque volume de contrôle
(ou cellule). Les symboles • désignent les points où sont interpolées les conditions aux limites (cf. Section 15.3.4).

Ensuite, on place le nœud xij = (xi , yj ) au centre de la cellule Ωij . L’espacement des nœuds dans chaque
direction est :
xi+1 − xi = ∆xi
+ ∆xi+1 f i,
≡ ∆x (15.12a)
2 2
∆yj ∆yj+1 f j.
yj+1 − yj = 2 + 2 ≡ ∆y (15.12b)
Une représentation de ce type de maillages est donnée dans la figure 15.1. Notons que dans un maillage centré
en cellule, il n’y a pas de nœuds sur la frontière Γ = ∂Ω du domaine de calcul.
Une autre conséquence de cette construction est que les frontières des cellules ne se situent pas à égale
distance des nœuds, sauf dans le cas particulier d’un maillage uniforme (i.e. , ∆xi = ∆xi+1 ).
Remarque 15.1. Inutile ? : En effet, si on désigne par xi+ 12 l’abscisse de la frontière séparant Ωij de Ωi+1,j , on
obtient :
∆xi+1
xi+ 12 ≡ xi + ∆x
2 = xi+1 −
i
2 .
Mais, en utilisant (15.12),
xi + xi+1
= xi+ 12 + 14 (∆xi+1 − ∆xi ) ,
2
ce qui montre que la frontière des cellules est située au milieu des nœuds uniquement dans le cas particulier
d’un maillage uniforme.
Pour ce type de maillages, il est possible de relier la valeur moyenne des variables (15.8) aux valeurs nodales : en
effet, en utilisant la formule de quadrature du point milieu pour approcher les intégrales de volume, on obtient :
Z
1
φij ≡ φ dV ∼
= φ(xij ), (15.13)
Vij Ωij
ce qui montre que l’inconnue φij est une approximation de la valeur de φ au centre de la cellule. De même, pour
le terme source, on obtient
qij ∼
= q(xij ). (15.14)

15.2.2 Maillages centrés en nœuds (“vertex centered”)

Il existe un autre type de maillage de volumes finis : les maillages centrés en nœuds. Dans cette variante, une
distribution de nœuds est créée en premier lieu, comme en différences finies. Ensuite, les volumes de contrôle
sont construits en plaçant chaque face au milieu de 2 nœuds voisins (cf. Fig. 15.2).
Ainsi, la différence majeure avec le maillage centré en cellule est qu’il y a des nœuds sur la frontière du
domaine, et que les noeuds ne sont plus situés a priori au centre des cellules. Cela entraîne des différences
mineures sur l’expression des flux discrets, qui disparaissent dans le cas de maillages uniformes, sauf près des
bords. Néanmoins, l’ordre de troncature des discrétisations sur maillages centrées en nœuds ou en cellules est
globalement similaire, même pour des maillages non-uniformes (cf. [14]).
Enfin, notons que c’est la discrétisation sur maillages centrés en nœuds qui a le plus de point commun
avec les différences finies : pour l’équation de convection-diffusion avec coefficients constants, elles sont même
rigoureusement identiques.
Dans la suite du cours, nous ne détaillerons pas cette variante des volumes finis.
108 15 La méthode des volumes finis

Fig. 15.2. Maillage non-uniforme centré en nœuds. Les volumes de contrôle sont crées en plaçant chacune de leurs
faces au milieu des nœuds. Les nœuds (◦) sont utilisés pour discrétiser les équations, tandis que les nœuds (•) servent à
implémenter les conditions aux limites.

15.3 Approximation des flux sur un maillage centré en cellules

Dans cette partie, nous allons décrire la marche à suivre pour discrétiser les flux convectifs et diffusifs dans le
cas de maillages centrés en cellules. Pour cela, il est plus aisé d’utiliser les notations géographiques introduites
dans le paragraphe 8.1. Ces notations sont résumées dans la figure 15.3, et on notera en particulier que :
• La cellule Ωij est renommée ΩP . Son volume est noté VP = ∆xP × ∆yP .
• ΩP a une frontière (ou “face”) commune avec 4 cellules voisines : ΩS , ΩW , ΩN et ΩE . Le centre de ces cellules
est indexé avec une majuscule : par exemple, xP est le centre de ΩP .
• La frontière de ΩP se décompose comme
ΓP = Γs ∪ Γw ∪ Γn ∪ Γe ,
tel que, par exemple, la frontière Γe représente la face ’est’ de ΩP , qui la sépare de ΩE .
• Le centre des faces est indicé en minuscule : par exemple, xe est le centre de Γe .
• De même, les normales extérieures à la surface de ΩP sont indexées en minuscule : par exemple, ne est la
normale extérieure sur la face ’est’.

15.3.1 Evaluation des intégrales surfaciques

Soit F = (Fx , Fy ) une fonction de flux convectif ou diffusif. On veut déterminer une approximation du flux F · n
au travers de la cellule élémentaire ΩP , soit :
Z ·Z Z Z Z ¸
1 1
F · n dΓ = F · ne dy + F · nw dy + F · nn dx + F · ns dx (15.15a)
VP ΓP VP Γ e Γw Γn Γs
·Z Z Z Z ¸
1
= Fx dy − Fx dy + Fy dx − Fy dx . (15.15b)
∆xP ∆y P Γe Γw Γn Γs

Pour approcher chacune de ces intégrales, on utilise la formule du point milieu. Par exemple, on obtient sur la
face Γe : Z
Fx dy ∼ = ∆y P (Fx )e ,
Γe
et on déduit alors l’approximation suivante du flux :
Z
1 1 h i 1 h i
F · n dΓ ∼
= Fh ≡ (Fx )e − (Fx )w + (Fy )n − (Fy )s . (15.16)
VP Γ P ∆xP ∆y P
La fonction Fh est appelée fonction de flux discret. L’approximation du flux en fonction de ses valeurs nodales
n’est pas encore achevée, car on vient simplement d’exprimer Fh en fonction des nœuds s, w, n et e de la surface
de ΩP .
Dans une seconde étape, nous allons interpoler ces valeurs surfaciques en fonction des valeurs nodales
S, W, N, E et P. Pour cela, il convient de distinguer le cas des flux convectifs et diffusifs.
 15.3 Approximation des flux sur un maillage centré en cellules 109

Fig. 15.3. Notations géographiques pour un maillage de quadrangles centré en cellules.

15.3.2 Interpolation des flux convectifs

On se place ici dans le cas où F = F c , soit d’après (15.5) :

Fx = uφ, Fy = vφ. (15.17)

Le flux convectif discret s’écrit, d’après (15.16) :

1 h i 1 h i
Fhc ≡ (uφ)e − (uφ)w + (vφ)n − (vφ)s . (15.18)
∆xP ∆y P
Approximation centrée

Le principe de cette approximation sera détaillée pour la face ’e’, qui sépare les nœuds P et E (Fig. 15.4). La
manière la plus simple (et la plus utilisée) est d’interpoler φe comme la moyenne algébrique :

φe ∼
= 21 φE + 12 φP . (15.19)

Le flux sur la face ’e’ s’écrit alors :

(uφ)e = ue φe (15.20a)
∼
= 21 ue φE + 12 ue φP . (15.20b)

Remarque 15.2. On rappelle que la vitesse de convection v = (u, v) est une donnée du problème. On a considéré
dans l’approximation ci-dessus qu’elle était connue sur les faces des volumes de contrôle, et donc au point xe .
Il n’y a donc pas besoin de l’interpoler en fonction de ses valeurs prises au centre des cellules. On retrouvera
ce type de raisonnement dans la discrétisation des équations de Navier-Stokes sur maillages décalés (Chap. 18),
où les vitesses normales sont stockées sur les faces des cellules.
De même pour la face ’w’, on obtient :

(uφ)w ∼
= 12 uw φP + 12 uw φW . (15.21)

Pour les faces Γn et Γs , il est possible d’exprimer d’une manière analogue (vφ)n et (vφ)s en fonction des valeurs
de φ aux nœuds N, S et P. En reportant ces diverses expressions dans (15.18), on obtient l’expression finale du
flux convectif Fhc en fonction des valeurs nodales. Ce calcul est laissé au lecteur.
Il est intéressant d’observer que dans le cas d’un maillage uniforme,

∆xW = ∆xP = ∆xE = ∆x, ∆y S = ∆y P = ∆y N = ∆y, (15.22)

et dans le cas d’une vitesse de convection constante au travers de la cellule (i.e. ue = uw = u, vn = vs = v),
l’expression du flux convectif se simplifie en
110 15 La méthode des volumes finis

Fig. 15.4. Volume élémentaire ΩP pour un maillage non-uniforme centré en cellules.

u h i v h i
Fhc = φE − φW + φN − φS , (15.23)
2∆x 2∆y

et on retrouve la formule centrée des différences finies. L’interpolation (15.19) conduit donc à une approximation
centrée du flux convectif1 .

Remarque 15.3. L’ordre de convergence d’un schéma de volumes finis utilisant des interpolations du type (15.19)
pour les flux convectifs est au second-ordre, même sur maillage non-uniforme [14].
Néanmoins, l’ interpolation (15.19) n’est en toute rigueur précise au second-ordre que sur des maillages
uniformes. Notez en effet que la face Γe n’est pas située à égale distance des nœuds P et E sur un maillage
non-uniforme : son abscisse xe ∈ [xP , xE ] prend la valeur

xe = λe xE + (1 − λe ) xP , (15.24a)

avec

xe − xP ∆xP f P ≡ ∆xP + ∆xE .

λe ≡ = , ∆x (15.24b)
xE − xP fP
2∆x 2 2

Il est alors possible d’utiliser cette formule pour améliorer l’erreur de troncature locale de l’interpolation des
flux. En effet, le flux sur la face ’e’ s’écrit alors :

(uφ)e = ue φ(λe xE + (1 − λe ) xP , yP ), (15.25)

et en supposant que φ varie linéairement dans [xP , xE ], on obtient :

(uφ)e ∼
= ue [ λe φ(xE , yP ) + (1 − λe )φ(xP , yP ) ] (15.26a)
= λe ue φE + (1 − λe )ue φP . (15.26b)

De meme pour la face Γw , l’abscisse xw s’écrit :

1
Aussi notée CDS : “Centered Differencing Scheme”.
 15.3 Approximation des flux sur un maillage centré en cellules 111

xw − xW ∆xW f W ≡ ∆xW + ∆xP .

xw = λw xP + (1 − λw ) xW , λw ≡ = , ∆x (15.27)
xP − xW fW
2∆x 2 2

L’interpolation linéaire entre les nœuds W et P conduit à :

(uφ)w ∼
= λw uw φP + (1 − λw )uw φW . (15.28)

Ces approximations se réduisent à (15.20) et (15.21) dans le cas uniforme, et conduisent à une approximation
centrée dont l’erreur de troncature reste au second-ordre sur tout type de maillage. Certains auteurs (cf. [14])
déconseillent néanmoins l’utilisation de ces formules, en raison de la rupture de symétrie du système algébrique
ainsi construit.

Approximations décentrées

Au lieu d’interpoler φe linéairement entre les valeurs φP et φE (Eq. (15.26)), nous allons l’approximer suivant
la direction de la convection v au travers de la face 0 e0 de la cellule, tel que :
½
φP si (v · n)e ≥ 0,
φe = (15.29)
φE si (v · n)e < 0.

Des approximations similaires sont utilisées pour les trois autres faces. Comme on le constatera ci-après, cette
approximation correspond à un décentrage upwind du terme convectif2 .
Par exemple, dans le cas où (v · n)e ≥ 0, l’approximation (15.29) conduit à :

(uφ)e ∼
= u e φP .

Si v · n ≥ 0 au travers de chaque face de ΩP , le flux discret s’écrit :

1 h i 1 h i
Fhc = ue φP − uw φW + vn φP − vs φS , (15.30)
∆xP ∆yP
et si, de plus, on suppose que v est constante, on obtient :

Fhc = u δx− φij + v δy− φij , (15.31)

et on reconnait l’approximation décentrée amont (5.14) des différences finies, valable sur maillages non-
uniformes.
Remarque 15.4. Astuce de programmation : la généralisation de (15.31) aux vitesses de convection > 0 ou < 0
s’écrit :
1h i 1h i
Fhc = (u − |u|) δx+ + (u + |u|) δx− φij + (v − |v|) δy+ + (v + |v|) δy− φij .
2 2
Plus généralement, le test (15.29) conduit à l’approximation

(uφ)e ∼
= 21 (ue + |ue |)φP + 12 (ue − |ue |)φE ,

pour ue de signe quelconque, et on obtient l’approximation suivante du flux convectif :

1 h1 i
Fhc = 2 (ue − |ue |)φE + 12 [ (ue + |ue |) − (uw − |uw |) ] φP − 21 (uw + |uw |)φW
∆xP
1 h1 1 1
i
+ (v n − |v n |)φN + [ (v n + |v n |) − (v s − |vs |) ] φ P − (v s + |v s |)φ S ,
∆yP 2 2 2

qui généralise (15.30) aux vitesses de convection quelconques.

On connait les désavantages de l’approximation upwind (15.30) : elle est précise au premier ordre uniquement,
et introduit une forte dissipation numérique. De nombreuses stratégies ont été proposées pour augmenter la
précision des formules décentrées. La plus populaire est l’approximation QUICK3 , qui interpole φe en utilisant
3 nœuds décentrés suivant la direction de la convection, tel que
½ +
λ φE + µ+ φP + (1 − λ+ − µ+ )φW si (v · n)e > 0,
φe = (15.32)
λ− φP + µ− φE + (1 − λ− − µ− )φEE si (v · n)e < 0,

qui est à comparer avec (15.29). Les coefficients d’interpolation λ+,− et µ+,− sont déterminés de façon que φ
possède un profil parabolique entre les nœuds d’interpolation (cf. [4]). Cette approximation est au second-ordre,
et dissipe beaucoup moins que l’approximation upwind.
2
Aussi notée UDS : “Upwind Differencing Scheme”.
3
Acronyme de “Quadratic Upwind Interpolation for Convection Kinematics”, proposée par Leonard (1979).
112 15 La méthode des volumes finis

Sur maillage uniforme avec v · n constant et positif, le schéma QUICK s’écrit

Fhc = u [δx− ]Q φij + v [δy− ]Q φij , (15.33a)

avec, dans la direction x par exemple

3φi+1,j + 3φi,j − 7φi−1,j + φi−2,j
[δx− ]Q φij = . (15.33b)
8∆x
Résumé des courses

L’approximation volumes finis du flux convectif s’écrit :

Z Z
1 1
∇ · φv dV = φv · n dΓ ∼
= Fhc (vφij ). (15.34a)
Vij Ωij Vij Γij

Pour un maillage uniforme et une vitesse constante, l’interpolation linéaire (CDS) conduit à :

Fhc (φij v) = uδx0 φij + vδy0 φij , (15.34b)

où δx0 et δy0 sont les opérateurs de différences finies centrés définis en (5.15). Ainsi, l’interpolation CDS est
équivalente à l’approximation centrée en différences finies. Toujours pour un maillage uniforme et une vitesse
constante, l’approximation upwind (UDS) s’écrit
1h i 1h i
Fhc (φij v) = (u − |u|) δx+ + (u + |u|) δx− φij + (v − |v|) δy+ + (v + |v|) δy− φij . (15.34c)
2 2

15.3.3 Interpolation des flux diffusifs

Dans le cas où F = F d ≡ ν∇φ, le flux diffusif discret s’écrit :
µ ¶ µ ¶ µ ¶ µ ¶
1 h ∂φ ∂φ i 1 h ∂φ ∂φ i
Fhd ≡ ν − ν + ν − ν . (15.35)
∆xP ∂x e ∂x w ∆y P ∂y n ∂y s
Comme précédemment, intéressons-nous à l’interpolation sur la face Γe : en faisant de nouveau l’hypothèse de
linéarité de φ dans [xP , xE ], on obtient :
µ ¶
∂φ ∼ φE − φP φE − φP
ν = νe = νe . (15.36)
∂x e xE − xP fP
∆x
De même, on obtient pour la face Γw : µ ¶
∂φ ∼ φP − φW
ν = νw , (15.37)
∂x w fW
∆x
et on laisse au lecteur la tâche de construire des relations analogues pour les faces Γs et Γn .
Remarque 15.5. Dans les approximations (15.36) et (15.37), on a considéré le cas général où la viscosité ν
dépendait de x. Dans le cadre des équations de Navier-Stokes, cela arrive par exemple lorsqu’on utilise un
modèle de turbulence : la viscosité turbulente est alors habituellement calculée au centre des faces. Pour (15.36)
par exemple, il se trouve alors nécessaire d’effectuer une interpolation de la forme :
νe ∼
= λe νE + (1 − λe )νP ,
avec λe = 1/2 ou bien déterminé comme dans la remarque 15.3.
Remarque 15.6. Dans le cas du maillage uniforme (15.22) et d’une viscosité constante, le flux diffusif s’écrit :
ν h i ν h i
Fhd = 2
φW − 2φP + φE + 2
φS − 2φP + φN . (15.38)
∆x ∆y
qui correspond au schéma de différences finies à 5 points du paragraphe 8.1.
En résumé, l’approximation volumes finis du terme diffusif s’écrit :
Z Z
1 1
∇ · (ν∇φ) dV = ν∇φ · n dΓ ∼
= Fhd (φij ), (15.39a)
Vij Ωij Vij Γij

avec, sur maillage uniforme et viscosité constante :

£ 0 ¤
Fhd (φij ) = ν δxx 0
φij + δyy φij . (15.39b)
0 0
où δxx et δyy sont les opérateurs de différences finies centrés définis en (5.17).
 15.3 Approximation des flux sur un maillage centré en cellules 113

15.3.4 Implémentation des conditions aux limites

Lorsque la cellule ΩP se trouve à proximité des parois du domaine, il est nécessaire de tenir compte des conditions
aux limites (15.2) pour discrétiser les flux.

Fig. 15.5. Cas où la cellule ΩP est proche de la paroi ’est’ du domaine.

Considérons la configuration représentée sur la figure 15.5, où ΩP est la cellule située à la paroi ’est’ de Ω
: dans ce cas, la face Γe fait partie de la paroi, et les conditions de Dirichlet et de Neumann (15.2) s’écrivent
respectivement

φ e = fe , (15.40a)
µ ¶
∂φ
= ge , (15.40b)
∂x e

et l’interpolation du flux discret (15.16) doit être modifiée afin de tenir compte de ces conditions aux limites.
Plus précisément, les flux au travers de Γw , Γn , et Γs ne nécessitent pas de modification, et sont interpolés
comme dans les paragraphes 15.3.2 et 15.3.3. Seule l’interpolation sur la face Γe du terme
Z
1 1
Fx dy ∼
= (Fx )e ,
VP Γ e ∆xP

nécessite d’être modifiée.

Cas des flux convectifs

Dans le cas où Fx ≡ uφ, la condition de Dirichlet (15.40a) est aisément imposée en remarquant que :

(Fx )e = (uφ)e = ue φe = ue fe . (15.41)

L’imposition de la condition de Neumann (15.40b) nécessite plus de travail. Dans un premier temps, on effectue
le développement de Taylor tronqué au point xe de la paroi :
µ ¶
∼ ∆xP ∂φ
φP = φe − 2 , (15.42)
∂x e

qui est ensuite utilisé pour obtenir l’interpolation suivante du flux :

(Fx )e = ue φe ∼
= ue φP + ue ∆x
2 ge .
P
(15.43)

Cas des flux diffusifs

Dans le cas où Fx ≡ ν∂φ/∂x, l’imposition d’une condition de Dirichlet repose sur un raisonnement analogue à
(15.42), i.e. : µ ¶
∂φ ∼ φe − φP
= , (15.44)
∂x e ∆xP /2
114 15 La méthode des volumes finis

et le flux s’interpole alors comme : µ ¶

∂φ ∼ fe − φP
(Fx )e = ν = νe . (15.45)
∂x e ∆xP /2
Le cas d’une condition de Neumann est plus élémentaire :
µ ¶
∂φ
(Fx )e = ν = νe ge . (15.46)
∂x e

15.3.5 Un mot sur la précision spatiale

Second-ordre, sauf l’approximation upwind du terme convectif qui est au premier ordre.

15.4 Intégration en temps

Une fois complétée la discrétisation de l’équation de convection-diffusion et l’imposition des conditions aux
limites, on obtient une équation différentielle algébrique de la forme :
d
φij + Fhc (vφij ) − Fhd (φij ) = qij
0
, i = 1, . . . , Nx , j = 1, . . . , Ny , (15.47)
dt
0
où qij correspond à la discrétisation du second-membre q, augmentée des conditions aux limites. Ce système
s’écrit alors sous la forme matricielle :
d
Φ(t) + Lh Φ(t) = Q(t), (15.48)
dt
où Φ représente le vecteur où sont rangées les valeurs nodales φij , par exemple en utilisant l’ordre lexicographique
(cf. Section 8.3), et Q le vecteur second-membre. La matrice Lh est telle que :

Lh = Ch (v) − Kh , (15.49)

où Ch (v) et Kh sont les représentations matricielles respectives du flux convectif et diffusif. Par exemple :

[Ch (v)Φ]ij ≡ Fhc (vφij ) = uδx0 φij + vδy0 φij , (15.50)

si on utilise une approximation centrée avec vitesse de convection constante sur un maillage uniforme.
La discrétisation de la condition initiale (15.3) est simplement :

[Φ(0)]ij = φ0 (xij ), i = 1, . . . , Nx , j = 1, . . . , Ny . (15.51)

L’intégration en temps du système différentiel (15.48) peut alors être effectué par des schémas en temps pour
EDO.

15.4.1 Schémas explicites et critères de stabilité

La plus simple des intégrations en temps utilise le schéma d’Euler progressif (EP1) :

Φn+1 − Φn
+ Lh Φn = Qn , (15.52)
∆t
dont l’erreur de troncature temporelle est en O(∆t). Il est qualifié d’explicite car la solution au temps (n + 1)∆t
s’écrit comme :
Φn+1 = Φn + ∆t (Qn − Lh Φn ) ,
et ne nécessite par d’inversion de système linéaire, seulement une multiplication matrice-vecteur.

Analyse de stabilité de von Neumann

L’étude de stabilité par la méthode de von Neumann se fait d’une manière analogue au cas 1D (cf. section
6.3) : on le fait pour un schéma non forcé (Q ≡ 0), sur un maillage uniforme (xi = i∆x, yj = j∆y), sans
considérer les conditions aux limites, et en “gelant” les divers coefficients de l’équation : dans le cas de l’équation
de convection-diffusion, la viscosité ν et la vitesse de convection v sont considérées comme constantes. Ainsi,
nous nous ramenons à étudier la stabilité des schémas dans le cas où le flux convectif prend la forme (15.34b)
ou (15.34c) dans le cas centré CDS ou décentréUDS respectivement, et la forme centrée (15.39b) pour le flux
diffusif.
On considère alors une solution de la forme
 15.4 Intégration en temps 115

φnij = φbnkl exp [ i(kxi + lyj ) ] , (15.53)

que l’on porte dans le schéma EP1 de façon à obtenir le facteur d’amplification g défini par :

φbn+1 = g(k∆x, l∆y) φbn . (15.54)

Le critère de stabilité de von Neumann s’écrit alors :

|g(k∆x, l∆y)| ≤ 1, −π ≤ k∆x, l∆y ≤ π. (15.55)

L’étude de stabilité des schémas de convection multidimensionnels est rendue très fastidieuse par la technicité des
calculs. Dans la suite, nous nous contenterons d’énoncer sans démonstration le critère de stabilité des schémas
importants.

Stabilité du schéma d’Euler progressif

Proposition 15.7. Le schéma EP1 avec approximation centrée des flux convectif et diffusif est stable si les
conditions suivantes sont vérifiées :

2dx + 2dy ≤ 1, (15.56a)

Cx2 Cy2
+ ≤ 1, (15.56b)
2dx 2dy

avec les nombres de diffusion définis par :

ν∆t ν∆t
dx ≡ , dy ≡ , (15.57)
∆x2 ∆y 2

et les nombres de Courant (unidirectionels) :

|u|∆t |v|∆t
Cx ≡ , Cy ≡ . (15.58)
∆x ∆y

Pour les détails de l’analyse de von Neumann, cf. [14]. Ce critère de stabilité généralise les résultats obtenus
précédemment dans le cas unidimensionel.
Remarque 15.8 (Stabilité de l’équation de la chaleur). Dans le cas purement diffusif (Eq. (15.1) avec ν > 0, et
v = 0), à caractère parabolique, la condition (15.56b) est trivialement vérifiée. Dans le cas 1D, obtenu en posant
dy = 0 dans (15.56a), on retrouve le critère :
∆x2
∆t ≤ , (15.59)
2ν
déjà vu dans la Section 6.3.2. Dans le cas 2D avec ∆x = ∆y = h, le critère (15.56a) conduit à :

h2
∆t ≤ , (15.60)
4ν
qui est une condition deux fois plus restrictive que dans le cas 1D ! En 3D, on trouve que le critère de stabilité
est 6 fois plus restrictif qu’en 1D.

Remarque 15.9 (Stabilité de l’équation de convection). Dans le cas purement convectif (ν = 0, v 6= 0), de nature
hyperbolique, le critère (15.56b) ne peut être vérifié. On retrouve là l’instabilité caractéristique du schéma
explicite centré appliqué à l’équation de transport, vu dans la section 9.3.

Proposition 15.10. Dans le cas où l’approximation upwind (15.34c) du terme advectif est utilisée, la condition
de stabilité du schéma EP1 s’écrit :
Cx + 2dx + Cy + 2dy ≤ 1. (15.61)

Remarque 15.11 (La condition CFL). Pour la convection pure, on retrouve l’expression multidimensionnelle du
critère CFL de la section 9.4, soit :
|u|∆t |v|∆t
+ ≤ 1. (15.62)
∆x ∆y
Dans le cadre d’un calcul pratique, où la vitesse de convection et le maillage ne sont pas uniformes, la condition
CFL s’écrit :
116 15 La méthode des volumes finis

max (Cx )l + (Cy )l ≤ 1, l = 1, . . . , Nx Ny , (15.63a)

l

où interviennent les nombres de Courant de la cellule Ωl définis par :

|ul |∆t |vl |∆t
(Cx )l ≡ , (Cy )l ≡ . (15.63b)
∆xl ∆yl
La généralisation de cette condition CFL à l’équation de convection-diffusion s’écrit :
max (Cx )l + 2 (dx )l + (Cy )l + 2 (dy )l ≤ 1, l = 1, . . . , Nx Ny ,
l

d’après (15.61).
Schémas explicites du second ordre

Une précision du premier ordre se révèle généralement insuffisante pour le calcul de solutions instationnaires,
et des schémas du second ou troisième ordre sont couramment employés dans les simulations numériques. Il
est possible d’améliorer la précision temporelle du schéma EP1 en augmentant le nombre de niveaux de temps
utilisés : ce type de schéma sera appelé schéma multi-niveaux.
Une technique générale de construction de schémas multi-niveaux sera présenté dans la section 15.4.4. Un
exemple important est le schéma d’Adams-Bashforth (AB2) :
Φn+1 − Φn ¡ ¢ ¡ ¢
+ 21 3Lh Φn − Lh Φn−1 = 12 3Qn − Qn−1 , (15.64)
∆t
qui est un schéma explicite en O(∆t2 ), possédant 3 niveaux de temps : Φn−1 , Φn et Φn+1 . Il nécessite donc
du stockage supplémentaire dans la mémoire d’un ordinateur. Les conditions de stabilité de ce schéma sont très
compliquées à écrire, cf. [14]. Il faut néanmoins retenir que ses propriétés de stabilité sont similaires à celles du
schéma EP1, et il se révèle plus attractif pour le calcul d’écoulements visqueux instationnaires.
Pas de temps critique pour un schéma explicite

On défini le pas de temps critique ∆tcr comme le pas de temps maximal qui assure la stabilité pour un maillage
uniforme ∆x = ∆y = h donné. Lorsqu’on exprime ∆tcr en fonction de h, cela permet d’évaluer la diminution
du pas de temps nécessaire pour la stabilité du calcul lorsque le maillage est raffiné (i.e. h est diminué).
Examinons l’ordre de grandeur de ∆tcr pour le schéma EP1. Dans le cas de l’équation de convection pure,
la condition CFL (15.62) du schéma EP1 upwind s’écrit :
(|u| + |v|) ∆t ≤ h,
et le pas de temps critique est tel que :
∆tcr ∼ h.
Pour l’équation de diffusion avec approximation spatiale centrée, la remarque 15.8 donne ∆t ≤ h2 /4ν et donc :
∆tcr ∼ h2 ,
ce qui est une restriction beaucoup plus sévère que la condition CFL lorsqu’on diminue h.
Qu’en est-il pour l’équation de convection-diffusion ? Pour l’approximation centrée, en multipliant (15.56b)
par (2dx + 2dy ) et en utilisant (15.56a), on obtient :
µ ¶ µ ¶
2 dy 2 dx
Cx 1 + + Cy 1 + ≤ 2dx + 2dy ≤ 1,
dx dy
et dans le cas ∆x = ∆y = h, cette expression devient :
¡ ¢
2 u2 + v 2 ∆t2 ≤ 4ν∆t ≤ h2 ,
et on obtient donc
∆tcr ∼ h2 .
On trouve un résultat similaire pour l’approximation upwind d’après la condition (15.61). Ainsi, la condition de
stabilité du schéma EP1 pour l’équation de convection-diffusion est essentiellement dictée par le terme diffusif,
et elle est beaucoup plus restrictive que le critère CFL du cas convectif. Ce phénomène est observé par tous les
schémas explicites, et on le retrouve pour les équations de Navier-Stokes.
En résumé, l’utilisation de schémas explicites pour l’équation de convection (ν = 0) nécessite de respecter
une condition de type CFL, telle que ∆tcr ∼ h. Cette condition n’est pas trop restrictive dans la pratique, et
les schémas explicites, tels le schéma de Lax-Wendroff (cf. Section 9.6), sont très couramment employés pour
les problèmes hyperboliques.
Par contre, en présence de diffusion (ν > 0), le pas de temps critique est tel que ∆tcr ∼ h2 , ce qui peut
représenter une contrainte prohibitive en dimensions supérieures à 2, comme l’énonce la remarque 15.8.
 15.4 Intégration en temps 117

15.4.2 Schémas implicites

Afin de se libérer de la sévère contrainte de stabilité que doivent observer les schémas explicites en présence de
diffusion, il est intéressant de considérer des schémas implicites, tels le schéma d’Euler rétrograde (ER1) :

Φn+1 − Φn
+ Lh Φn+1 = Qn+1 . (15.65)
∆t
Ces schéma est en O(∆t), et l’analyse de von Neumann montre qu’il est stable sans conditions.

Schémas implicites du second ordre

Une façon simple d’obtenir un schéma implicite du second-ordre est de prendre la demi-somme des schémas
ER1 et EP1, ce qui conduit au schéma de Crank-Nicolson (CN) :

Φn+1 − Φn ¡ ¢ ¡ ¢
+ 1
2 Lh Φn+1 + Lh Φn = 1
2 Qn+1 + Qn , (15.66)
∆t
qui est stable sans condition.
D’une manière totalement différente, on peut améliorer la précision du schéma ER1 en considérant une
approximation en O(∆t2 ) de la dérivée en temps, ce qui conduit au schéma d’Euler rétrograde du second-
ordre (ER2) :
3Φn+1 − 4Φn + Φn−1
+ Lh Φn+1 = Qn+1 . (15.67)
2∆t
Ce schéma4 possède trois niveaux, et il est stable sans conditions. Voici un nouvel exemple de schéma multi-
niveaux, où l’élévation de l’ordre de précision passe par l’ajout de niveaux de temps.

Intérêts et désavantages

On le répète, le principal intérêt des schémas implicites ER1, CN et ER2 est d’être stable sans conditions. Leur
grand désavantage est de devoir résoudre au temps (n + 1)∆t un système linéaire de la forme :

( I + α∆t [Ch (v) − Kh ] ) Φn+1 = F n+1 , (15.68)

où α est un réel dépendant du schéma utilisé, et F n+1 un vecteur regroupant les termes déjà connus à l’intant
(n + 1)∆t. Par exemple, pour le schéma CN :
¡ n+1 ¢
α = 12 , F n+1 = Φn + ∆t 2 Q + Qn − Lh Φn .

Le système (15.68) peut être très délicat à inverser à cause du terme convectif : dans le cas le plus simple ou
la convection est constante et l’approximation centrée (15.34b) est utilisée, la matrice Ch (v) n’est pas définie
positive : cela entraîne des difficultés pour la résolution itérative de (15.68). Dans les application pratiques,
i.e. les équations de Navier-Stokes, le flux convectif Ch (v)Φn+1 est un terme non-linéaire : le système (15.68)
doit donc être résolu par des méthodes (quasi-)Newton. L’utilisation de schéma implicite peut donc être très
pénalisante sur le plan du temps de calcul et de la complexité des méthodes de résolution à utiliser.
Il faut donc trouver un juste milieu entre les questions de stabilité du schéma et l’efficacité du temps de
calcul. Ce compromis est représenté par les schémas semi-implicites.

15.4.3 Les schémas semi-implicites

Ces considérations motivent donc le développement de schémas qui combinent la prise en compte implicite du
terme diffusif, pour lever la restriction sur la stabilité ∆tcr ∼ h2 , et la prise en compte explicite du terme
convectif Ch (v), pour éviter la résolution de systèmes mal fichus : ce sont les schémas semi-implicites, aussi
appelés IMEX (implicites-explicites).

Construction

Il est possible de bâtir des schémas IMEX à partir de schémas implicites, en considérant l’extrapolation :

Φn+1 ∼
= 2Φn − Φn−1 , (15.69)

qui est précise en O(∆t2 ).

4
Il est aussi connu sous le nom de schéma de Gear du second-ordre dans le domaine des équations différentielles. En
CFD, il est aussi appelé BDF (“Backward Differentiation Formula”).
118 15 La méthode des volumes finis

Remarque 15.12. Si on appliquait cette linéarisation à la totalite du schéma CN, on obtiendrait le schéma
explicite AB2, Eq. (15.64). En référence à cette remarque, on désignera (15.69) comme l’extrapolation d’Adams-
Bashforth dans la suite de ces notes.
Si on applique cette extrapolation au terme convectif uniquement du schéma CN, on obtient :
Φn+1 − Φn ¡ ¢ ¡ ¢ ¡ ¢
+ 1
2 3Ch (v)Φn − Ch (v)Φn−1 − 1
2 Kh Φn+1 + Kh Φn = 1
2 Qn+1 + Qn , (15.70)
∆t
qui est le schéma IMEX à 3 niveaux connu sous le nom d’Adams-Bashforth/Cranck-Nicolson (AB/CN). Ce
schéma conduit à la résolution du système :
(I − α∆tKh ) Φn+1 = F n+1 , (15.71)
avec ¡ n+1 ¢
α = 12 , F n+1 = Φn + ∆t
2 Q + Qn + [Kh − 3Ch (v)]Φn + Ch (v)Φn−1 .
Si on compare avec le système (15.68) du schéma CN, on n’a plus qu’à résoudre une équation de Helmholtz
discrète, de type elliptique, et qui correspond donc à un système linéaire défini positif, relativement “facile” à
inverser.
Lorsqu’on applique la linéarisation (15.69) à la partie convective du schéma ER2 (15.67), on obtient le
schéma AB/ER2 :
3Φn+1 − 4Φn + Φn−1 ¡ ¢
+ 2Ch (v)Φn − Ch (v)Φn−1 − Kh Φn+1 = Qn+1 , (15.72)
2∆t
qui est précis en O(∆t2 ). Comme pour le schéma AB/CN, on se ramène à résoudre un système de Helmholtz
de la forme (15.71).

Stabilité

L’une des motivations du développement d’un schéma IMEX est d’essayer de se débarasser de la contraine
∆tcr ∼ h2 du terme diffusif en le considérant implicitement, et récupérer une contrainte du type CFL ∆tcr ∼ h
causée par le caratère explicite du terme convectif. A-t’on atteint ce but avec les schémas AB/ER2 et AB/CN ?
C’est compliqué, et les résultats précis sont fragmentaires. Dans la pratique, on essaie simplement d’utiliser
un pas de temps tel que la condition CFL (15.63) est satisfaite.
On montre néanmoins que, lorsque des approximations centrées sont utilisées :
∆tcr (AB/CN) < ∆tcr (AB/ER2) < ∆tcr (LF/CN),
avec ∆tcr (LF/CN) = CFL, c.à.d. que le pas de temps critique du schéma LF/CN (défini dans la section suivante)
est donné par la condition CFL (15.62). Ce dernier schéma est à déconseiller à cause du découplage pair/impair
de la solution sur la grille de calcul, donnant lieu à des oscillations lorsque la viscosité est faible. Le schéma
AB/ER2 est supérieur au schéma AB/CN, de part sa plus grande stabilité, et son meilleur amortissement des
fréquences élevées de la solution dans le cas non-linéaire. De plus le schéma AB/ER2 est tel que ∆tcr ∼ h, ce
qui ne semble pas être obtenu avec le schéma AB/CN.

15.4.4 Construction de schémas multi-niveaux d’ordre élevé

Technique générale

Dans cette partie, nous nous limiterons à détailler la technique de construction de schémas du second ordre au
plus. Pour cela, on considère le schéma général à trois niveaux pour l’équation de convection-diffusion (15.48) :
(1 + ²)Φn+1 − 2² Φn − (1 − ²)Φn−1 £ ¤
+ Ch γ1 Φn+1 + γ2 Φn + (1 − γ1 − γ2 )Φn−1
2∆t £ ¤
− Kh θ1 Φn+1 + θ2 Φn + (1 − θ1 − θ2 )Φn−1 (15.73)

= α1 Qn+1 + α2 Qn + (1 − α1 − α2 )Qn−1 .
où ² et les γi , θi , αi sont des paramètres à déterminer de façon que l’erreur de troncature soit au second-ordre.
L’erreur de troncature de ce schéma s’obtient à partir des développements de Taylor suivants autour du temps
n∆t :
n n
dΦn ∆t2 d2 Φ ∆t3 d3 Φ
Φn+1 = Φn + ∆t + + + O(∆t4 ), (15.74a)
dt 2 dt2 6 dt3
n n
dΦn ∆t2 d2 Φ ∆t3 d3 Φ
Φn−1 = Φn − ∆t + − + O(∆t4 ). (15.74b)
dt 2 dt2 6 dt3
 15.4 Intégration en temps 119

En premier lieu, on utilise ces développements pour obtenir l’erreur de troncature de la dérivée temporelle :
n n
(1 + ²)Φn+1 − 2² Φn − (1 − ²)Φn−1 dΦ n ²∆t d2 Φ ∆t2 d3 Φ
= + 2
+ + O(∆t3 ),
2∆t dt 2 dt 6 dt3
et de l’extrapolation :
n
dΦn ∆t2 d2 Φ
γ1 Φn+1 + γ2 Φn + γ3 Φn−1 = (γ1 + γ2 + γ3 )Φn + ∆t(γ1 − γ3 ) + (γ1 + γ3 ) 2
dt 2 dt (15.75)
n
∆t3 d3 Φ
+ (γ1 − γ3 ) 3 + O(∆t4 ).
6 dt
Cette extrapolation est utilisée pour les termes convectifs, diffusifs, et le second-membre du schéma général.
A partir de ces relations, on obtient l’erreur de troncature du schéma en ordonnant les différents termes en
puissance de ∆t :
d n
Schéma (15.73) = Φ + Ch Φn − Kh Φn − Qn
dt · ¸
d ² dΦn
+ ∆t + (2γ1 + γ2 − 1)Ch Φn − (2θ1 + θ2 − 1)Kh Φn − (2α1 + α2 − 1)Qn
dt 2 dt (15.76)
· ¸
∆t2 d2 1 dΦn n n n
+ + (1 − γ2 )Ch Φ − (1 − θ2 )Kh Φ − (1 − α2 )Q
2 dt2 3 dt
+ O(∆t3 ).

On déduit de cette identité :

• Le schéma est consistant quelle que soit l’expression de ses coefficients5 .
• Pour obtenir une précision du second-ordre, il faut annuler le terme en O(∆t) de l’erreur de troncature : les
coefficients doivent donc vérifier les 3 relations suivantes :
²
= 2γ1 + γ2 − 1 = 2θ1 + θ2 − 1 = 2α1 + α2 − 1. (15.77)
2
• Pour une précision en O(∆t3 ), les coefficients doivent vérifier les relations (15.77) ainsi que :
1
= 1 − γ2 = 1 − θ2 = 1 − α2 . (15.78)
3
Exemples classiques de schémas au second-ordre

Les relations (15.77) définissent l’expression des schémas présentés dans les sections précédentes. Il suffit de
poser γ1 = θ1 = 0 pour obtenir des schémas explicites. Le cas γ1 = 0 et θ1 6= 0 correspond aux schémas
semi-implicites.
Ces schémas peuvent être rassemblés en trois classes, suivant les niveaux utilisés par l’approximation de la
dérivée temporelle.
• Pour ² = 1, on obtient la famille des θ−schémas. Le schéma implicite CN (15.66) correspond au choix
1 1 1
γ1 = γ2 = , θ1 = θ2 = , α1 = α2 = .
2 2 2
C’est le seul schéma à deux niveaux précis en O(∆t2 ). Si, au contraire, on pose γ1 = 0, alors nécessairement
γ2 = 3/2 et on obtient le schéma semi-implicite AB/CN (15.70).
• Pour ² = 0, on obtient les schémas de type Leap-Frog. Un exemple classique est le schéma semi-implicite
LF/CN défini par :
1
γ1 = α1 = 0, γ2 = α2 = 1, θ1 = , θ2 = 0.
2
5
Parmi les schémas à deux niveaux, on reconnait les schémas du premier ordre EP1 (15.52) pour

² = 1, γ1 = θ1 = α1 = 0, γ2 = θ2 = α2 = 1.

et ER1 (15.65) pour

² = 1, γ1 = θ1 = α1 = 1, γ2 = θ2 = α2 = 0,
120 15 La méthode des volumes finis

• Pour ² = 2, on obtient les schémas d’Euler retardés du second-ordre. Le schéma totalement implicite ER2
(15.67) est donné par le jeu de coefficients :

γ1 = θ1 = α1 = 1, γ2 = θ2 = α2 = 0.

Le schéma semi-implicite AB/ER2 (15.72) correspond à utiliser plûtot :

γ1 = 0, γ2 = 2, θ1 = α1 = 1, θ2 = α2 = 0.

Enfin, il faut noter que les coefficients de ces schémas ne vérifient pas les conditions (15.78) : ils sont donc du
second-ordre en temps.

15.4.5 Compléments sur les schémas multi-niveaux

Schémas de démarrage

Prenons l’exemple du schéma à trois niveaux AB2 (15.64). La condition initiale (15.51) fournit Φ0 , et la calcul
de Φ1 avec le schéma AB2 est impossible, car on ne connait pas Φ−1 (ni Q−1 ). Ainsi, le calcul de Φ1 doit être
effectué avec un autre schéma, le “schéma de démarrage”, par exemple le schéma EP1 (15.52), soit :
¡ ¢
Φ1 = Φ0 + ∆t Q0 − Lh Φ0 .

Les calculs suivants (n ≥ 2) pourront alors être réalisés sans problèmes par le schéma AB2.
Remarque 15.13. L’utilisation d’un schéma de démarrage en O(∆t) dans le cadre d’un calcul au second-ordre ne
dégrade pas la précision globale de la méthode, qui reste en O(∆t2 ). C’est la conséquence d’un résultat classique
sur les méthodes numériques pour EDO (e.g. [6]).
Les schémas à 3 niveaux présentés ci-dessus nécessitent donc des schémas de démarrage, d’ordre O(∆t) au
moins. Par exemple, le schéma EP1 (ou ER1) pourra être choisi pour démarrer un calcul qui sera réalisé avec
les schémas ER2, AB/CN et AB/ER2.

Pas de temps variable

pas fait
 Partie I

Annexes et compléments
A Travaux Dirigés

A.1 Rappel : Méthodes numériques de base pour les equations

diérentielles ordinaires
Exercice A.1. (Résolution d'une E.D.O. non-linéaire)
On considère le problème de Cauchy :

ẏ(t) = −t y 2 (t),

t ∈ [0, 1],
y(0) = 2,
dont la solution exacte est y(t) = 2/(1 + t2 ).
1) Cette E.D.O. est elle linéaire ? Pour le schéma EP1, établir la relation de récurrence donnant
y n+1 en fonction de yn . Pour h = 1/10, calculez la solution aux temps t1 et t2 .
2) Déterminez une condition nécessaire de stabilité pour EP1. Vériez que la méthode est stable
pour le pas de temps utilisé dans la question précédente.
3) Appliquez ER1 à la résolution de cette EDO. Montrez que ce schéma se met sous la forme de
récurrence :
G(tn+1 , y n+1 ) = H(tn , y n ),
et explicitez les fonctions G et H . Pouvez-vous facilement calculer y n+1 ? Pourquoi ? Décrivez
n+1
une méthode itérative permettant de calculer y .

Exercice A.2. (Schéma de Crank-Nicolson)

On considère le problème de Cauchy :

ẏ(t) = f (t, y(t)), t ∈ [0, T ] (A.1a)

y(0) = y0 . (A.1b)

1) En appliquant la quadrature du trapèze entre les instants tn et tn+1 = tn + h, montrez qu'on

dénit le schéma :
y n+1 − y n 1
f (tn , y n ) + f (tn+1 , y n+1 ) .

= (A.2)
h 2
qui est connu sous le nom de schéma de Crank-Nicolson (C-N). Cette méthode est-elle explicite
ou implicite ?
2) En utilisant des développements de Taylor, montrez que le schéma ER1 est précis au premier
ordre, tandis que C-N est au second-ordre.
3) On applique le schéma C-N au problème test

ẏ(t) = λy(t), λ < 0.

Calculez son facteur d'amplication. Vériez que la méthode est bien au second-ordre. Etudiez
la stabilité de la solution numérique.
4) On a utilisé deux schémas X et Y pour intégrer numériquement une E.D.O. de la forme (A.1).
Les résultats sont reportés dans le tableau A.1. L'une de ces méthodes est ER1, l'autre C-N :
pouvez-vous les identier ?
5) Complétez le tableau A.2 avec les propriétés du schéma CN.
4 A Travaux Dirigés

h Schéma X Schéma Y
1 3.05 × 10−2 5.22 × 10−3
0.5 1.46 × 10−2 1.03 × 10−3
0.25 7.16 × 10−3 2.50 × 10−4
0.125 3.55 × 10−3 6.19 × 10−5
0.0625 1.77 × 10−3 1.54 × 10−5
Tableau A.1. Erreurs commises à l'instant t = t? pour diérentes valeurs du pas de temps.

Stabilité pour ẏ = λy
Schémas Type Ordre q σ(λh)
λ ∈ R−,? λ ∈ i R?

2
EP1 Expl. 1 1 + λh h < Non
−λ
1
ER1 Impl. 1 Oui Oui
1 − λh

CN

Tableau A.2. Principales propriétés des schémas numériques pour EDO.

A.2 Méthode des diérences nies

A.2.1 Diérences nies I : généralités sur les diérences nies
Exercice A.3 (Formules de diérences nies sur maillage uniforme).
On considère dans cet exercice un maillage uniforme {xi = i∆x} de pas unifome ∆x.
1) Utilisez les formules centrées du cours pour calculer l'approximation de u0 (xi ) et u00 (xi ) pour
2
les fonctions suivantes : u(x) = 1, x et x . Comparez avec les résultats analytiques. Qu'en
serait-il pour les polynômes de degré supérieur ?

2) En utilisant des développements de Taylor, montrez que la formule centrée à trois points du
cours :
0 ui+1 − 2 ui + ui−1
δxx ui = , (A.3)
∆x2
 2 
∂ u
est une approximation au second ordre de .
∂x2
i
3) Utilisez la technique générale de construction de formule de diérences nies pour établir la
formule décentrée du second ordre qui exprime (∂u/∂x)i en fonction de ui , ui+1 , et ui+2 .
Montrez que l'erreur de troncature s'écrit :

∆x2 ∂3u
 
τi = − .
3 ∂x3 i

Exercice A.4 (Formules de diérences nies sur maillage non-uniforme).

Sur le maillage non-uniforme {xi }, tel que

∆xi = xi − xi−1 , on considère la formule de diérences-
∂u
nies centrée suivante pour calculer
∂x i :
 
1 ∆xi ∆xi+1
δx0 ui = (ui+1 − ui ) + (ui − ui−1 ) (A.4)
∆xi+1 + ∆xi ∆xi+1 ∆xi
 A.2 Méthode des diérences nies 5

1) Montrez que sur maillage uniforme on retrouve la formule centrée usuelle.

2) Montrez que la formule (A.4) est une approximation au second-ordre sur maillage non uni-
forme.

Exercice A.5 (Conditionement des opérateurs de diérences nies).

On considère l'équation de Laplace :

−u00 (x) = f (x) dans Ω =]0, 1[, (A.5)

avec conditions aux limites de Dirichlet homogènes. On dénit le maillage Ωi = {xi = i∆x, i =
0, . . . , N + 1} avec le pas uniforme ∆x = 1/(N + 1).
1) Montrez que la résolution de(A.5) par méthode de diérences nies centrées conduit à la
résolution du système linéaire :
1
− KU = F (A.6)
∆x2
avec U et F respectivement les vecteur des inconnues {ui } et du second-membre {f (xi )} aux
points intérieurs du maillage, et K = TridiagN [1, −2, 1] la matrice tridiagonale de taille N × N
telle que :
 
 −2 1 
 
 
 
 
 
 1 −2 1 
 
 
 
K= .. .. .. .
 
 . . . 
 
 
 
 

 1 −2 1 

 
 
 
1 −2

2) Montrez que les valeurs propres de K s'écrivent

1 :

λk = 2 − 2 cos kπ∆x, pour k = 1, . . . , N, (A.7)

Enoncez les propriétés principales de la matrice K.

3) En remarquant que :

cos N π∆x = − cos π∆x,

déduire de la question précédente que le nombre de conditionnement de K s'écrit

1 + cos π∆x
χ(K) = ,
1 − cos π∆x

et que, lorsque le maillage devient très n (∆x → 0),

4
χ(K) ' .
π 2 ∆x2
1
On rappelle que les valeurs propres de la matrice TridiagN [a, b, c] s'écrivent :

√ kπ
λk = b − 2 ac cos , k = 1, 2, . . . , N.
N +1
6 A Travaux Dirigés

4) Que pouvez-vous déduire de ces questions sur la performance des méthodes itératives pour
résoudre le système linéaire (A.6) ?
5) De la même façon, montrez que la matrice de masse (cf. cours sur les éléments nis) D =
TridiagN [1, 4, 1]/6∆x est bien conditionnée.

A.2.2 Diérences nies II : Analyse des schémas instationnaires, consistance,

stabilité.
Exercice A.6 (Le θ-schéma pour les équations paraboliques).
Cette exercice concerne l'approximation par diérences nies de l'équation de la chaleur unidi-
mensionnelle
∂u ∂2u
− ν 2 = f (x, t), (A.9)
∂t ∂x
n
avec ν ≥ 0. On notera ∆t le pas de temps, ∆x le pas d'espace, et ui l'approximation diérences
nies de la solution u(tn , xi ) au temps tn = n∆t et au n÷ud xi = i∆x. L'approximation sera basée
sur des diérences nies centrées en espace, et l'approximation temporelle utilisera le θ -schéma
qui s'écrit, pour l'EDO ẏ(t) = f (y(t), t) :

y n+1 − y n
= θf (tn+1 , y n+1 ) + (1 − θ)f (tn , y n ), (A.10)
∆t
avec θ un paramètre réel compris dans l'intervalle [0, 1]. Les conditions aux limites sont négligées.

1) Eectuez une semi-discrétisation de l'équation de la chaleur en utilisant des diérences nies

centrées en espace. Ecrivez le systèmes d'EDO obtenu sous une forme matricielle, et de manière
0
indicielle en utilisant le symbole de dérivation discrète du cours : δxx .

2) Achevez la discrétisation en utilisant le θ-schéma. Combien de niveaux de temps possède ce

schéma ? Pour quelles valeurs de θ le schéma est il : explicite ? implicite ? A quels schémas
les valeurs θ=0 et θ=1 correspondent-elles ?

3) A l'aide de la formule de Taylor, déterminez une approximation de θφn+1 + (1 − θ)φn au

temps n∆t, jusqu'à l'ordre ∆t2 inclus.

4) Exprimez les deux premiers termes de l'erreur de troncature de (un+1 − un )/∆t et

0
δxx ui ,
2 2
approximation respective des opérateurs ∂u/∂t et ∂ u/∂x .

5) En utilisant les questions précédentes, exprimez l'erreur de troncature du schéma au point

(n∆t, i∆x), jusqu'à l'ordre ∆t2 , ∆x2 inclus.

6) Etudiez l'ordre de consistance du schéma suivant la valeur de θ.

7) Le schéma construit dans la question 2) avec θ = 1/2 est connu sous le nom de schéma de
Crank-Nicolson (CN). Ecrivez-le sous la forme itérative :

i+1
X i+1
X
1 n+1 2 n
Bij ui = Bij uj ,
j=i−1 j=i−1

où les Bl sont les matrices d'itération du schéma. Exprimez les coecients de ces matrices en
fonction du nombre de diusion r = ν∆t/∆x2 .
8) Exprimez le facteur d'amplication g(ξ) du schéma
k∆x
2 =
CN en fonction de r et de ξ = sin2
1−cos k∆x
. Rappelez la condition de stabilité de Von Neumann et analysez la stabilité du
2
schéma.

9) Question subsidiaire : même question pour θ quelconque. Montrez que si θ ≥ 1/2, le schéma
est inconditionnellement stable et, si θ < 1/2, il est stable sous la condition

∆x2
∆t ≤ .
2ν(1 − 2θ)
 A.2 Méthode des diérences nies 7

10) Que peut-on dire sur la convergence du theta-schema ?

11) Faites un bilan des propriétés du schema suivant la valeur de θ = 0, 1 ou 1/2.

Exercice A.7 (Le schéma de Richardson).

Pour discrétiser l'équation de la chaleur 1D, on considère maintenant l'approximation temporelle
du point milieu :
y n+1 − y n−1
= f (tn , y n ), (A.11)
2∆t
qui est construit en utilisant la formule de quadrature de la moyenne pour approcher l'intégrale
Z tn+1
f (t, y(t)) dt. L'approximation spatiale utilisera des diérences nies centrées.
tn−1

1) Eectuez une semi-discrétisation de l'équation de la chaleur, puis achevez la discrétisation

avec la méthode du point milieu. Le schéma résultant est connu sous le nom de schéma de
Richardson pour l'équation de la chaleur.

2) Combien de niveaux de temps possède ce schéma ? Quel est son ordre de troncature ?

3) Exprimez la matrice d'amplication du schéma. (On introduira la variable auxiliaire vin+1 = uni
pour se ramener à un schéma à 2 niveaux.)

4) Le schéma de Richardson est-il stable ?

A.2.3 Diérences nies III : schémas pour équations elliptiques

Exercice A.8 (Approximation d'un problème de convection-diusion stationnaire).

Fig. A.1. A gauche : Solution exacte de l'équation de convection diusion stationnaire pour diverses
valeur du nombre de Péclet. A droite : Solution par diérences nies centrées pour n = 13.

On considère le problème de convection-diusion avec conditions de Dirichlet :

au0 (x) − νu00 (x) = 0 dans Ω =]0, 1[, (A.1a)

u(0) = 0, u(1) = 1, (A.1b)

8 A Travaux Dirigés

avec la viscosité ν >0 et la vitesse de convection a supposée positive. La solution exacte de ce

problème aux limites est :
exp νa x − 1
u(x) = . (A.2)
exp νa − 1
1) Adimensionnez le problème (A.1) avec comme longueur de référence L, la taille du domaine,
a la vitesse de référence, et exhibez le nombre adimensionnel de Péclet noté Pe. Montrez
graphiquement que la solution exacte possède une couche limite lorsque la convection domine,
i.e. Pe  1.
2) Faites une approximation de (A.1) avec la méthode des diérences nies centrées sur la maillage
{xi = i∆x, i = 0, . . . , N + 1} avec ∆x = 1/(N + 1) le pas de maillage. Ecrivez les coecients
du schéma centré en fonction du nombre de Péclet de maille dénit par :

a∆x
Peh = . (A.3)
ν
Ecrivez le schéma complet sous sa forme matricielle et donnez son ordre de troncature.

3) Cherchez la solution discréte du schéma centré sous la forme du mode normal ui = q i . Déter-
minez q en fonction des coecients du schéma, et montrez que la solution s'écrit :

 i
2 + Peh
−1
2 − Peh
ui =  N +1 . (A.4)
2 + Peh
−1
2 − Peh

En déduire que la solution numérique présente des oscillations lorsque Peh > 2.
−3
4) Application numérique : on considère le problème (A.1) avec ν = 10 et a = 1. Calculez le
nombre de Péclet. Quel est nombre minimal de points N pour que la méthode de diérences
nies centrée ne présente pas d'oscillations ?

5) Remplacez l'approximation centrée du terme convectif de (A.1) par une approximation dé-
centrée (ou upwind ) dans le sens de la vitesse de convection. Utilisez la méthode des modes
normaux pour montrer que la solution discréte s'écrit maintenant :

i
(1 + Peh ) − 1
ui = N +1
, (A.5)
(1 + Peh ) −1

et qu'elle ne présente pas d'oscillations pour toutes valeurs de Peh .

Exercice A.9 (N÷ud fantôme et condition de Neumann)).

Pour l'équation de Helmholtz

λ u(x) − u00 (x) = f (x) dans Ω =]a, b[, (A.1)

on va mettre au point une nouvelle technique pour implémenter la condition de Neumann

u0 (a) = g− , (A.2)

dans un schéma de diérences nies centrées. Comme dans le cours, nous considèrerons le maillage
uniforme
b−a
Ωh = {xi = a + (i − 1)h, i = 1, . . . , N } avec h= ,
N −1
et noterons ui et fi les approximations respectives de u(x) et f (x) aux points de ce maillage.
 A.2 Méthode des diérences nies 9

1) En introduisant le point fantôme

x0 = a − h,
discrétisez la condition aux limites (A.2) par diérences nies centrées du second ordre.

2) Eliminez u0 du système à l'aide de la formule courante pour (A.1), écrite au bord x=a du
domaine. Comparez avec la formule obtenue avec la première méthode du cours.

3) Etudiez l'erreur de troncature du schéma, notamment au point x2 (faire le cas λ = 0 pour

alléger les calculs).

4) Au bord x = b, on imposera la condition de Dirichlet

u(b) = g+ . (A.3)

Ecrivez le système complet, de dimension N × N, sous forme tridiagonale tel qu'il soit prêt à
être résolu par l'algorithme de Thomas.

Exercice A.10 (Condition de solvabilité d'un problème de Neumann).

On considère l'équation de Poisson dans Ω ⊂ R2 ,

−∇2 u = f, (A.4)

muni des conditions de Neumann

∂u
(x) = g(x), ∀ x ∈ ∂Ω. (A.5)
∂n
1) Vériez que si ce problème admet une solution, elle est dénie à une constante près.

2) Par application du théoreme de la divergence, montrez que le problème est solvable si les
données vérient la condition de compatibilité :
Z Z
g dS + f dV = 0. (A.6)
∂Ω Ω

On considère maintenant la version 1D du problème précédant, soit :

−u00 (x) = f (x) pour x ∈ Ω =]0, 1[, (A.7a)

avec les conditions de Neumann non homogènes

u0 (0) = g− et u0 (1) = g+ . (A.7b)

Pour discrétiser ce problème par diérences nies, introduisons le maillage

Ωh = {xi = (i − 1) h, 1 ≤ i ≤ N } avec h = 1/(N − 1). (A.8)

3) Montrez que la solution analytique est dénie à une constante près et exhibez l'analogue 1D
de la condition de compatibilité (A.6).

4) Discrétisez (A.7) en utilisant la méthode du point fantôme. Ecrivez le système discret sous la
forme matricielle de taille N ×N :
A U = F, (A.9)

avec U= (u1 , . . . , uN )T . Explicitez A et F.

5) Montrez que la solution de (A.9) est dénie à une constante près. En utilisant l'alternative de
3
Fredholm , énoncez une condition de solvabilité de ce système linéaire. Comparez-la avec le
cas continu.
3
Un théorème fondamental de l'algèbre linéaire, connu sous le nom de l' alternative de Fredholm, établit
que le système linéaire
10 A Travaux Dirigés

6) Supposons que la donnée f (x) vérie la condition de compatibilité du problème continu (A.7).
On dénit les fi du second-membre de (A.9) comme :

Z h
2 2
f1 = f (x) dx,
h 0
h
Z xi + 2
1
fi = f (x) dx pour i = 2, . . . , N − 1,
h h
xi − 2
Z 1
2
fN = f (x) dx.
h h
1− 2

Montrez alors que la condition de solvabilité du système linéaire est vériée. Cette approxi-
mation des fi est-elle semblable à l'interpolation usuelle
4 ?

7) On cherche à résoudre le système linéaire (A.9) par l'algorithme de Thomas. Dans la suite de
cette question, on notera :

f1 g− fN g+
fe1 = − , feN = + .
2 h 2 h
Montrez que, à l'issue de la première étape de cet algorithme, les (N − 1) premières lignes du
système ont la forme :

 
  
h2 fe1
1 −1   u1  
 

    


  h2 fe1 + f2
    
0 1 −1   u2  
  

    
    
    · 

 · · ·  ·  
  


     i 
  2 e X 
 =  h f1 + fk  ,
 
0 1 −1   ui  
 
 
    k=2 
    
    

 · · ·  ·  
   · 

N −1
     
    X 
0 1 −1 u   2
  N −1   h fe1 + fk 

 
k=2
    
    
• • uN
 
•

et que la dernière ligne, qui s'écrit :

 N
X −1 
pivot(N ) uN = h2 fe1 + fk + feN ,
k=2

AX = B,
est solvable si et seulement si, pour tout Y tel que

AT Y = 0,

la condition de solvabilité
B T Y = 0,
est vériée. En d'autres termes, ce théorème formule des conditions de compatiblité sur le second
membre B : pour que le système linéaire possède des solutions (solution particulière + les éléments du
noyau de A), il faut que B soit orthogonal à tous les éléments du noyau de AT .
4
Conseil : utilisez la formule de quadrature du rectangle.
 A.2 Méthode des diérences nies 11

possède un pivot nul. En déduire que, si la condition de compatibilité discrète est vériée, le
système admet néanmoins une solution dénie à une constante près. Proposez une stratégie
algorithmique pour contourner la division par zéro de la dernière ligne.

A.2.4 Diérences nies IV : schémas pour équations hyperbliques

Exercice A.11 (Implémentation Matlab du θ-schéma pour l'équation de transport).
Pour discrétiser l'équation hyperbolique du transport

∂u ∂u
+a = 0, dans Ω =]0, 1[ (A.10)
∂t ∂x
avec a>0 et la condition d' inow u(ax , t) = g− (t), on considère le θ−schéma :
un+1 − uni
i
+ a θ δx0 un+1 + (1 − θ) δx0 uni = 0,
 
i pour θ ∈ [0, 1] (A.11)
∆t
sur le maillage Ωi = {xi = (i − 1)∆x, i = 1, . . . , N } avec le pas uniforme ∆x = 1/(N − 1).
1) Reconnaissez les diérents schémas du cours pour θ = 0, 1/2 et 1.
2) Pour imposer la condition au limite numérique à l'outow x = 1 on utilisera le schéma upwind
en i = N . Ecrivez-le θ−schéma sous sa forme itérative

i+1
X i+1
X
1 n+1 2 n
Bij uj = Bij uj , (A.12)
j=i−1 j=i−1

et exprimez les coecients des matrices B1 et B2 en fonction de θ et du nombre de Courant

C = a∆t/∆x.
3) Ecrivez les fonctions matlab M=cree_matrixB1(c,theta,n) et unew=theta_schema(c,theta,ug,uold,B1)
qui créent la matrice B1 et calculent la solution un+1
i en fonction des données du problème et
n
de la solution au temps précédent ui .

Exercice A.12 (Etude de la dispersion et dissipation numérique d'un schéma aux

diérences).
On considère l'équation hyperbolique du transport

∂u ∂u
+a = 0, (A.13)
∂t ∂x
avec a > 0. On a vu en cours que le schéma explicite décentré upwind

uin+1 − uni
+ a δx− uni = 0, (A.14)
∆t
correspondait à une stabilisation du schéma explicite centré par l'introduction de diusion
numérique.
Le but de cet exercice est d'introduire de façon plus rigoureuse le concept de diusion
numérique, et d'introduire celui de dispersion numérique. L'approche que nous utiliserons repose
sur l'étude de Fourier de l' équation équivalente du schéma numérique (A.14). Ces notions sont bien
sûr valables pour d'autres EDP (pas nécessairement hyperboliques) et d'autres schémas.
L'équation équivalente d'un schéma est une généralisation des notions de consistance et d'erreur
de troncature. Si on note
A u(t, x) = 0, (A.15)

l'équation à résoudre (qui correspond ici à (A.13)), et

12 A Travaux Dirigés

Ah uni = 0, (A.16)

le schéma numérique (ici, Eq. (A.14)) alors, par l'utilisation de développements de Taylor, il est
possible d'écrire que la solution exacte u(t, x) vérie aux points du maillage


Ah u(tn , xi ) = A u(tn , xi ) + Eh u(tn , xi ) , (A.17)

où Eh correspond à l'erreur de troncature du schéma. C'est le terme dominant

τin ∼ Eh u(tn , xi ) ,


pour ∆t, ∆x petit,

que vous exhibez lorsque on vous demande de calculer l'ordre de troncature du schéma.
L'équation équivalente au schéma numérique (A.16) est dénie comme l'équation aux dérivées
partielles


A u(t, x) + Eh u(t, x) = 0. (A.18)

Puisque Eh (u) comporte un nombre inni de termes qui sont obtenus par développement de Taylor,
Eq. (A.18) a peu d'utilité pratique. Elle a, par contre, un grand intérêt théorique car elle permet
d'évaluer quantitativement la précision du schéma numérique. Outre la détermination de l'ordre
de troncature du schéma, elle permet, si l'on conserve les termes d'ordre supérieur de Eh (u),
d'observer les phénomènes de dissipation et dispersion numérique du schéma, ce qui est le but du
TD.

1) Pour le schéma upwind (A.14), montrez que l'équation équivalente (A.18) s'écrit
5 :

∂u ∂u ∂2u ∂3u
+a + η 2 + ∆ 3 + O(∆t3 , ∆x3 ) = 0, (A.19)
∂t ∂x ∂x ∂x
où η et ∆ dépendent de a, ∆x et C le nombre de Courant. Pour cela, on utilisera des développe-
ments de Taylor de la solution autour de (xi = i∆x, tn = n∆t), qu'on aura pris soin d'exprimer
en fonction des derivées spatiales uniquement. Mettez en évidence le terme dominant de l'erreur
de troncature du schéma. Vériez qu'il est du premier ordre en espace et en temps.

2) Que se passe t'il pour C =1 ? Pour quelles valeurs de C 6= 1 le terme dominant de l'erreur
est le plus petit possible ?

3) Nous savons que la solution de l'équation de transport (A.13) cherchée sous la forme du mode
de Fourier
b(t) eikx ,
u(x, t) = u (A.20)

s'écrit comme
b0 eik(x−at) ,
u(x, t) = u avec u
b0 = u
b(0), (A.21)

c'est-à-dire comme une onde se propageant à la vitesse a>0 sans déformation.

D'une façon analogue, nous pouvons déterminer la représentation de Fourier de la solution
numérique, dans le but d'examiner la capacité du schéma à représenter la propagation de
cette onde. Pour cela, nous allons déterminer une solution sous la forme (A.20) de l'équation
équivalente (A.19) tronquée au second-ordre, i.e. solution de l'équation
∂u ∂u ∂2u ∂3u
+a + η 2 + ∆ 3 = 0. (A.22)
∂t ∂x ∂x ∂x
Montrez que la solution s'écrit sous la forme

0
u(x, t) = u00 (t) eik(x−a t) , avec a0 ∈ R, (A.23)

et explicitez u00 (t) en fonction de u

b0 , η et k , et a0 en fonction de a, ∆ et k.
4) Que se passe-t-il pour η>0 ? Montrez qu'on retrouve la condition CFL assurant la stabilité
du schéma upwind.

5
On dit aussi : le schema upwind (A.14) est consistant en O(∆t2 , ∆x2 ) avec l'EDP (A.19).
 A.3 Méthode des éléments nis 13

5) Tracez sur un même graphe la solution exacte (A.21) et la solution numérique (A.23) de
l'équation de transport. Explicitez les termes liés a la dissipation numérique (amortissement
de la solution) et de dispersion numérique (altération de la vitesse de propagation) du schéma
upwind. En déduire que ce schema est dissipatif au premier ordre et dispersif au second ordre.

6) Si vous avez du courage, ou de bonne connaisances en Maple, montrez que l'équation équiva-
lente du schéma de Lax-Wendro appliqué à (A.13) s'écrit :

∂u ∂u ∆x2 ∂3u ∆x3 ∂4u

+a +a (1 − C 2 ) 3 + a C (1 − C 2 ) 4 + O(∆t4 , ∆x4 ) = 0, (A.24)
∂t ∂x 6 ∂x 8 ∂x
et donc que le schéma est précis au second ordre en espace et en temps, dispersif au second
ordre et dissipatif au troisième ordre. On remarque que le coecient du premier terme dissipatif
permet de retrouver la condition de stabilité du schéma.

A.3 Méthode des éléments nis

Exercice A.13 (Fléchissement d'une poutre).
On reprend l'exemple de la poutre du cours de diérences nies, c'est-à-dire le problème de
Helmholtz

c(x) u(x) − u00 (x) = f (x) dans Ω =]0, 1[, (A.25a)

avec les conditions aux limites homogènes

u(0) = u(1) = 0. (A.25b)

1) Montrez que l'approximation par éléments nis IP1 s'écrit :

Mu + Ku = f

où K et f sont respectivement la matrice de rigidité et le vecteur charge. La matrice M,

appellée matrice de masse, a pour composante :
Z
Mij = c(x)ψi (x)ψj (x) dx.
Ω

2) Dans le cas c(x) = 1 et pour un maillage uniforme de pas h, calculez M et écrivez le système
complet.

3) Pour c(x) 6= 1 et un maillage uniforme, établir la forme du système linéaire dans le cas où
l'intégration numérique du trapèze est utilisée. Eectuez-vous un crime variationnel ? Com-
parez avec l'approximation par diérences-nies.

4) On considère maintenant c(x) = λ et les conditions mixtes

u0 (0) = 0 et u(1) = 0.

Etablir la formulation variationnelle du problème puis le système discret complet en utilisant

la formule du trapèze.

Exercice A.14 (Eléments nis pour équation d'évolution).

On considère l'approximation par éléments nis IP1 de l'équation de la chaleur 1D avec conditions
de Dirichlet homogènes.
14 A Travaux Dirigés

1) En utilisant le principe de semi-discrétisation ( cf. cours Equations paraboliques), écrivez la

version GFEM (Galerkin Finite Element Method) du schéma d'Euler progressif. Ce schéma
est-il explicite ? Qu'en est-il si on utilisait la formule de quadrature du trapèze ?

2) On peut, comme en diérences nies, eectuer une analyse de stabilité du schéma. Montrez
que le facteur d'amplication s'écrit

12 r ξ
g =1− ,
3 − 2ξ

avec ξ = (1 − cos k∆x)/2 et r = ν∆t/∆x2 . En déduire la condition sur la stabilité du schéma.

Comparez avec la stabilité obtenue par diérences nies.

Remarque A.15 (Sur les éléments nis pour les équations d'évolution).
Pour une vraie méthode GFEM (intégration exacte), il n'est guère possible d'obtenir un schéma
en temps qui soit explicite, et donc très économique en temps de calcul, tel qu'il n'y ait pas de
système linéaire de la forme
Mun+1 = F (un , un−1 , . . . ), (A.26)

à resoudre à chaque pas de temps.

Il est à remarquer que pour tout N, le nombre de conditionnement de M est

γN (M) ∼ C

où C est une constante bornée, proche de l'unité et indépendante de N. En conséquence, la

résolution de (A.26) n'est pas très coûteuse, car l'algorithme du gradient conjugué converge en un
petit nombre d'itérations indépendant de N, à la diérence de la résolution d'une équation de la
forme
Ku = F, (A.27)

comme vous l'avez vu en TP.

Néanmoins, la recherche d'une ecacité de calcul maximale des codes d'éléments nis pousse à
l'utilisation d'algorithmes explicites. C'est pour cela qu'il est courant de s'eorcer à ce que à ce que
M ∼ Id , par exemple par l'utilisation d'une sous intégration de la formulation variationnelle. On
parle alors d'éléments nis avec approximation condensée de la matrice de masse (mass lumping).
On gagne alors en temps de calcul et, parfois, en stabilité. Par contre, on perd sur le plan de la
précision de la méthode, notamment les phénomènes de dispersion et dissipation numérique sont
plus prononcés.

A.4 Méthode des volumes nis

Exercice A.16. On considère l'équation d'advection-diusion stationnaire
 
duφ d dφ
− ν = q(x), (A.28)
dx dx dx

dans le domaine Ω =]0, 1[. La viscosité ν = ν(x) est strictement positive pour tout x ∈ Ω . La
vitesse d'advection u = u(x) est de signe quelconque. On suppose néanmoins que u(x = 0) > 0 et
u(x = 1) > 0 : c'est-à-dire qu'au bord x = 0 le uide entre, et en x = 1 le uide sort. Pour que le
problème soit mathématiquement bien posé, on choisit donc les conditions aux limites suivantes :
on impose une condition d'entrée (inow) en x = 0,

φ(0) = g− , (A.29)
 A.4 Méthode des volumes nis 15

et une condition de sortie (outow) en x = 1, qui est généralement de type Neumann :

dφ
(1) = g+ . (A.30)
dx
On se propose de résoudre ce problème par la méthode des volumes nis sur un maillage non-
uniforme centré en cellules.

1) Construction du maillage. Soit {xi+ 12 , i = 0, . . . , N } une suite de réels tels que :

0 = x 21 < x 32 < · · · < xN + 12 = 1.

Construisez un maillage centré en cellules tels que les xi+ 21 représentent les frontières des
cellules élémentaires Ωi . On notera xi le centre de Ωi , et ∆xi son volume. L'espacement entre
les n÷uds xi et xi+1 sera noté ∆x
f i. Combien y-a-t'il de cellules ? Représentez ce maillage sur
un schéma. Soignez ce schéma en particulier autour d'une cellule Ωi quelconque, et près des
bords x=0 et x = 1.
2) Dénition des inconnues. On dénit les valeurs moyennes dans une cellule :
Z Z
1 1
φi = φ(x) dx, qi = q(x) dx, i = 1, . . . , N.
∆xi Ωi ∆xi Ωi

Eectuez une approximation des intégrales permettant de relier ces moyennes aux valeurs
nodales.
3) Formulation volumes nis. En intégrant l'équation (A.28) dans le volume Ωi , montrez que la
discretisation volumes nis s'écrit :

Fi+ 21 − Fi− 12 = ∆xi qi , i = 1, . . . , N, (A.31)

avec F = F c −F d , les fonctions de ux conservatif et dissipatif dont vous donnerez l'expression.
Le système discret est-il clos ? Quelle diérence notable observez-vous par rapport au cas
bidimensionnel ?
4) Interpolation des ux. Il faut maintenant exprimer les ux discrets Fi+ 12 , Fi− 12 en fonction des
valeurs nodales φi .
a) Dans le cas d'une vitesse d'advection quelconque, établissez les interpolations UDS et
CDS du ux convectif
c
Fi+ pour i = 1, . . . , N − 1, c'est-à-dire loin des bords. En déduire
1
2
l'expression de F
c . Retrouvez la formule de diérences nies upwind dans le cas u(x) =
i− 1 2
u > 0.
b) Eectuez l'interpolation du ux diusif
d
Fi+ pour i = 1, . . . , N − 1. Sur maillage uniforme
1
2
et ν(x) = ν > 0, retrouvez la formule de diérences nies :

d − Fd
Fi+ 1
2 i− 1 2 0
= νδxx φi .
∆x
5) Implémentation des conditions aux limites. Il faut maintenant tenir compte des conditions aux
limites dans l'interpolation des ux discrets près des parois.
a) Représentez les C.L. (A.29) et (A.30) sur le dessin du maillage. D'après Eq. (A.31), quels
sont les ux discrets dont l'interpolation nécessite d'être modiée a priori pour tenir
compte des conditions aux limites ?
b) Implémentez les C. L. dans les ux convectifs pour l'approximation CDS.
c) Même question pour l'approximation UDS. D'après l'énoncé, u 12 > 0 : la condition de

sortie (A.30) est-elle nécessaire ?

d) Implémentez les C. L. dans les ux diusifs.
e) On se place maintenant dans le cadre d'un problème purement diusif (u = 0 dans (A.28)),
à viscosité constante. Exprimez le schéma (A.31) pour i=1 et i=N dans le cas d'un
maillage uniforme. Comparez avec les diérences nies. Peut-on interpréter cette implé-
mentation des C. L. comme le retour du point fantôme ?
16 A Travaux Dirigés

6) Representation matricielle du schéma. Montrez que la formulation volumes nis s'écrit sous la
forme matricielle, de taille N ×N :

( Ch (u) − Dh ) Φ = Q0 .

Dans le cas CDS, donnez l'expression des matrices tridiagonales Ch (u) et Dh , ainsi que le
second-membre Q0 .
Bibliographie

En ce qui concerne l’analyse numérique “classique”, niveau Licence ou Maîtrise, citons le livre de Quarteroni [9],
et celui de Ciarlet [2] portant notamment sur l’optimisation. Citons aussi sur le même sujet le récent ouvrage
de Bergonioux [1]. On rappelle que le grand classique présentant le coté pratique de l’analyse numérique est
le “Numerical Recipes” de Press et al. [8], qui regroupe de nombreux programmes en C ou Fortran, et qui est
disponible gratuitement sur le Net3 .

Au sujet des méthodes numeriques, le bouquin de Sainsaulieu [11] est un excellent complément de ces notes
de cours, bien qu’il possède un formalisme mathématique légérement plus élevé. Pour le niveau d’une maîtrise,
citons aussi le livre de Rappaz and Picasso [10], et surtout celui d’Euvrard [3], qui possède une grande sensibilité
“mécanique”.

Pour ceux qui veulent connaître plus sur la mécanique des fluides numérique, un livre en anglais est incon-
tournable : Ferziger and Perić [4]. Citons aussi les “classiques” de Fletcher [5] et Peyret and Taylor [7], et celui
de Versteeg and Malalasekera [13] portant sur la méthode des volumes finis et son utilisation dans les codes de
calcul commerciaux tels Fluent. Pour les gens vraiment accrochés : le récent bouquin de Wesseling [14].

[1] M. Bergonioux. Optimisation et Contrôle des Systèmes Linéaires. Dunod, Paris, 2001.
[2] P. G. Ciarlet. Introduction à l’Analyse Numérique Matricielle et à l’Optimisation. Masson, Paris, 1990.
[3] D. Euvrard. Résolution Numérique des Equations aux Dérivées Partielles. Masson, Paris, 1988.
[4] J. H. Ferziger and M. Perić. Computational Methods for Fluid Dynamics. Springer, Berlin, 1996.
[5] C. A. J. Fletcher. Computational Techniques for Fluid Dynamics. Springer, Berlin, 1991.
[6] R. Peyret. Spectral Methods for Incompressible Viscous flow. Springer, New York, 2002.
[7] R. Peyret and T. D. Taylor. Computational Methods for Fluid Flow. Springer, New York, 1983.
[8] W. H. Press, B. P. Flannery, S. A. Teukolski, and W. T. Vettering. Numerical Recipes. Cambridge Univ.
Press, Cambridge, 1992.
[9] A. Quarteroni. Méthodes Numériques pour le Calcul Scientifique. Springer, Paris, 2000.
[10] J. Rappaz and M. Picasso. Introduction à l’Analyse Numérique. Presses polytechniques et universitaires
romandes, Lausanne, 1998.
[11] L. Sainsaulieu. Calcul Scientifique. Masson, Paris, 1996.
[12] G. Strang. Introduction to Applied Mathematics. Wellesley-Cambridge Press, Wellesley, MA, 1986.
[13] H. K. Versteeg and W. Malalasekera. An Introduction to Computational Fluid Dynamics. The Finite
Volume Method. Longman, Harlow, 1995.
[14] P. Wesseling. Principles of Computational Fluid Dynamics. Springer, Berlin, 2001.

3
[Link]