cours-edp

Informations

• Publié par: Ustlo
• Nombre de lectures: 22
• Langue: Français

Extrait

M´ethodes num´eriques avanc´es
pour la r´esolution
des ´equations aux d´eriv´ees partielles
[e-mail]
15 d´ecembre 2006
x
tn+1
n uX (x)
tn
1Copyright (c) 2006 Pierre Saramito
Copyleft : cette œuvre est libre, vous pouvez la redistribuer et/ou la modifier selon les termes
de la Licence Art Libre. Vous trouverez un exemplaire de cette Licence sur le site d’Art Libre
http://www.artlibre.orgainsi que sur d’autres sites.
2Chapitre 1
Discr´etisation et algorithmes
1.1 Introduction
Nous allons introduire ici quelques ´equations aux d´eriv´ees partielles tr`es classiques qui seront
par la suite ´etudi´ees en d´etail. Ces ´equations sont s´el´ectionn´ees parce qu’elles introduisent des
difficult´es sp´ecifiques, par ordrecroissant.Ces´equationssont g´en´eralementissues de probl`emesde
dphysique, et le contexte est bri`evement rapell´e. Dans la suite Ω d´esigne un ouvert born´e de R ,
avecd = 1,2,3.Eng´en´erald = 3:c’estnotreespacephysique.Pourdes raisonsde sym´etries,ilest
parfois possible de r´eduire la dimension des calculs a` d = 1 ou 2. D’autre part, T > 0 repr´esente
le temps final pour lequel nous souhaitons r´ealiser les calculs.
1.1.1 L’´equation de la chaleur
Le probl`eme s’´enonce :
trouver φ d´efini de ]0,T[×Ω dans R tel que
∂φ
ρC +u.∇φ −kΔφ = f dans ]0,T[×Ωp
∂t
φ(t=0,x) = φ (x) dans Ω0
φ(t,x) = φ (t,x) sur ]0,T[×∂ΩΓ
Dans ce syst`eme, nous avons introduit un certain nombre de notations :
φ(t,x) : la temp´erature, inconnue du probl`eme.
φ (x) : la temp´erature initiale, fonction connue.0
φ (t,x) : la temp´erature aux bords, fonction connue.Γ
u(t,x) : la vitesse du milieu, fonction connue. On peut avoir u = 0 pour un solide.
f(t,x) : la source de chaleur, fonction connue.
ρ, C , k : des constantes strictement positives, connues.p
On peut sans perte de g´en´eralit´e, supposer que ρ, C et k = 1 sont ´egales `a 1. Dans le cas ou` lap
vitesse du milieu u = 0, on obtient le probl`eme :
(C) : trouver φ d´efini de ]0,T[×Ω dans R tel que
∂φ−Δφ = f dans ]0,T[×Ω
∂t
φ(t=0,x) = φ (x) dans Ω0
φ(t,x) = φ (t,x) sur ]0,T[×∂ΩΓ
36
Le mˆeme probl`eme (C) se pr´esente dans des situations diverses :
– probl`eme de concentration des poluants en faible quantit´e dans des ´ecoulements (ex. dans les
nappes phr´eatiques) :φ repr´esente la concentration du poluant.
– probl`eme de Navier-Stokes, ou` φ = u repr´esente la vitesse elle-mˆeme. Le probl`eme est alors
non-lin´eaire.
Dans le cas ou` la vitesse du milieu u = 0 (milieu solide) et ou` on recherche l’´etat stationnaire
∂φ
= 0, obtenu a` l’´equilibre, au bout d’un temps infini, le probl`eme se r´eduit `a :
∂t
(C ) : trouver φ d´efini de Ω dans R tel que0
−Δφ = f dans Ω
φ = φ sur ∂ΩΓ
C’est le probl`eme mod`ele pour la pr´esentation des m´ethodes num´eriques.
Dans le cas ou` la vitesse du milieu u =0 (milieu fluide) et que cette vitesse devient grande, il faut
faire face `a des difficult´es :
– couche limite, ou` la solution varie tr`es brutalement
– les m´ethodes num´eriquesles plus simples sont instables : la solution approch´ee oscille fortement
alors que la solution exacte ne pr´esente pas ces oscillations.
On sait aujourd’hui rem´edier a` cela.
1.1.2 L’´equation de l’´elasticit´e
Le probl`eme s’´enonce :
dtrouver u d´efini de ]0,T[×Ω dans R tel que
2∂ u
ρ −div(λdivu.I +με(u)) = f dans ]0,T[×Ω
2∂t
u(t=0,x) = u (x) dans Ω0
∂u
(t=0,x) = g (x) dans Ω0
∂t
u = u (t,x) sur ]0,T[×∂ΩΓ
avec les notations :
u = (u ,u ,u ) :le champs de vecteur de d´eplacement1 2 3
Tε(u) = (∇u+∇u )/2 :le tenseur de d´eformation :
1 ∂u ∂ui j
ε (u) = +i,j
2 ∂x ∂xj i
λ≥0,μ> 0 : constantes positives, connues.
1.1.3 Les ´equations de Navier-Stokes
Le probl`eme s’´enonce :
trouver u et p d´efinis dans ]0,T[×Ω tels que
∂u
+u.∇u −div(νD(u))+∇p = f dans Ω
∂t
divu = 0 dans Ω
u(t=0,x) = u (x) dans Ω0
u(t,x) = u (t,x) sur ]0,T[×∂ΩΓ
avec les notations :
4u = (u ,u ,u ) :le champs de vecteur vitesse1 2 3
p :le champs de pression
TD(u) =(∇u+∇u )/2 :le tenseur des taux de d´eformation :
ν : viscosit´e, constante stritement positive, connue.
f : champs de force exterieures, connue.
1.2 Probl`emes elliptiques
On consid`ere le probl`eme mod`ele :
(P) : trouver u d´efini de Ω dans R tel que −Δu = f dans Ω
u = g sur ∂Ω
dou` Ω est un ouvert r´egulier deR ,d =1,2,3, et f et g sont des fonctions donn´ees.
1.2.1 Approximation par diff´erences finies
LorsqueΩestunsegment(d = 1),uner´egionrectangulaire(d = 2)oucubique(d =3),lam´ethode
la plus simple pour approcher la solution de (P) est la m´ethode des diff´erences finies.
La dimension un
Quitte `a changer de syst`eme de coordonn´ee, nous pouvons supposer sans perte de g´en´eralit´e que
Ω =]0,1[. L’intervalle [0,1] est d´ecompos´een une subdivision r´eguli`erede pash = 1/(N+1), avec
N ≥ 0 le nombre de points internes de la subdivision. L’id´ee de la m´ethode est de remplacer les
d´eriv´ees par des diff´erences, en s’appuyant sur un d´eveloppement limit´e de u :
2 3 4h h h0 00 (3) (4)u(x +h) =u(x )+hu (x )+ u (x )+ u (xi )+ u (x )+...i i i i i i
2 6 24
En changeanth en −h dans le d´eveloppement pr´ec´edent, il vient :
2 3 4h h h0 00 (3) (4)u(x −h) =u(x )−hu (x )+ u (x )− u (xi )+ u (x )+...i i i i i i
2 6 24
En faisant la diff´erence des deux relations pr´ec´edentes, nous obtenons :
4h2 00 (4)−h u (x +h/2)=−u(x +h)+2u(x )−u(x −h)+ u (x )+...i i i i i
12
Posons x = ih et introduisons u ≈ u(x ) l’approximation de u aux sommets de la subdivision,i i i
d´efinie par :
(P) : trouver (u ) d´efini de Ω dans R tel queh i 1≤i≤N
2 −u +2u −u = h f(x ), 1≤i≤Ni−1 i i+1 i
u = g(x )0 0
u = g(xN+1 N+1
Pour passer du probl`eme (P) au probl`eme approch´e(P) , nous avons remplac´ela d´eriv´eesecondeh
2 2 2−Δu(x ) = (−u(x +h)+2u(x )−u(x −h))/h +O(h ) parla diff´erence(−u +2u −u )/h .i i i i i+1 i i−1
5Le probl`eme approch´e s’´ecrit encore comme un syst`eme lin´eaire de N ´equations `a N inconnues :    
22 −1 0 ... 0 u h f(x )+g(x )1 1 0
2    −1 2 −1 0 u h f(x )2 2    
2    0 −1 2 −1 u h f(x )3 3=        . . ... . . .    . . . .
2−1 2 u h f(x )+g(x )N N N+1
Il s’agit d’un syst`eme lin´eaire tridiagonal : il se r´esout par une m´ethode directe (par exemple par
pivotage de Gauss) en O(N) op´erations [Sar04], ce qui est optimal.
La dimension deux ou sup´erieure
Nous supposons ici que le domaine Ω est rectangulaire. Sans perte de g´en´eralit´e, nous pouvons
2 2alors supposer que Ω =]0,1[ . Le domaine [0,1] est d´ecompos´e en une partition r´eguli`ere de pas
h = 1/(N +1), avec N ≥ 0 le nombre de points internes de la subdivision. Posonsx = (ih,jh)i,j
et introduisonsu ≈u(x ) l’approximation de u aux sommets de la subdivision, d´efinie par :i,j i,j
2−u −u +4u −u −u = h f(x ), 1≤i,j≤N i−1,j i,j−1 i,j i+1,j i,j+1 i,j
u = g(x ), i∈{0,N +1} et 0≤j ≤N +1i,j i,j
ou 1≤i≤N et j ∈{0,N +1}
Le probl`eme se met sous la forme d’un syst`eme lin´eaire dont la matrice pr´esente des ´el´ements non
nuls sur la diagonale ainsi que sur quatre parall`eles, cinq diagonales en tout. Ce type de syst`eme
2lin´eaire se r´esout par une m´ethode directe en un temps de calcul en O(N log N). L’id´ee est un
m´elange de la m´ethode pr´ec´edente et d’analyse de Fourier : voir [Sar04] pour un code c++ ou bien
la librairie fishpack pour un code fortran. Cette approche s’´etend sans difficult´e `a la dimension
trois ou plus.
1.2.2 Approximation par ´el´ements finis
Triangulation
Fig. 1.1 – Triangulation d’un cercle.
Cette m´ethode permet d’´etendre la m´ethode des diff´erences finies au cas ou` le domaine Ω n’est
dpasrectangulaire(cf Fig. 1.1). Ledomaine Ω⊂R ,d =1,2,3est d´ecompos´een une triangulation,
6not´eeT , et compos´ee de simplexes ferm´es :h
Ω =∪K∈Th
CelasupposeenparticulierqueΩsoitpolygonal,bienqueleprincipedelam´ethodepuisses’´etendre
aux domaines non polygonaux. Les simplexes sont des segments pour d = 1, des triangles pour
d = 2 et des t´etra`edres pour d = 3. On note diam(K) la plus longue arrˆete d’un simplexeK. Par
analogieavec le cas des diff´erences finies. on introduith, appel´e le pas de la triangulation, et d´efini
par :
h = max diam(K)
K∈Th
Formulation variationnelle
Afin d’introduire le probl`eme approch´e nous avons besoin d’utiliser la formulation variationnelle
du probl`eme (P). Pour cela, multiplions la premi`ere ´equation par une fonction quelconque v et
int´egrons sur Ω. Nous obtenons : Z Z
− Δuvdx = fvdx
Ω Ω
Uti