0% ont trouvé ce document utile (0 vote)
45 vues98 pages

Cours Lahm

Transféré par

Aissatou Ba
Copyright
© © All Rights Reserved
Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.
Formats disponibles
Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd
0% ont trouvé ce document utile (0 vote)
45 vues98 pages

Cours Lahm

Transféré par

Aissatou Ba
Copyright
© © All Rights Reserved
Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.
Formats disponibles
Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

Analyse théorique et numérique des équations de

Saint-Venant

Arnaud Duran
Institut Camille Jordan
Université Claude Bernard Lyon1
[Link]@[Link]
Khaled Saleh
Institut Camille Jordan
Université Claude Bernard Lyon1
[Link]@[Link]
Contents

1 Analyse théorique des lois de conservation scalaires 3


1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Etude des solutions régulières (d=1): méthode des caractéristiques . . . . . . 5
1.3 Solutions faibles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.4 Solutions entropiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.5 Résultats d’existence et d’unicité . . . . . . . . . . . . . . . . . . . . . . . . 19
1.6 Le problème de Riemann en 1d . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.6.1 La condition d’Oleinik . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.6.2 Résolution pratique des problèmes de Riemann . . . . . . . . . . . . 23

2 Schémas numériques pour les lois de conservation scalaires 29


2.1 Rappels sur les différences finies . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 Stabilité L2 dans le cas linéaire . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3 Schémas conservatifs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.4 Schéma de Godunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.5 Variantes du schéma de Godunov . . . . . . . . . . . . . . . . . . . . . . . . 43
2.5.1 Schéma de Roe-Murman. . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.5.2 Schéma de Enquist-Osher. . . . . . . . . . . . . . . . . . . . . . . . . 44
2.5.3 Solveurs de Riemann approchés . . . . . . . . . . . . . . . . . . . . . 44
2.6 Schémas monotones et entropiques . . . . . . . . . . . . . . . . . . . . . . . 47
2.6.1 Schémas monotones - schémas TVD . . . . . . . . . . . . . . . . . . . 47
2.6.2 Schémas entropiques . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

3 Problème de Riemann pour les systèmes hyperboliques 56


3.1 Systèmes hyperboliques unidimensionnels . . . . . . . . . . . . . . . . . . . . 56
3.1.1 Hyperbolicité, symétrisabilité et solutions fortes . . . . . . . . . . . . 56
3.1.2 Solutions faibles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.1.3 Inégalités d’entropie . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.2 Problème de Riemann . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.2.1 Champs caractéristiques . . . . . . . . . . . . . . . . . . . . . . . . . 63
3.2.2 Courbes intégrales - Ondes de détente . . . . . . . . . . . . . . . . . . 64
3.2.3 Discontinuités de contact . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2.4 Solution du problème de Riemann . . . . . . . . . . . . . . . . . . . . 66

1
4 Résolution pratique du problème de Riemann 67
4.1 Equations de Saint-Venant 1d . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.1.1 Caractérisation des chocs . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.1.2 Ondes de détente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.1.3 Bilan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.1.4 Un exemple : rupture de barrage . . . . . . . . . . . . . . . . . . . . 70
4.2 Le cas 2d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.2.1 Généralités . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.2.2 Champs caractéristiques . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.2.3 Problème de Riemann . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.2.4 Caractérisation des chocs et des détentes . . . . . . . . . . . . . . . . 74
4.2.5 Discontinuité de contact . . . . . . . . . . . . . . . . . . . . . . . . . 75

5 Résolution numérique des équations de Saint-Venant 76


5.1 Exemples de schémas numériques . . . . . . . . . . . . . . . . . . . . . . . . 76
5.1.1 Schéma de Godunov . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
5.1.2 Schéma de Rusanov . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.1.3 Schéma HLL (Harten - Lax - Van Leer) . . . . . . . . . . . . . . . . . 77
5.2 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.2.1 Positivité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.2.2 États stationnaires . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.2.3 Stabilité L2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

6 Complément - Aspects physiques - Dérivation des équations 87


6.1 Dérivation des équations SW . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.1.1 Conditions aux limites . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.1.2 Adimensionnement des équations . . . . . . . . . . . . . . . . . . . . 89
6.1.3 Équation de la masse . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.1.4 Equation de quantité de mouvement . . . . . . . . . . . . . . . . . . 91
6.1.5 Équation d’énergie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.2 Analyse linéaire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.2.1 Relation de dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.2.2 Ondes de gravité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.2.3 Ondes inertie-gravité . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

2
Chapter 1

Analyse théorique des lois de


conservation scalaires

1.1 Introduction
On considère des problèmes d’évolution d’inconnue upx, tq, où t P R` désigne la variable de
temps et x P Rd la variable d’espace. On a u : Rd ˆ R` Ñ K avec K Ă R un sous ensemble
convexe de R. K est communément appelé l’ensemble des états admissibles. On considère
l’EDP du premier ordre:

ÿ
d
Bt u ` divpf puqq “ 0 ô Bt u ` Bxi fi puq , t P R` , x P Rd , (1.1)
i“1

où f : R ÝÑ Rd est de classe C 8 .

Exemple 1.1.1. Obtention d’une loi de conservation. Considérons un fluide compressible


de densité ρpx, tq, se déplaçant à la vitesse upx, tq. Considérons un volume matériel Ω se
déplaçant à la vitesse du fluide, délimité par x1 ptq et x2 ptq (voir Figure 1):

Ωptq

x1 ptq x2 ptq
Figure 1: Volume matériel Ω se déplaçant à la vitesse du fluide.

La masse totale de ce volume est donnée par


ż x2 ptq
mptq “ ρpx, tqdx . (1.2)
x1 ptq

3
Remarque 1.1.1. La formule (1.2) peut être vue comme la limite de la somme de Riemann
associée à une subdivision x0 ă ¨ ¨ ¨ ă xn de l’intervalle px1 ptq, x2 ptqq. Pour n P N et t P R`
fixé, la masse totale est la somme des masses dans chaque rectangle, c’est à dire ρpxi , tq∆x.
Cette quantité converge vers mptq lorsque n tend vers `8.
dm
La conservation de la masse s’écrit “ 0. Or on a:
dt
ż x2 ptq
dm
“ Bt ρpx, tqdx ` Bt x2 ptqρpx2 ptq, tq ´ Bt x1 ptqρpx1 ptq, tq ,
dt x1 ptq

et en utilisant Bt x2 ptq “ upx2 ptq, tq et Bt x1 ptq “ upx1 ptq, tq:


ż x2 ptq ż x2 ptq
Bt ρpx, tqdx ` Bx pρuqpx, tqdx “ 0 .
x1 ptq x1 ptq

Par suite, l’intervalle rx1 ptq, x2 ptqs étant arbitraire:

Bt ρ ` Bx pρuq “ 0 . (1.3)

Un autre point de vue consiste à fixer deux points x1 , x2 (ils ne dépendent plus de t), et
d’étudier la variation de la masse dans le volume V délimité par x1 et x2 (voir Figure 2)
induite par une petite avancée en temps ∆t:

x1 ` ∆tu1 x2 ` ∆tu2

V
upx1 , tq upx2 , tq

x1 x2
Figure 2: Variation de masse dans le volume fixe V .

ż x2
En notant Mptq “ ρpx, tqdx la masse dans le volume fixe V , on a : Mpt ` ∆tq “
x1
∆tupx2 , tqρpx2 , tq . On en déduit:
∆tupx1 , tqρpx1 , tq ´ looooooooomooooooooon
Mptq ` looooooooomooooooooon
Variation de masse en x1 Variation de masse en x2
ż x2
dM Mpt ` ∆tq ´ Mptq
“ lim “ ρupx1 , tq ´ ρupx2 , tq “ ´ Bx pρuqpx, tqdx ,
dt ∆tÑ0 ∆t x1

et l’on retrouve (1.3).


Terminologie
‚ ρ est appelée variable conservative.

‚ f pρq “ ρu est appelé le flux.


4
‚ On dit que Bt u ` Bx f puq “ 0 est la forme conservative de l’équation.

‚ On dit que Bt u ` f 1 puqBx u “ 0 est la forme non conservative de l’équation.

Remarque 1.1.2.

‚ Ces équations seront dans la suite associées à un problème de Cauchy à travers la


donnée d’une condition initiale upx, 0q “ u0 pxq.

‚ Dans toute la suite, on choisira x P R ou x P Rd de sorte que nous d’aborderons pas le


problème des conditions aux limites.

Exemple 1.1.2. Autres exemples.

‚ Transport linéaire. On se donne a P R˚ , et on considère f puq “ au. L’équation (1.1)


s’écrit:
Bt u ` aBx u “ 0 ô Bt u ` Bx au “ 0 . (1.4)

‚ Equation de Burgers. La fonction de flux est définie par f puq “ u2 {2, ce qui donne lieu
aux équations suivantes:
ˆ 2˙
u
Bt u ` Bx “ 0 pforme conservativeq
2
Bt u ` uBx u “ 0 pforme non conservativeq

‚ Modèle de trafic routier. On cherche à modéliser un flot de véhicules sur une route
rectiligne sans intersections. On introduit upx, tq la densité linéique de véhicules, et
ˆ la vitesse
V puq ˙ de circulation. Un choix possible pour définir V est le suivant: V puq “
u
Vm 1 ´ , avec Vm la vitesse maximale et um une densité critique. La conservation
um
des véhicules peut alors s’écrire:

Bt u ` Bx pV puquq “ 0 pforme conservativeq


ˆ ˙
u
B t u ` Vm 1 ´ Bx u “ 0 pforme non conservativeq
2um

A noter que les formes conservative et non conservative ne sont équivalentes que pour des
solutions régulières. Or, en général, les solutions ne le sont pas.

1.2 Etude des solutions régulières (d=1): méthode des


caractéristiques
On considère le problème de Cauchy suivant:
"
Bt u ` Bx f puq “ 0 ,
(1.5)
upx, 0q “ u0 pxq ,
5
avec u0 de classe C 1 .
Définition 1.2.1: Solution régulière

Une solution régulière est une solution de classe C 1 sur R ˆ R` .

Dans toute la suite, on notera apuq “ f 1 puq.


Définition 1.2.2: Courbe caractéristique

Soit upx, tq une solution régulière de (3.2). On appelle courbe caractéristique la


fonction t ÞÝÑ xptq où xptq est solution du problème de Cauchy:
" 1
x ptq “ apupxptq, tqq , t ą 0 ,
(1.6)
xp0q “ y .

y est appelé pied de la caractéristique.

Remarquons que si u est solution régulière, on a px, tq ÞÑ apupx, tqq P CpR ˆ R` q. Le


théorème de Cauchy-Lipschitz assure l’existence d’une solution maximale unique sur un
intervalle du type r0, T ˚r.
Proposition 1.2.1

Soit x : r0, T ˚rÝÑ R une courbe caractéristique. Alors:

1. Le long de la courbe caractéristique t ÞÑ xptq, la solution u est constante.

2. La courbe t ÞÑ xptq est une droite.

Proof. On s’intéresse au premier point, et on va montrer que t ÞÑ upxptq, tq est constante.


d dxptq
upxptq, tq “ Bx upxptq, tq ` Bt upxptq, tq
dt dt
“ a pupxptq, tqq Bx upxptq, tq ` Bt upxptq, tq
“ Bx f pupxptq, tqq ` Bt upxptq, tq “ 0 .
On aborde le second point. Pour tout t P r0, T ˚r, upxptq, tq “ upxp0q, 0q “ upy, 0q “
u0 pyq. Donc l’équation de la caractéristique devient: x1 ptq “ apu0 pyqq, et donc xptq “
y ` tapu0 pyqq.
Ainsi, pour connaitre upx, tq, il faut trouver y, le pied de la caractéristique qui passe par
le point xptq, comme illustré dans la Figure 3.
Exemple 1.2.1. Transport linéaire. On considère le problème de Cauchy:
"
Bt u ` Bx pcuq “ 0 ,
upx, 0q “ u0 pxq ,
avec f puq “ cu et c P R˚ . On a apuq “ c. L’équation des caractéristiques s’écrit x “
y ` ct ô y “ x ´ ct. Pour tout px, tq P R ˆ R` , la caractéristique qui passe par px, tq
6
t
px, tq
b

1
pente
apu0 pyqq

y x

Figure 3: Caractéristique de pied y associée au point px, tq: on a upx, tq “ u0 pyq.

a pour pied y “ x ´ ct. La solution étant constante le long de la caractéristique, on a


upx, tq “ u0 pyq “ u0 px ´ ctq. La condition initiale est transportée à la vitesse c (voir Figure
4).
t

px, tq
b

y “ x ´ at x

u0 pxq u0 px ´ ctq
ct

Figure 4: Caractéristique de pied y associée au point px, tq: on a upx, tq “ u0 pyq.

Exemple 1.2.2. Equation de Burgers. On considère l’équation de Burgers avec la condition


initiale suivante:
$ ˆ 2˙

’ u

’ Bt u ` Bx “ 0,
& $2
& 1 si x ď 0 ,



’ upx, 0q “ 1 ´ x si x P r0, 1s ,
% %
0 si x ě 1 .
7
On a f puq “ u2 {2, apuq “ f 1 puq “ u. L’équation des caractéristiques est donnée par
xptq “ y ` apu0 pyqqt “ y ` u0 pyqt.

‚ Si y ď 0, alors xptq “ y ` t.

‚ Si y P r0, 1s, alors xptq “ y ` p1 ´ yqt.

‚ Si y ě 1, alors xptq “ y.

Par conséquent, les pentes des droites caractéristiques dépendent de y, et peuvent se couper
(voir Figure 5). A leur intersection, la solution régulière n’est plus définie (elle est multivar-
iée).

t
1 b

1 x

u
u0 pxq
upx, tq

1 x

Figure 5: Caractéristiques pour l’équation de Burgers.

Le problème illustré dans l’exemple précédent est générique. En réalité, dès lors que f
n’est pas linéaire, et même si la donnée initiale est aussi régulière que l’on veut, les solutions
régulières ont en général une durée d’existence finie.

8
Théorème 1.2.1
On se place dans le cas d “ 1. Soit u0 P C 1 pRq une condition initiale bornée avec u10
bornée.

i) Si a ˝ u0 est croissante, alors il existe une unique solution régulière sur R ˆ R` .

ii) Si a ˝ u0 n’est pas croissante, il existe y tel que pa ˝ u0 q1 pyq ă 0. On note alors
1
T˚ “ ´ . (1.7)
inf yPRpa ˝ u0q1 pyq

Alors il existe une unique solution classique au problème de Cauchy sur R ˆ


r0, T ˚r.

iii) Pour t ě T ˚ , il n’y a pas de solution régulière.

Proof. i) : Supposons a ˝ u0 croissante. Soit px, tq P R ˆ R` . La fonction y ÞÑ y ` a ˝ u0 pyqt


est strictement croissante sur R, et donc bijective. Il existe donc un unique y P R tel que
xptq “ y ` a ˝ u0 pyqt “ x, que l’on notera ypx, tq. En revenant à l’équation y ` a ˝ u0 pyqt “ x,
on a:

Bt y ` Bt ypa ˝ u0 q1 pyqt ` a ˝ u0 pyq “ 0 ñ Bt y r1 ` pa ˝ u0 q1 pyqts “ ´a ˝ u0 pyq


Bx y ` Bx ypa ˝ u0 q1 pyqt “ 1 ñ Bx y r1 ` pa ˝ u0 q1 pyqts “ 1

Par suite, en définissant u par upx, tq “ u0 pypx, tqq pour tout px, tq P R ˆ R` , on obtient:

Bt u ` Bx f puq “ Bt y u10pyq ` a ˝ u0 pyq Bx pu0 pyqq


(1.8)
“ Bt y u10pyq ` a ˝ u0 pyq Bxy u10 pyq ,

et en multipliant par 1 ` pa ˝ u0 q1 pyqt ą 0:

r1 ` pa ˝ u0 q1 pyqts pBt u ` Bx f puqq


´ ¯
1 1 1
“ Bloooooooooooomoooooooooooon loooooooooooomoooooooooooon u0 pyq (1.9)
t y r1 ` pa ˝ u0 q pyqts `a ˝ u0 pyq B x y r1 ` pa ˝ u 0 q pyqts
´a˝u0 pyq 1

“ 0.

Ceci garantit l’existence d’une solution globale u sur R ˆ R` .


Concernant l’unicité, si u est une solution régulière, la Proposition 1.2 garantit que upx, tq “
u0 pyq, où y est le pied de la courbe caractéristique qui passe par le point px, tq, c’est à dire
ypx, tq.

ii) : Soit t P r0, T ˚r. On considère la fonction ft : y ÞÑ y ` a ˝ u0 pyqt. On a lim ft pyq “ ˘8


yÑ˘8
car u0 est bornée. Par ailleurs, en vertu de la Définition 1.7:

f 1 pyq “ 1 ` pa ˝ u0 q1 pyqt ą 0 ,
9
et donc ft est strictement croissante. On en déduit qu’il existe un unique réel y “ ypx, tq tel
que ft pyq “ x, c’est à dire tel que xptq “ y ` a ˝ u0 pyqt “ x. Muni de ce résultat, on conclut
exactement comme dans le point précédent, à la différence près que la solution u est cette
fois définie sur R ˆ r0, T ˚r.
1
iii) : Considérons à présent t ą T ˚ . On a donc t ą ´ , et il existe y ˚ P R tel
inf yPR pa ˝ u0 q1 pyq
que
1 ` pa ˝ u0 q1 py ˚qt ă 0 .
Revenons à la fonction ft : y ÞÑ y ` a ˝ u0pyqt. Nous avons vu que si la solution est régulière,
elle s’écrit forcément upx, tq “ u0 pyq où y est le pied de la caractéristique passant par px, tq.
On a toujours lim ft pyq “ ˘8, avec cette fois ft1 py ˚ q ă 0. Ces éléments fournissent des
yÑ˘8
indications utiles sur le comportement de ft , illustré à travers la représentation graphique
proposée dans la Figure 6:

x0

y1 y˚ y2 y

Figure 6: Représentation graphique de la fonction ft .

Pour x0 bien choisi, l’équation ft pyq “ x0 admet donc trois racines distinctes, dont l’une
est égale à y ˚ . Ainsi: "
x0 “ y1 ` pa ˝ u0 qpy1 qt ,
(1.10)
x0 “ y ˚ ` pa ˝ u0 qpy ˚qt ,
ce qui implique y1 ´ y ˚ “ pa ˝ u0 qpy ˚qt ´ pa ˝ u0 qpy1 qt. Par conséquent:
y1 ‰ y ˚ ñ pa ˝ u0 qpy ˚q ‰ pa ˝ u0 qpy1 q ñ u0 py ˚q ‰ u0 py1 q ,
ce qui implique que upx0 , tq est multivaluée.
Exercice 1.2.1. Explosion en temps fini (1). Soit u “ upx, tq une solution régulière définie
sur R ˆ r0, T r avec T P R` . Montrer que
Bx upx, tq r1 ` pa ˝ u0 q1 pyqts “ u10 pyq ,
où y désigne le pied de la caractéristique passant par px, tq. En déduire que si a ˝ u0 n’est
pas croissante, alors T ď T ˚ .
On a upx, tq “ u0 pyq, où x “ y ` pa ˝ u0 qpyqt, avec y “ ypx, tq. Il en résulte:
Bx upx, tq “ Bx rx ´ pa ˝ u0 qpyqts u10 pyq .
“ r1 ´ Bx y pa ˝ u0 q1 pyq ts u10 pyq
“ u10 pyq ´ pa ˝ u0 q1 pyq t Blooomooon
1
x y u0 pyq ,

Bx pu0 pyqq

10
et donc
Bx upx, tq p1 ` pa ˝ u0 q1 pyqtq “ u10 pyq , @px, tq P R ˆ r0, T r . (1.11)
1
Supposons que a ˝ u0 n’est pas croissante. On a donc T ˚ “ ´ ă `8.
inf yPR pa ˝ u0 q1 pyq
1
Supposons T ą T ˚ . Alors il existe un réel t tel que T ˚ ď t ď T et t “ ´ avec
pa ˝ u0 q1 py0 q
y0 P R, et donc le membre de droite dans (1.11) est nul. Notons que pa ˝ u0 q1 py0 q ă 0 ñ
u10 py0 q ‰ 0, ce qui est incompatible avec (1.11).
Exercice 1.2.2. Explosion en temps fini (2). On considère une solution régulière de classe
C 2 définie sur R ˆ r0, T r avec T P R` . On pose:

vpx, tq “ pa1 ˝ uqpx, tq Bx upx, tq , @px, tq P R ˆ r0, T r ,

et pour toute courbe caractéristique Xptq de pied y, on définit wy ptq “ vpXptq, tq.
1. Montrer que wy satisfait l’équation de Ricatti:
d
z “ ´z 2 , (1.12)
dt

Rappelons que la caractéristique Xptq de pied y est définie par Xptq “ y ` pa ˝ u0 qpyqt.
On a:
d d
wy “ rpa1 ˝ uqpXptq, tq BxupXptq, tqs
dt dt
d
“ XptqBx rpa1 ˝ uq ˆ Bx us pXptq, tq ` Bt rpa1 ˝ uq ˆ Bx us pXptq, tq
dt “ ‰
pa ˝ u0 qpyq a2 ˝ u ˆ pBx uq2 ` pa1 ˝ uq ˆ Bx2 u pXptq, tq
“ loooomoooon
pa˝uqpXptq,tq
“ ‰
` a2 ˝ u ˆ pBx uBt uq ` pa1 ˝ uq ˆ Bxt 2
u pXptq, tq
“ ‰
“ pa2 ˝ uq Bx u rpa ˝ uqBx u ` Bt us ` pa1 ˝ uq pa ˝ uqBx2 u ` Bxt
2
u .

On a pa ˝ uqBx u ` Bt u “ 0 car u est solution régulière, ce qui entraîne la nullité du


membre de gauche dans la dernière égalité. Concernant l’autre membre:
2
Bxt u “ Bx pBt uq “ ´Bx ppa ˝ uqBx uq
“ ´pa ˝ uqBx2 u ´ pa1 ˝ uq pBx uq2 ,
d
et donc wy “ ´pa1 ˝ uq ˆ pa1 ˝ uq pBx uq2 “ ´pwy q2 .
dt
v0 pyq
2. Montrer que wy ptq “ , où v0 pyq “ vpy, 0q.
1 ` tv0 pyq
α
Un calcul élémentaire donne que toute fonction de la forme zptq “ avec α ą 0,
1 ` αt
t P r0, 1{αr est solution de (1.12). On conclut avec wy p0q “ vpXp0q, 0q “ vpy, 0q “
v0 pyq.
11
3. En déduire que si a ˝ u0 n’est pas croissante, alors T ď T ˚ . Comme précédemment, si
a ˝ u0 n’est pas croissante, on a T ˚ ă `8 et pour tout T ą T ˚ , il existe y0 P R tel que
1
pa ˝ u0 q1 py0 q ă 0 et T ď t “ ´ ď T . Considérons la caractéristique X0 ptq
pa ˝ u0 q1 py0 q
issue de y0 . En revenant à (1.11), on a Bx upX0p0q, 0q “ u10 py0q, et donc

vpX0 p0q, 0q “ pa1 ˝ uqpX0 p0q, 0q ˆ Bx upX0 p0q, 0q


“ pa1 ˝ u0 qpy0 q u10py0 q “ pa ˝ u0 q1 py0 q ,

de sorte tv0 py0 q “ ´1, et on conclut comme dans l’exercice précédent.

Dès lors le problème est de trouver une stratégie pour définir les solutions au delà du
temps T ˚ . A cet effet, on introduit la notion de solution faibles, c’est à dire de solution au
sens des distributions.

1.3 Solutions faibles


Revenons au problème "
Bt u ` Bx f puq “ 0 ,
(1.13)
upx, 0q “ u0 pxq .
Pour simplifier l’analyse et les notations, on se placera dans le cadre d “ 1, en gardant à
l’esprit que les prochains développements se généralisent en dimension supérieure.
Définition 1.3.1: Solution faible
Soit u0 P L8 loc pRq. On dit que u est solution faible de (1.13) sur R ˆ r0, T r si u P
8
Lloc pRˆs0, T rq et
żTż ż
puBt φ ` f puqBxφq dxdt ` u0 pxqφpx, 0qdx “ 0 , @φ P DpR ˆ r0, T rq .
0 R R

Remarque 1.3.1.

‚ Si T “ `8, on dit que la solution est globale sur R ˆ R` .

‚ Dans la définition Supppφq Ă R ˆ r0, T r et non pas Supppφq Ă Rˆs0, T r, sinon la


condition initiale ne serait pas prise en compte à travers le terme:
ż
u0 pxqφpx, 0qdx .
R

Proposition 1.3.1
Si u est régulière, alors u est solution faible si et seulement si u est solution classique.

12
Proof. Supposons u régulière. Pour tout φ P DpR ˆ r0, T rq, on a:
żTż ż
puBtφ ` f puqBx φq dxdt ` u0 pxqφpx, 0qdx
0 R R
żTż ż
“ pBt u ` Bx f puqq φdxdt ` pu0 pxq ´ upx, 0qq φpx, 0qdx .
0 R R

ð: Supposons que u soit solution classique. Alors le membre de droite dans l’égalité précé-
dente est nul pour toute fonction φ P DpR ˆ r0, T rq. Ceci implique la nullité du membre de
gauche, et donc u est solution faible.
ñ: Supposons que u soit solution faible. Cette fois, le membre de gauche est nul pour tout
φ P DpR ˆ r0, T rq. En particulier, pour tout φ P DpRˆs0, T rq, on a:
żTż
pBt u ` Bx f puqq φdxdt ,
0 R

żce qui implique Bt u ` Bx f puq “ 0 pour tout x P R et tout t ą 0. Par suite, on a


pu0 pxq ´ upx, 0qq φpx, 0qdx “ 0 pour tout φ P DpR ˆ r0, T rq, et donc:
R
ż
pu0 pxq ´ upx, 0qq ψpxqdx “ 0 , @ψ P DpRq . (1.14)
R

Par conséquent u0 pxq “ upx, 0q pour tout x P R, ce qui permet de conclure que u est solution
régulière de (1.13).
La définition précédente donne la possibilité de considérer des solutions faibles discontin-
ues. Cependant, ces solutions n’admettent pas n’importe quelle discontinuité. La théorème
qui suit précise les conditions que doit vérifier une solution faible régulière de part et d’autre
d’une courbe régulière de discontinuité.
Théorème 1.3.1
Soit u P L8 loc pRˆs0, T rq de classe C en dehors d’une courbe régulière Γ séparant
1

Rˆs0, T r en deux composantes connexes Ω´ et Ω` . Alors il s’agit d’une solution


faible si et seulement si:

‚ C’est une solution régulière dans Ω´ et Ω` .

‚ Le long de Γ, les sauts rus de u et rf puqs de f puq sont reliés par la relation de
Rankine-Hugoniot:
nt rus ` nx rf puqs “ 0 , (1.15)
où ~n “ pnx , nt q est un vecteur unitaire normal à Γ qui pointe vers la droite
pnx ą 0q.

Proof. Voir cas système

13
t
Ω´ Γ

~n
nt
Ω`
nx

Figure 7: Relation de Rankine-Hugoniot: séparation du domaine en deux composantes con-


nexes Ω´ , Ω` par une courbe de discontinuité Γ.

Remarque 1.3.2. Comme f est de classe C 1 , si u est bornée (à valeurs dans un intervalle
ra, bs par exemple), alors il existe M ą 0 tel que pour tout u P ra, bs, on a |f puq| ď M|u|.
Alors, par la relation (1.15), en tout point px, tq P Γ, on a

|nt px, tq| ď M|nx px, tq| , (1.16)

Si la courbe est définie par l’ensemble des points du plan Mpx, tq tels que Γpx, tq “ 0,
la normale au point px, tq est donnée par ~n “ t pnx , nt q “ ∇Γpx, tq. Or ~n, la normale au
point px, tq, étant un vecteur unitaire, la relation (1.16) assure que Bx Γpx, tq ‰ 0. Une
application directe du théorème des fonctions implicites assure que Γpx, tq “ 0 ô x “ Xptq
nt
avec X 1 ptq “ ´ . La relation (1.15) s’écrit alors:
nx
nt
rf puqs “ σptqrus , σptq “ X 1 ptq “ ´ . (1.17)
nx
ˆ ˙
u2
Exemple 1.3.1. Equation de Burgers. On considère à nouveau l’équation Bt u`Bx “ 0.
2
On cherche une solution sous la forme:
"
uL si x ă σt ,
upx, tq “
uR si x ą σt ,

avec uL , uR P R et σ P R, c’est à dire une discontinuité entre deux états constants se


propageant à une vitesse constante σ. La relation (1.17) donne:
ˆ 2˙ ˆ 2˙
uR uL uL ` uR
´ “ σ puR ´ uL q ñ σ “ .
2 2 2
uL ` uR
La solution est donc un choc se propageant à vitesse σ “ .
2

14
Le principal problème avec l’introduction de solutions faibles est que l’on
ˆ perd
˙ l’unicité
u2
des solutions. Par exemple, si l’on revient à l’équation de Burgers Bt u ` Bx “ 0, avec
2
la condition initiale upx, 0q “ 0, on peut aisément vérifier que pour tout s ą 0, la fonction
$

’ 0 si x ă ´2s ,
&
´2s si ´ 2s ă x ă 0 ,
upx, tq “

’ 2s si 0 ă x ă 2s ,
%
0 si x ą 2s ,

est solution faible. On démontre ainsi qu’il y a une infinité (non dénombrable) de solutions
faibles. Nous allons voir que l’on peut récupérer l’unicité des solutions à travers le concept
de solution entropique.

1.4 Solutions entropiques


On cherche un critère qui va permettre de sélectionner l’unique solution physique parmi
toutes les solutions faibles. Les lois de conservation du type Bt u ` Bx f puq “ 0 sont souvent
des approximations de lois du type

Bt u ` Bx f puq “ εBxx u , (1.18)

où l’on a négligé le terme de diffusion εBxx u (viscosité physique). On peut considérer que les
"bonnes" solutions faibles sont celles qui sont limites de solutions uε de (1.18) lorsque ε tend
vers 0. Il semble donc pertinent d’identifier une information propre aux solutions de (1.18)
qui persiste lors du passage à la limite uε Ñ u.

Soit E : R ÝÑ R une fonction convexe de classe C 2 (on a donc E 2 ą 0). En multipli-


ant (1.18) par E 1 puε q:

E 1 puε qBt uε ` E 1 puε qf 1 puε qBx uε “ εE 1 puε qBxx uε ,

puis en définissant F puq telle que F 1 puq “ E 1 puqf 1puq, on obtient:

Bt Epuε q ` Bx F puεq “ εBxx pEpuε qq ´ εpBx uε q2 E 2 puε q , (1.19)

et donc Bt Epuε q`Bx F puεq´εBxx pEpuε qq ď 0. Cette forme conservative nous permet d’espérer
passer à la limite dans une forme faible. Formellement, en passant à la limite ε Ñ 0:

Bt Epuq ` Bx F puq ď 0 .

Définition 1.4.1: Entropie

Soit E une fonction strictement convexe de R dans R de classe C 1 . On dit que E est
une entropie de l’équation Bt u ` Bx f puq “ 0 s’il existe une fonction F de R dans R de
classe C 1 telle que F 1 “ E 1 f 1 . La fonction F est appelée flux d’entropie.

15
Remarque 1.4.1. Cette définition s’étend aux fonctions E qui sont seulement continues.
Dans ce cas, on peut définir un flux F via l’introduction de suites régularisantes.
Exemple 1.4.1. Entropies de Kruzkov. Les entropies de Kruzkov sont définies par
Epuq “ |u ´ k| pour żtout u P R, où k est une constante réelle. Si on pose F puq “
u
f 1 puqEpuq´f 1pkqEpkq´ f 2 pyqEpyqdy, on peut montrer que F puq “ pf puq´f pkqqsgnpu´kq.
k

Définition 1.4.2: Solution entropique

Soit u P L8
loc pRˆs0, T rq. On dit que c’est une solution entropique du problème de
Cauchy "
Bt u ` Bx f puq “ 0 ,
(1.20)
upx, 0q “ u0 pxq ,
si pour toute entropie continue convexe E associée à un flux F , on a:
żTż ż
pEpuqBt φ ` F puqBxφq dxdt` Epu0 pxqqφpx, 0qdx ě 0 , @φ P DpRˆr0, T rq , φ ě 0 .
0 R R

La notion de solution entropique contient la notion de solution faible:


Proposition 1.4.1

Soit u P L8 8
loc pRˆs0, T rq une solution entropique de (1.20), avec u0 P L pRq. Alors u
est solution faible de (1.20).

Proof. Supposons que u et u0 soient à valeurs dans un intervalle borné ra, bs Ă R. Ecrivons
l’inégalité d’entropie pour Epuq “ |u ´ a|.
żTż ż
ppu ´ aqBt φ ` pf puq ´ f paqqBx φq dxdt ` pu0 pxq ´ aqφpx, 0qdx ě 0
0 R R
żTż ż
ñ puBt φ ` f puqBx φq dxdt ` u0 pxqφpx, 0qdx
0 R R
¨ ˛
ż ż żTż ż
˚ T ‹
˚
´ ˝a Bt φdxdt ` f paq Bt φdxdt ` a φpx, 0qdx‹ ‚ě 0 .
looooooooooooooooooooooooooooooooomooooooooooooooooooooooooooooooooon
0 R 0 R R
“0

Un calcul analogue impliquant Epuq “ |u ´ b| donne l’inégalité dans l’autre sens. On obtient
donc:
żTż ż
puBt φ ` f puqBxφq dxdt ` u0 pxqφpx, 0qdx “ 0 , @φ P DpR ˆ r0, T rq , φ ě 0 ,
0 R R

Pour étendre ce résultat à φ P DpR ˆ r0, T rq quelconque, on écrit φ “ φ` ´ φ´ avec φ` “


|φ| ` φ |φ| ´ φ
ě 0 et φ´ “ ě 0, et on conclut avec la linéarité de l’intégrale.
2 2
16
Remarque 1.4.2.

‚ Les inégalités d’entropie pour les seules fonctions de la forme Epuq “ |u ´ k| suffisent
à obtenir toutes les inégalités d’entropie. Nous ne démontrerons pas ce résultat ici.

‚ La définition des solutions entropiques introduit de l’irréversibilité: si upx, tq est une


solution faible entropique sur R ˆ r0, T s, alors vp´x, T ´ tq est solution faible de
"
Bt v ` Bx f pvq “ 0 ,
(1.21)
vpx, 0q “ upx, ´T q ,

mais ce n’est pas une solution entropique, car on a Bt Epuq ` Bx F puq ě 0 au sens faible
(sauf si bien sûr ces inégalités sont des égalités).

Théorème 1.4.1: Solutions entropiques discontinues

Soit u P L8loc pRˆs0, T rq de classe C en dehors d’une courbe régulière Γ. Alors u est
1

solution faible entropique si et seulement si:

‚ C’est une solution classique

‚ Le long de Γ on a l’inégalité:

σrEpuqs ě rF puqs (1.22)


nt
pour toute entropie E et flux associé F , avec σ “ ´ et ~n qui pointe vers la
nx
droite pnx ą 0q.

Exercice 1.4.1. Cas où f est strictement convexe. Considérons le cas d’un flux f strictement
convexe. Pour un couple entropie/flux pE, F q, on considère la fonction définie sur Rz tu´ u.

f puq ´ f pu´ q ` ´
˘ ` ´
˘
Zpuq “ Epuq ´ Epu q ´ F puq ´ F pu q . (1.23)
u ´ u´
1. Montrer que Z peut être prolongée par continuité en u´ .
On a:
f puq ´ f pu´ q
lim´ “ f 1 pu´ q ,
uÑu u ´ u´
Epuq ´ Epu´ q
lim´ “ E 1 pu´ q
uÑu u ´ u´
F puq ´ F pu´ q
lim´ ´
“ F 1 pu´ q “ f 1 pu´ qE 1 pu´ q ,
uÑu u´u
ce qui permet de conclure en posant Zpu´ q “ 0.

17
2. Montrer que Z est décroissante. Le calcul donne:
f 1 puq pu ´ u´ q ´ pf puq ´ f pu´ qq ` ˘
Z 1 puq “ Epuq ´ Epu ´
q
pu ´ u´ q2
f puq ´ f pu´ q 1
` E puq ´ F 1 puq .
u ´ u´
f 1 puq pu ´ u´ q ´ f puq ` f pu´ q ` ´
˘
“ Epuq ´ Epu q
pu ´ u´ q2
f puq ´ f pu´ q ´ f 1 puqpu ´ u´ q 1
` E puq .
u ´ u´
pf 1 puq pu ´ u´ q ´ f puq ` f pu´ qq pEpuq ´ Epu´ q ´ pu ´ u´ q E 1 puqq
“ .
pu ´ u´ q2
Par suite, f et E étant strictement convexes:
` ˘ ` ˘
f pu´ q ď f puq ` u´ ´ u f 1 puq et Epu´ q ď Epuq ` u´ ´ u E 1 puq ,
de sorte que les deux termes du produit au numérateur sont de signe opposé, d’où le
résultat.
3. Montrer qu’un choc impliquant deux états pu` , u´ q est une solution faible entropique
si et seulement si Zpu` q ě 0. En déduire que ce choc est admissible si et seulement si
u´ ą u` .
Nous avons vu via (1.22) qu’un tel choc correspond à une solution faible entropique si
et seulement si
f pu` q ´ f pu´ q ` ` ´
˘ ` ` ´
˘
“ σ et σ Epu q ´ Epu q ě F pu q ´ F pu q ,
u` ´ u´
où σ est la vitesse du choc. En d’autres termes, un choc est admissible si et seulement
si:
f pu` q ´ f pu´ q ` ` ´
˘ ` ` ´
˘
Epu q ´ Epu q ´ F pu q ´ F pu q ě 0,
u` ´ u´
c’est à dire Zpu` q ě 0. Par suite, Z étant décroissante et Zpu´ q “ 0, on aboutit à la
condition nécessaire et suffisante u` ă u´ . A noter que cette condition est indépen-
dante du choix du couple pE, F q dans la définition de Z.
4. En déduire que l’inégalité d’entropie au sens faible pour un seul couple entropie/flux
suffit à obtenir toutes les inégalités d’entropie.
Supposons que la solution u soit entropique au sens faible pour un couple entropie/flux
pE0 , F0 q. Si u est régulière, alors elle vérifie toutes les inégalités d’entropie au sens fort,
c’est à dire qu’on a Bt Epuq ` Bx pF puqq pour tout couple entropie/flux. Supposons que
la solution consiste en un choc séparant deux états u` , u´ à travers une courbe Γ.
Considérons alors la fonction Z0 définie par (1.23) construite sur la base du couple
pE0 , F0 q. L’inégalité d’entropie étant vérifiée au sens faible pour ce couple particulier,
on a u´ ą u` . Par suite, en vertu de ce qui précède, Zpu` q ě 0 pour tout couple
entropie/flux pE, F q, c’est à dire que (1.22) est vérifiée pour tout couple entropie/flux
pE, F q.
18
Nous avons donc la propriété suivante:
Proposition 1.4.2: Cas d’un flux strictement convexe
Si f est strictement convexe, l’inégalité d’entropie au sens faible pour une seule entropie
strictement convexe E suffit à obtenir toutes les inégalités d’entropie. En particulier,
la condition de choc entropique est équivalente à u´ ą u` le long de toute courbe de
discontinuité Γ.
ˆ ˙
u2
Exemple 1.4.2. Equation de Burgers. Revenons à l’équation de Burgers Bt u ` Bx ,
2
avec des solutions faibles de la forme:
"
uL si x ă σt ,
upx, tq “
uR si x ą σt ,
uL ` uR
avec σ “ . D’après le résultat précédent, on sait déjà que la solution est entropique
2
si et seulement si uL ą uR . On peut vérifier ce résultat directement, en revenant à (1.22),
2
avec par exemple Epuq “ u2 et F puq “ u3 . On a:
3
uL ` uR ` 2 ˘ 2` 3 ˘
´σrEpuqs ` rF puqs ď 0 ô ´ uR ´ u2L ` uR ´ u3L ď 0
2 3
1` 3 3
˘
ô u ´ uL ď 0 .
6 R
Cette inégalité est vérifiée ssi uL ą uR .

1.5 Résultats d’existence et d’unicité


Commençons par introduire l’espace fonctionnel des fonctions à variation bornée: Soit Ω un
ouvert de Rd . Pour u P L1loc pRq, on définit la variation totale de u par:
"ż *
1 d
T VΩ puq “ sup u divpϕq , ϕ P Cc pΩq , }ϕ}L8pΩqd ď 1 .

Définition 1.5.1
Une fonction u P L1loc pΩq est dite à variation bornée sur Ω si T VΩ puq ă `8. On note:
(
BV pΩq “ u P L1loc pΩq , T VΩ puq ă `8 . (1.24)

Nous admettrons les trois résultats fondamentaux suivants:


Théorème 1.5.1: Unicité
Pour toute donnée initiale u0 P L8 pRd q, il existe une unique solution entropique de
(1.13) u P L8 pRd ˆ R` q X C 0 pr0, T s, L1loc pRd qq.

19
Théorème 1.5.2: Existence
Soit u0 P pL8 X L1 X BV q pRd q. Alors il existe une solution faible entropique u du
problème (1.13), avec u P L8 pRd ˆ R` q X C 0 pr0, T s, L1locpRd qq pour tout T ą 0. Cette
solution vérifie les estimations suivantes:

i) }up., tq}L8pRd q ď }u0}L8 pRd q p.p. t ě 0.

ii) up., tq P BV pRd q et T V pup., tqq ď T V pu0 q p.p. t ě 0.

iii) }up., tq ´ up., sq}L1pRd q ď C T V pu0 q|t ´ s|, avec t, s ě 0.


ż ż
iv) upx, tqdx “ u0 pxqdx p.p. t ě 0.
Rd Rd

Théorème 1.5.3: Unicité (Kruzkov)

Soient u et v deux solutions entropiques de donnée initiale u0 et v0 dans


pL8 X L1 X BV q pRd q. On pose M “ max|λ|ď}u}L8 |f 1 pλq|. Alors, pour tout R ą 0
|λ|ď}v}L8
et tout t ě 0:
ż ż
|upx, tq ´ vpx, tq|dx ď |u0 pxq ´ v0 pxq|dx . (1.25)
xďR |x|ďR`M t

Remarque 1.5.1.
‚ Si u0 “ v0 , alors u “ v p.p. : c’est une conséquence importante du théorème, qui assure
l’unicité des solutions (p.p.) pour une condition initiale donnée.
‚ Les domaines d’intégration dans (1.25) font ressortir le cône d’influence : la solution au point
px, tq ne peut être impactée que par les valeurs comprises dans l’intervalle sx ´ Mt, x ` Mtr
au temps 0 (M majore la vitesse de propagation de l’information).
‚ Ce résultat se démontre en utilisant les entropies de Kruzkov.
‚ L’existence est assurée par u0 P L8 pRd q, et l’unicité par u0 P pL8 X L1 X BV q pRd q.
‚ L’opérateur solution est contractant dans L1 pRd q. Si u0 , v0 P L8 pRd q et u0 ´ v0 P L1 pRd q,
alors up., tq ´ vp., tq P L1 pRd q et:
}up., tq ´ vp., tq}L1 pRd q ď }u0 ´ v0 }L1 pRd q .
Rappelons aussi que dans ces conditions on a la conservation de la masse (il suffit d’intégrer
l’équation (1.13) en espace):
ż ż
pupx, tq ´ vpx, tqq dx “ pu0 pxq ´ v0 pxqq dx,
Rd Rd

ce qui implique, en particulier que si u0 P L1 pRd q, alors up., tq P L1 pRd q et


ż ż
}up., tq}L1pRd q ď }u0}L1 pRd q , upx, tqdx “ u0 pxqdx .
Rd Rd
20
‚ On déduit de ce qui précède le principe du maximum:
u0 ď v0 p.p. ñ uptq ď vptq p.p. . (1.26)
En effet:
ż ż
}up., tq ´ vp., tq}L1 ď }u0 ´ v0 }L1 “ pv0 ´ u0 qdx “ pvp., tq ´ up., tqqdx ,
Rd Rd
ż
et d’autre part pvp., tq ´ up., tqqdx ď }vp., tq ´ up., tq}L1 . On en déduit que
ż ” Rd ı
|vp., tq ´ up., tq| ´ pvp., tq ´ up., tqq “ 0, c’est à dire vp., tq ´ up., tq “ |vp., tq ´ up., tq| p.p.,
Rd
ou encore vp., tq ě up., tq p.p.

1.6 Le problème de Riemann en 1d


On s’intéresse à l’unique solution entropique du problème de Cauchy suivant (appelé prob-
lème de Riemann): $
& Bt u ` Bx f puq
" “ 0,
uL si x ă 0 , (1.27)
% upx, 0q “
uR si x ą 0 .

Définition 1.6.1
‚ On appelle onde de choc une solution de (1.27) composée de deux états con-
stants uL , uR séparés par une discontinuité de vitesse σ P R:
"
uL si x ă σt ,
upx, tq “
uR si x ą σt .

‚ On appelle onde de détente une solution continue de (1.27) et C 1 par morceaux


du type upx, tq “ vpx{tq.

Lemme 1.6.1
Les ondes de détente sont du type upx, tq “ cte ou bien upx, tq “ vpx{tq, où vpξq “
pf 1 q´1 pξq.

Proof. Soit upx, tq “ vpx{tq une solution. On injecte dans l’équation (1.27):
x 1
´ 2 v 1 px{tq ` f 1 pvpx{tqq v 1 px{tq “ 0 , @px, tq P pR ˆ R˚` q .
t t
par suite, soit v est constante, auquel cas upx, tq “ cte, ou bien pf 1 ˝ vqpξq “ ξ pour tout
ξ.
Remarque 1.6.1. Une solution qui ne dépend que de x{t est dite auto-semblable ou
auto-similaire.
21
1.6.1 La condition d’Oleinik
Il s’agit d’un critère graphique pour savoir si une discontinuité entre deux états distincts u´
et u` se propageant à la vitesse σ est bien entropique. Rappelons que par la condition de
Rankine-Hugoniot (1.17), on a:

´ σrus ` rf puqs “ 0 . (1.28)

D’autre part, la solution est entropique si et seulement si

@k P R , ´σr|u ´ k|s ` rsgnpu ´ kqpf puq ´ f pkqqs ď 0 . (1.29)

Considérons alors l’intervalle I d’extrémités u´ et u` .

‚ Si k R I, alors on a (1.29) grâce à (1.28).

‚ Si k P I, alors k “ au´ ` p1 ´ aqu` avec a P r0, 1s. On a:

r|u ´ k|s “ |u` ´ k| ´ |u´ ´ k| “ |apu´ ´ u` q| ´ |p1 ´ aq pu` ´ u´ q |


“ p2a ´ 1q|u´ ´ u` | .

D’autre part,

rsgnpu ´ kqpf puq ´ f pkqqs “ sgnpu` ´ kqpf pu` q ´ f pkqq ´ sgnpu´ ´ kqpf pu´ q ´ f pkqq
“ sgnpu` ´ u´ qpf pu` q ´ f pkqq ´ sgnpu´ ´ u` qpf pu´ q ´ f pkqq
“ sgnpu` ´ u´ q pf pu` q ` f pu´ q ´ 2f pkqq .

Ainsi, la solution est entropique si et seulement si:

rf puqs
´ p2a ´ 1q |u` ´ u´ | ` sgnpu` ´ u´ q pf pu` q ` f pu´ q ´ 2f pkqq ď 0 ,
rus
ou encore:

sgnpu` ´ u´ q p2a ´ 1q rf puqs ě sgnpu` ´ u´ q pf pu` q ` f pu´ q ´ 2f pkqq .

On vérifie que cette condition s’écrit:


Si u´ ă u` : f pau´ ` p1 ´ aqu` q ě af pu´ q ` p1 ´ aqf pu` q.
Si u` ă u´ : f pau´ ` p1 ´ aqu` q ď af pu´ q ` p1 ´ aqf pu` q.

Ces deux configurations sont représentés dans la Figure 8. On en déduit la proposition


suivante:
Proposition 1.6.1: Critère d’Oleinik

Une discontinuité pu´ , u` q est entropique ssi:


Si u´ ă u` : le graphe de f|ru´ ,u` s est au dessus de sa corde.
Si u` ă u´ : le graphe de f|ru` ,u´ s est au dessous de sa corde.

22
f puq non admissible f puq
admissible

admissible non admissible

u´ u´ u` u u` u´ u´ u

Figure 8: Critère d’Oleinik. Les cas u´ ă u` (gauche) et u` ă u´ (droite)

On déduit de ce critère le critère de Lax (exercice):


Proposition 1.6.2: Critère de Lax.

Si Γ est une discontinuité entre u´ px, tq et u` px, tq se propageant à la vitesse σptq,


alors:
f 1 pu` px, tqq ď σptq ď f 1 pu´ px, tqq . (1.30)
Ces conditions sont communément appelées inégalités de Lax.

1.6.2 Résolution pratique des problèmes de Riemann


Théorème 1.6.1
L’unique solution faible entropique des problèmes de Riemann (1.27) est constituée des
deux états constants uL et uR , séparés par des chocs et des détentes. Cette solution
est auto-similaire.

Proof.

i) Cas d’un flux strictement convexe.

‚ Si uR ă uL , f étant convexe, le graphe de f|ruR ,uL s est en dessous de sa corde. Il


f puR q ´ f puL q
s’agit donc d’un choc admissible qui se propage à la vitesse σ “ :
uR ´ uL
c’est l’unique solution.
‚ Si uL ă uR , on établit que la solution est une onde de détente. Notons que comme
f 2 ą 0, f 1 est monotone et est donc inversible. On pose alors:
$
& uL si x´ă
’ ¯σL t ,
1 ´1 x
upx, tq “ pf q si σL t ă x ă σR t ,

% t
uR si x ą σR t ,

23
avec la condition de continuité uL “ pf 1 q´1 pσL q et uR “ pf 1 q´1 pσR q, c’est à dire
σL “ f 1 puL q et σR “ f 1 puR q. Notons que l’on a bien σL ă σR car f 1 est croissante.

ii) Cas d’un flux strictement concave. Si f est strictement concave (donc f 2 ă 0, f 1
décroissante), un raisonnement analogue donne un choc admissible dans le cas uL ă uR
et une onde de détente si uL ą uR .

iii) Cas général. On supposera ici que f admet un nombre fini de points d’inflexion.

‚ Si uL ă uR , on construit l’enveloppe convexe inférieure fc de f sur le segment


I “ ruL , uR s. Par définition:

fc “ sup g.
g convexe sur I
gďf sur I

Là où f ” fc , on a une détente, et là où f ą fc , on a un choc admissible. Une


illustration est donnée en Figure (9).
u
f puq uR
b
u3
b

u2
b

b
b

u1
uL
uL u1 u2 u3 uR u x
choc choc
détente détente

Figure 9: Résolution pratique du problème de Riemann dans le cas général (uL ă uR ). Con-
struction de l’enveloppe convexe inférieure fc (gauche), et connexion des états intermédiaires
(droite).

‚ Si uR ă uL , on construit l’enveloppe concave supérieure f c de f sur le segment


I “ ruR, uL s:
fc “ sup g.
g concave sur I
gěf sur I

Là où f ” f c , on a une détente, et là où f ą f c , on a un choc admissible. Une


illustration est donnée en Figure (10).

Exemple 1.6.1.

24
u
f puq uL
b
u3
b
b
u2
b

u1
uR
uR u1 u2 u3 uL u x
choc choc
détente détente

Figure 10: Résolution pratique du problème de Riemann dans le cas général (uR ă uL ). Con-
struction de l’enveloppe concave supérieure f c (gauche), et connexion des états intermédiaires
(droite).

Cas d’un flux strictement convexe : équation de Burgers


On considère l’équation de Burgers associée à la donnée initiale suivante:
$ ˆ 2˙

’ u
& Bt u ` Bx “ 0,
"2 (1.31)

’ 2 si x ă 1 ,
% upx, 0q “
´1 si x ą 1 .

Un calcul préliminaire donne que les caractéristiques sont données par


"
y ` 2t si y ă 1 ,
xptq “ y ` u0 pyqt “
y ´ t si y ą 1 .

Le flux f est strictement convexe, et on a uL “ 2 ą uR “ ´1. On sait donc que l’on peut
2 ` p´1q
connecter ces deux états par un choc admissible de vitesse “ 1{2, centré en x “ 1.
2
1
L’équation de l’onde de choc est donc donnée par γptq “ 1 ` t. La solution est donnée par:
2
$
& 2 si x ă 1 ` 1 t ,

upx, tq “ 2
’ 1
% ´1 si x ą 1 ` t .
2
La Figure (11) illustre les caractéristiques de part et d’autre de la courbe de choc.
On considère à présent la condition initiale:
"
´1 si x ă 1 ,
upx, 0q “
2 si x ą 1 .

25
γptq “ 1 ` t{2
t

1 x

Figure 11: Caractéristiques et courbes de choc.

On est dans la configuration uL ă uR avec f convexe: la solution est donc une onde de
détente reliant uL à uR . On a f puq “ u, de sorte que pf 1 q´1 puq “ u, σL “ f 1 puL q “ ´1 et
σR “ f 1 p2q “ 2. La solution est donc:
$
& ´1 si x ă 1 ´ t ,

x´1
upx, tq “ si 1 ´ t ă x ă 1 ` 2t ,

% t
2 si x ą 1 ` 2t .
Cas d’un flux strictement concave.
On considère l’équation suivante:
$
& Bt u ` Bx pf"puqq “ 0 ,
1 si x ă 1 , (1.32)
% upx, 0q “
0 si x ą 1 ,
up2 ´ uq
avec f puq “ . Le flux f est concave pf 1 puq “ 1 ´ u, f 2 “ ´1q. Les caractéristiques
2
sont données par: "
1 y si y ă 1 ,
xptq “ y ` f pu0 pyqqt “ (1.33)
y ` t si y ą 1 .
On a uL “ 1 ą uR “ 0 avec f concave: la solution est une onde de détente centrée en x “ 1.
On a pf 1 q´1 puq “ 1 ´ u, σL “ f 1 p1q “ 0 et σR “ f 1 p0q “ 1, de sorte que:
$
& 1 si x ă 1 ,

x´1
upx, tq “ 1´ si 1 ă x ă 1 ` t ,

% t
0 si x ą 1 ` t .
Considérons à présent la condition initiale
"
0 si x ă 1 ,
upx, 0q “ (1.34)
1 si x ą 1 .
On vérifie que la solution est un choc de vitesse σ “ 1{2 reliant uL “ 0 à uR “ 1.
26
Cas d’un flux ni convexe ni concave. On considère l’équation Bt u ` Bx f puq “ 0, avec
f puq “ sinpπuq. On considère dans un premier temps la condition initiale suivante:
"
1 si x ă 0 ,
upx, 0q “
´1 si x ą 0 ,

fc

´u˚ u˚
b
fc
Figure 12: Enveloppe convexe inférieure fc et enveloppe concave supérieure f c associée à f
sur r´1, 1s.

On a uL “ 1 ą uR “ ´1, donc on considère l’enveloppe concave supérieure f c associée


à f sur r´1, 1s. En notant u˚ le point où f c coïncide avec f (caractérisé par f 1 pu˚ q “
f p´1q ´ f pu˚ q
, c’est à dire π cospπu˚ q p1 ` u˚ q “ sinpπu˚ q), on a f ą f c sur l’intervalle
´1 ´ u˚
r´1, u˚s et f ” f c sur ru˚, 1s. Dans un premier temps, un choc admissible de vitesse
f p´1q ´ f pu˚q
σ˚ “ va connecter uR “ ´1 à u˚ . Par suite, on connecte l’état u˚ à uL “ 1
´1 ´ u˚
par une détente délimitée par les ondes σ ˚ “ f 1 pu˚ q et σL “ f 1 puL q “ π cospπq “ ´π. A
noter que l’on a pf 1 q´1 pξq “ arccospξ{πq{π pour tout ξ Ps0, πr. Au final, la solution est
donnée par: $

& 1 si x ă ´πt ,
1 ´x¯
upx, tq “ arccos si ´ πt ă x ă σ ˚ t ,

% π tπ
´1 si x ą σ ˚ t ,

˚ ˚ f p´1q ´ f pu˚ q 1
avec σ “ π cospπu q “ ˚
. On vérifie que arccos p´1q “ 1 “ uL et
ˆ ˚˙ ´1 ´ u π
1 σ ˚
arccos “ u , de sorte que la relation de continuité est satisfaite pour l’onde de
π π
détente. Enfin, on a que uR “ ´1 et u˚ sont bien connectés par un choc admissible de
vitesse σ ˚ .
Considérons à présent la condition initiale:
"
´1 si x ă 0 ,
upx, 0q “
1 si x ą 0 ,

On a uL “ ´1 ă uR “ 1, et on considère donc à présent l’enveloppe convexe inférieure


fc sur r´1, 1s (la position du point ´u˚ se déduit de l’imparité de f ). La solution sera
constituée dans un premier temps d’une détente entre uL “ ´1 et ´u˚ , impliquant les
vitesses σL “ f 1 p´1q “ π cosp´πq “ ´π et σ ˚ “ f 1 p´u˚ q “ π cospπu˚q. On a ensuite un

27
f p´u˚q ´ f p1q sinpπu˚ q
choc reliant ´u˚ à uR “ 1 de vitesse “ “ σ ˚ . En notant que
´u˚ ´ 1 u˚ ` 1
pf 1 q´1 pξq “ ´ arccospξ{πq{π pour tout ξ Ps ´ π, 0r, la solution est donc:
$

& ´1 si x ă ´πt ,
1 ´x¯
upx, tq “ ´ arccos si ´ πt ă x ă σ ˚ t ,

% π tπ
1 si x ą σ ˚ t ,
avec σ ˚ “ π cospπu˚q.

28
Chapter 2

Schémas numériques pour les lois de


conservation scalaires

2.1 Rappels sur les différences finies


On considère le problème suivant:
"
Bt u ` Bx f puq “ 0 ,
(2.1)
upx, t “ 0q “ u0 pxq .

On introduit un maillage de R ˆ R` en notant ∆t ą 0 le pas de temps et ∆x ą 0 le pas


d’espace. Pour tout n P N, on note tn “ n∆t le temps discret et pour tout j P Z, on note
xj “ j∆x les noeuds du maillage. On note unj l’approximation de upxj , tn q : unj « upxj , tn q.
Dans un premier temps, une solution simple consiste à approcher les dérivées par différences
finies, par exemple:

ujn`1 ´ unj
Bt upxj , tn q « ,
∆t
n
f punj`1q ´ f punj´1q
Bx f pupxj , t qq « ,
2∆x
ce qui donne le schéma explicite centré:

ujn`1 ´ unj f punj`1q ´ f punj´1q


` “ 0, (2.2)
∆t 2∆x
mais nous verrons plus tard que cette option n’est pas convenable. Le schéma implicite
centré fournit un autre exemple

ujn`1 ´ unj n`1


f puj`1 n`1
q ´ f puj´1 q
` “ 0,
∆t 2∆x
qui possède de meilleures propriétés, mais plus difficile à mettre en oeuvre.

29
Définition 2.1.1: Schéma numérique
` ˘
Un schéma numérique est une formule F pum j q 0ďmďn`1 “ 0 qui permet de calculer les
jPZ ` ˘
valeurs pujn`1 qjPZ en fonction des valeurs connues aux temps pum j q 0ďmăn`1 . On dit
jPZ
que le schéma est consistant au sens des différences finies si l’erreur de troncature
E “ F puptm , xj qq tend vers 0 lorsque ∆t et ∆x tendent vers 0 (éventuellement avec
des conditions du type ∆t{∆x fixé) si et seulement si u est solution régulière de (2.1).

On se restreindra ici à l’étude des schémas à deux niveaux en temps, c’est à dire de la
forme:
` m
˘ ujn`1 ´ unj
F puj q 0ďmďn`1 “ “ Gpunj qjPZ (schémas explicites) ,
jPZ ∆t
ou
` m
˘ ujn`1 ´ unj
F puj q 0ďmďn`1 “ “ Gpujn`1qjPZ (schémas implicites) .
jPZ ∆t
Si G ne dépend que des punj1 q pour |j ´ j 1 | ď K, on dit que le schéma est à 2K ` 1 points.
Par convention, on écrira les schémas sous la forme:
´ ¯
n n`1 n npou n`1q
F puj q “ uj ´ uj ` ∆tG uj 1 . (2.3)
j 1 PZ

1
Attention, avec cette écriture, l’erreur de troncature devient E “ F pupxj , tn qq.
∆t
Exemple 2.1.1. Autres exemples de schémas numériques:
‚ Schéma de Lax-Friedrichs.
ˆ n ˙ ˆ ˙
n`1
uj´1 ` unj`1 f punj`1q ´ f punj´1q
uj ´ ` ∆t “ 0.
2 2∆x
A titre d’exercice, on pourra mettre ce schéma sous la forme (2.3).
„ 
f punj`1q ´ f punj´1q unj´1 ´ 2unj ` unj`1
ujn`1 ´ unj ` ∆t ´ “ 0.
2∆x 2∆t

‚ Schéma de Lax-Wendroff. Soit u, la solution de (2.1). On a:


1
upx, t ` ∆tq “ upx, tq ` ∆tBt upx, tq ` p∆tq2 Bt2 upx, tq ` Op∆t3 q .
2
1
“ upx, tq ´ ∆tBx f pupx, tqq ` p∆tq2 Bt p´Bx f pupx, tqqq ` Op∆t3 q .
2
On a: Bt p´Bx f puqq “ ´Bx pf puqBt uq “ Bx pf 1 puqBx f puqq. On en déduit le schéma
1

suivant:
ˆ ˙
n`1 n
f punj`1q ´ f punj´1q
uj “ uj ´ ∆t
2∆x
f pun n
j`1 q´f puj q f pun n
j q´f puj´1 q
1 f 1 puj`1{2q ` f 1 puj´1{2q
` p∆tq2 ∆x ∆x
,
2 ∆x
30
avec plusieurs choix possibles pour f 1 puj`1{2q, par exemple:
ˆ ˙
1 1
unj`1 ` unj f 1 punj`1q ` f punj q
f puj`1{2 q “ f , f 1 puj`1{2q “ .
2 2

Définition 2.1.2: Ordre d’un schéma


On dit que le schéma est d’ordre p en temps et q en espace si pour la solution régulière
de l’équation on a:
E “ Op∆tp q ` Op∆xq q .

Exemple 2.1.2. Montrons que le schéma explicite centré (2.2) est d’ordre 1 en temps et 2
en espace. L’erreur s’obtient en remplaçant unj par upxj , tn q dans le schéma:

upxj , tn`1 q ´ upxj , tn q f pupxj`1, tn qq ´ f pupxj´1, tn qq


E“ ` “ 0.
∆t 2∆x
On effectue alors des développements de Talyor:

upxj , tn`1 q “ upxi , tn q ` ∆tBt upxj , tn q ` Op∆t2 q ,


∆x2 2
f ˝ upxj`1 , tn q “ f ˝ upxj , tn q ` ∆xBx pf ˝ uqpxj , tn q ` Bx pf ˝ uqpxj , tn q ` Op∆x3 q ,
2
2
∆x
f ˝ upxj´1 , tn q “ f ˝ upxj , tn q ´ ∆xBx pf ˝ uqpxj , tn q ` B 2 pf ˝ uqpxj , tn q ` Op∆x3 q ,
2 x
de sorte que:
E “ pBt u ` Bt pf puqqq pxj , tn q ` Op∆tq ` Op∆x2 q .
On en déduit que E ÝÑ 0 si et seulement si u est solution et dans ce cas E “ Op∆tq `
∆t,∆xÑ0
Op∆x2 q.

A titre d’exercice, on pourra montrer que le schéma de Lax-Wendroff est d’ordre 2 en


temps et en espace.
Voir examen.

2.2 Stabilité L2 dans le cas linéaire


Considérons le cas particulier de l’équation linéaire
"
Bt u ` aBx u “ 0 , a P R ,
(2.4)
upx, t “ 0q “ u0 pxq .
Dans cette section on se limite à l’analyse de schémas explicites à trois points:
n`1
uj`1 “ c´1 unj´1 ` c0 unj ` c1 unj`1 .

31
L’erreur de troncature est donnée par:

upxj , tn`1 q ´ c´1 upxj´1 , tn q ´ c0 upxj , tn q ´ c1 upxj`1, tn q


E“ . (2.5)
∆t
Après développement de Taylor:
ˆ ˙
1 ´ c´1 ´ c0 ´ c1 ∆x ∆x2
E“ upxj , tn q ` Bt upxj , tn q ` pc´1 ´ c1 q Bx upxj , tn q ` O .
∆t ∆t ∆t

∆t
On supposera dans la suite le rapport ∆t{∆x fixe et on notera λ “ . On déduit de
∆x
l’égalité précédente que :
ˆ ˙ #
c´1 ` c0 ` c1 “ 1
E ÝÑ 0 ssi u est solution ô ∆t
∆t,∆xÑ0 c´1 ´ c1 “ a “ aλ .
∆x
Ceci permet de paramétrer les schémas consistants par q P R en posant:
1 1
c0 “ 1 ´ q , c´1 “ pq ` aλq , c1 “ pq ´ aλq , (2.6)
2 2
de sorte que:
aλ ` n ˘ q` n ˘
ujn`1 “ unj ´ uj`1 ´ unj´1 ` uj`1 ´ 2unj ` unj´1 (2.7)
2 2
On se pose alors la question de savoir si cette classe de schémas fournit des approximations
convenables de la solution, le sens du mot "convenable" étant bien sûr à définir. Une premier
critère essentiel est la notion de stabilité Lp . qui garantit un contrôle de la norme Lp de la
solution discrète. Nous nous intéresserons ici à la stabilité L2 . Comme nous le verrons
plus tard, ce critère peut généralement s’interpréter d’un point de vue physique par une
décroissance de l’énergie discrète associée au schéma.
Définition 2.2.1: Stabilité L2
˜ ¸1{2
ÿ
On définit }unj }L2 “ ∆x|unj |2 la norme L2 de la solution approchée (constante
jPZ
par morceaux) au temps tn . On dit que le schéma est stable en norme L2 s’il existe
une constante c ą 0 indépendante de ∆t et ∆x telle que:

}unj }L2 ď c }u0j }L2 , (2.8)

pour tout ∆t, ∆x (avec ∆t{∆x fixé), tout n P N, et toute donnée initiale pu0j q.

Proposition 2.2.1

Le schéma (2.7) est L2 stable ssi λa2 ď q ď 1.

32
ż
1
Proof. On utilise la transformée de Fourier ϕpξqp “ ? e´ixξ ϕpxqdx. Rappelons que si
2π R
l’on note τh l’opérateur de translation défini par pτh φqpxq “ φpx` hq, alors y p
τh φpξq “ φpξqeihξ
.
Notons u pxq la solution constante par morceaux égale à uj si x Psxj´1{2 , xj`1{2r, où l’on a
n n

posé xj`1{2 “ pj ` 1{2q∆x pour tout j P Z. Le schéma se réécrit sous la forme:

aλ q
un`1 pxq “ un pxq ´ pτ∆x un pxq ´ τ´∆x un pxqq ` pτ∆x un pxq ´ 2un pxq ` τ´∆x un pxqq .
2 2
En passant en Fourier:
aλ ` iξ∆x ˘ n q ` iξ∆x ˘ n
pn`1 pξq “ u
u pn pξq ´ e ´ e´iξ∆x u
p pξq ` e ´ 2 ` e´iξ∆x u
p pξq ,
2 2
c’est à dire:
pn`1 pξq “ hpξqp
u un pξq , (2.9)
avec hpξq “ 1 ´ aiλ sinpξ∆xq ` q pcosp∆xq ´ 1q. Par le théorème (égalité) de Plancherel, on
a:
}unj }L2 “ }un }L2 pRq “ }p
un }L2 pRq ď max |hpξq|n }p
u0 }L2 pRq
ξPR

ď max |hpξq|n }u0j }L2 .


ξPR

On en déduit la condition suffisante de stabilité: maxξPR |hpξq| ď 1. On a:

|hpξq|2 “ λ2 a2 sin2 pξ∆xq ` p1 ` q pcospξ∆xq ´ 1qq2 .


En posant y “ sin2 pξ∆x{2q, on obtient |hpξq|2 “ 4λ2 a2 yp1 ´ yq ` p1 ´ 2qyq2. Après une étude
de fonctions sur r0, 1s, on aboutit à:
` ˘
|hpξq|2 ď 1 , @ξ P R ô λ2 a2 ď q ď 1 . (2.10)

Remarque 2.2.1. En fait, cette condition est nécessaire car s’il existe ξ0 tel que |hpξ0 q| ą 1,
on peut construire une condition initiale u0 pxq dont l’un des modes de Fourier est aussi
proche de ξ0 que l’on veut pour ∆x petit.

Remarque 2.2.2.
‚ La condition nécessaire de stabilité maxξPR |hpξq| ď 1 est communément appelée con-
dition CFL (Courant, Friedrichs, Lewy). Cette condition s’écrit:

∆x
λ|a| ď 1 ô ∆t ď .
|a|

Elle peut s’interpréter de la manière suivante. En un point pxj , tn`1 q, on définit le cône
de dépendance numérique comme le plus petit cône de sommet pxj , tn`1 q contenant les
points pxj , tn qjPZ impliqués dans le calcul de ujn`1 (voir Fig. 1).

33
xj´1 xj xj`1
tn`1 b b b b b

∆t

tn b b b b b

∆x

Figure 1: Cône de dépendance numérique en pxj , tn`1 q.

Dans ce cas, s’agissant d’un schéma à trois points, ce cône doit avoir pour autres som-
∆t
mets pxj´1 , tn q et pxj`1 , tn q. La pente du cône numérique est donnée par ˘ .
∆x
Rappelons que la solution exacte en pxj , tn`1 q n’est affectée que par les valeurs com-
1
prises dans le cône de dépendance exact, de pente ˘ (le cône de dépendance exact
|a|
se déduit directement la méthode des caractéristiques; en effet, on a upx, tq “ f px´atq,
de sorte que l’information se propage à la vitesse a - voir Fig. 3).

xj´1 xj xj`1
tn`1

1 1 ∆t
´
|a| |a|

tn

∆x
Figure 2: Cône de dépendance exact en pxj , tn`1 q.

La condition CFL exprime que le cône de dépendance numérique doit contenir le cône
de dépendance exact. En d’autres termes, pour calculer une approximation ujn`1 de
upxj , tn`1 q, les approximations disponibles au temps tn doivent recouvrir la zone qui
affecte la solution exacte en pxj , tn`1 q

‚ Le schéma centré explicite (2.2) correspond au cas q “ 0: ce schéma est donc instable.

‚ Le cas q “ 1 correspond au schéma le plus diffusif (ou visqueux) possible: c’est le


schéma de Lax-Friedrichs, qui est donc stable ssi λ|a| ď 1.

34
‚ Le schéma de Lax-Wendroff correspond au cas q “ λ2 a2 . Il s’agit du schéma le moins
diffusif que l’on puisse construire dans ce contexte. Ce schéma est également stable
sous la condition λ|a| ď 1.
‚ Le schéma décentré amont (ou upwind) en anglais s’écrit:
` ˘
ujn`1 “ unj ´ λa unj ´ unj´1 si a ě 0 ,
` ˘
ujn`1 “ unj ´ λa unj`1 ´ unj si a ď 0 ,
correspond au cas q “ λ|a|. Il est stable sous la condition λ|a| ď 1.
Remarquons que l’on peut écrire:
ujn`1 ´ unj unj`1 ´ unj´1 ∆x unj`1 ´ 2unj ` unj´1
`a “q . (2.11)
∆t ∆x 2λ ∆x2
On peut dégager de cette analyse qu’on a toujours besoin d’introduire un peu de viscosité
numérique afin d’assurer la stabilité des schémas.
L’analyse de stabilité L2 que nous venons de voir, appelée analyse de stabilité au sens
de Von Neumann ne se limite pas au cas du transport linéaire et aux schémas explicites à
trois points. Nous verrons plus tard dans quelle mesure cette analyse peut s’étendre au cas
des systèmes d’équations non linéaires.

2.3 Schémas conservatifs


Pour traiter le cas non linéaire, une première approche naïve consiste à généraliser les schémas
qui précèdent. Considérons par exemple l’équation de Burgers:
pupx, tqq2
Bt upx, tq ` Bx “ 0, (2.12)
2
"
1 si x ă 0 ,
avec donnée initiale u0 pxq “ Ecrivons formellement Bt upx, tq`upx, tqBx upx, tq “
0 si x ě 0 .
0, et essyons le schéma upwind pour cette équation. On a:
∆t n ` n ˘
ujn`1 “ unj ´ uj uj ´ unj´1 .
∆x
On vérifie alors que:
" "
0 1 si j ă 0 n 1 si j ă 0
uj “ ùñ uj “ , @n P N .
0 si j ě 0 0 si j ě 0
La solution numérique consiste donc en un choc stationnaire (i.e. de vitesse σ “ 0). Or, la
1
solution exacte est un choc de vitesse σ “ . Tentons alors une autre approche, en fixant
2 "
∆t 1 si x ď 0 ,
λ“ “ 1 et avec la donnée initiale u0 pxq “
∆x 0 si x ą 0 ,
∆t n ` n ˘
ujn`1 “ unj ´ uj´1 uj ´ unj´1 .
∆x
35
on vérifie (exercice) que `la solution˘ est un choc se propageant à la vitesse σ “ 1.
On a u10 “ u00 ´ u0´1 u00 ´ u0´1 “ 0 ´ 1p0 ´ 1q “ 1 et u11 “ u01 ´ u00 pu01 ´ u00 q “ 1. On
montre par récurrence que:
" "
0 1 si j ă 0 n 1 si j ă n
uj “ ùñ uj “ , @n P N .
0 si j ě 0 0 si j ě n

Conclusion : Pour calculer les solutions faibles (discontinues), il faut abandonner cer-
taines idées des différences finies, valables uniquement pour les solutions régulières.
Définition 2.3.1: Schéma conservatif
Un schéma explicite à trois points est dit conservatif s’il existe une fonction continue
g : R ˆ R Ñ R telle que le schéma s’écrive:
∆t ` n ˘
ujn`1 “ unj ´ n
gj`1{2 ´ gj´1{2 , (2.13)
∆x
avec gj`1{2 “ gpunj , unj`1q. On dit qu’un tel schéma est consistant au sens faible si
gpu, uq “ f puq pour tout u P R. On dit alors que le flux numérique g est consistant.

Remarque 2.3.1.
‚ Le flux numérique g est défini à une constante additive près.
ÿ ÿ ÿ
‚ ujn`1 “ unj “ u0j : cette quantité est conservée. Ce résultat est le pendant
jPZ jPZ jPZ
discret du résultat: ż ż
upx, tqdx “ u0 pxqdx , (2.14)
R R
"
Bt u ` Bx pf puqq “ 0
pour toute solution u du problème avec u0 P L1 pRq.
upx, 0q “ u0 pxq
‚ La formule (2.13) est une formule de Volumes Finis. Le principe de construction
d’un schéma de type Volumes Finis est le suivant: on intègre l’équation sur la maille
sxj´1{2 , xj`1{2rˆstn , tn`1 r
ż tn`1 ż xj`1{2
pBt u ` Bx f puqq dxdt “ 0
tn xj´1{2
ż xj`1{2 ż xj`1{2
n`1
ô upx, t qdx ´ upx, tn qdx
xj´1{2 xj´1{2
ż tn`1 ż tn`1
` f pupxj`1{2, tqqdt ´ f pupxj´1{2, tqqdt “ 0 .
tn tn
ż xj`1{2
1
On pose alors unj “ upx, tn qdx et on définit gj`1{2
n
comme une approximation
∆t xj´1{2
ż tn`1
1
de f pupxj`1{2, tqqdt calculée à partir de unj et unj`1. On obtient alors:
∆t tn
` ˘ ` n ˘
∆x ujn`1 ´ unj ` ∆t gj`1{2 n
´ qj´1{2 “ 0.
36
Notons qu’aucune hypothèse sur la régularité de u n’est nécessaire pour dériver un tel
schéma (si ce n’est l’intégrabilité locale de u).

Définition 2.3.2: Définition plus générale

Un schéma explicite à p2K ` 1q points est conservatif s’il s’écrit

∆t ` n ˘
ujn`1 “ unj ´ gpuj´K`1, ¨ ¨ ¨ , unj`K q ´ gpunj´K , ¨ ¨ ¨ , unj`K´1q ,
∆x
Le schéma est dit consistant si gpu, ¨ ¨ ¨ , uq “ f puq pour tout u P R.

Proposition 2.3.1
Tout schéma explicite consistant est au moins d’ordre 1 au sens des différences finies.

Proof.

upxj , tn`1 q ´ upxj , tn q g pupxj , tn q, upxj`1q, tn q ´ g pupxj´1 , tn q, upxj , tn qq


E“ ` .
∆t ∆x
Notons que

g pupxj , tn q, upxj`1, tn qq “ g pupxj , tn q, upxj ` ∆x, tn qq


` ˘
“ g upxj , tn q, upxj , tn q ` ∆xBx upxj , tn q ` Op∆x2 q
` ˘
“ gpu, uq ` ∆xBx u ` Op∆x2 q By gpu, uq
“ gpu, uq ` ∆xBx uBy g ` Op∆x2 q ,

et un calcul similaire donne:

g pupxj´1 , tn q, upxj , tn qq “ gpu, uq ´ ∆xBx uBy g ` Op∆x2 q

On en déduit
E “ Bt u ` pBx g ` By gq Bx u ` Op∆xq .
Or, gpu, uq “ f puq ùñ f 1 puq “ Bx g ` By g, et par conséquent E “ Op∆tq ` Op∆xq si u est
solution.
Le théorème suivant montre que les schémas conservatifs sont adaptés au calcul des
solutions faibles.

37
Théorème 2.3.1: Théorème de Lax-Wendroff
Si un schéma conservatif consistant converge, alors il converge vers une solution faible
de l’équation (en particulier, les conditions de Rankine-Hugoniot sont vérifiées). Plus
précisément:
Soit pεk q une suite de R` qui tend vers 0. On considère une suite de solutions discrètes
uk px, tq calculées avec le maillage ∆x “ εk et ∆t “ λ∆x avec λ ą 0 fixé. On a:

uk px, tq “ unj pour xj´1{2 ă x ă xj`1{2 et tn ă t ă tn`1 , .


ż
1 xj`1{2
La donnée initiale discrète est uj “
0
u0 pxqdx. On suppose que:
∆x xj´1{2

i) }uk }L8 pR,R` q ď c pour tout k P N.

ii) Du t.q. uk ÝÑ u dans L1loc pR, R` q et p.p. dans R ˆ R` .

Alors u est solution faible.

Proof. Nous aurons besoin des deux Lemmes suivants:


Lemme 2.3.1
Soit u˘ ˘
k px, tq “ uk px ˘ ∆x{2, tq. Alors uk ÝÑ u dans Lloc .
1

Proof. Soit K un compact de R ˆ R` . On a:


ż ż ż ż
˘
|uk ´ u|dxdt “ |uk px ˘ ∆x{2, tq ´ upx, tq|dxdt
K ż ż K

“ |uk px1 tq ´ upx1 ¯ ∆x{2, tq|dx1 dt


1
ż żK :“K¯∆x
` ˘
ď |uk ´ u| ` |upx1 , tq ´ upx1 ¯ ∆x{2, tq| dx1 dt .
K1

Le premier terme tend vers 0 car uk ÝÑ u dans L1loc pR, R` q. Concernant le 2e terme, on sait
que u P L1loc pR, R` q, et par un argument de densité de Cc1 pR, R` q dans L1loc pR, R` q, on obtient:
ż ż
|upx1 , tq ´ upx1 ¯ ∆x{2, tq|dx1 dt ÝÑ 0 ,
K1

d’où le résultat.

38
Lemme 2.3.2
Soit ϕ P Cc1 pR, R` q. On pose, pour k fixé, ϕnj “ ϕpxj , tn q. Soit ϕk px, tq la fonction
constante par morceaux définie par ϕk px, tq “ ϕnj pour xj´1{2 ă x ă xj`1{2 et tn ă t ă
tn`1 . Alors:

i) ϕk ÝÑ ϕ uniformément sur tout compact de R ˆ R` .


ϕk px, t ` ∆tq ´ ϕk px, tq
ii) ÝÑ Bt ϕpx, tq uniformément sur tout compact de Rˆ R` .
∆t
ϕk px ` ∆x{2, tq ´ ϕk px ´ ∆x{2, tq
iii) ÝÑ Bx ϕpx, tq uniformément sur tout compact
`
∆t
de R ˆ R .

Proof. La démonstration est laissée à titre d’exercice.


Passons à la preuve du théorème. Considérons une fonction ϕ P Cc1 pR, R` q quelconque.
On définit ϕnj “ ϕpxj , tn q comme dans le Lemme précédent. On multiplie le schéma (2.3)
par ϕnj puis on somme sur j et sur n:
ÿ” ı
∆xϕnj pujn`1 ´ unj q ` ∆tϕnj pgj`1{2
n n
´ gj´1{2 q “ 0.
nPN
jPZ

On a :
ÿ ÿ ÿ
∆xϕnj unj “ ∆xϕjn`1 ujn`1 ` ∆xu0j ϕ0j ,
nPN nPN jPZ
jPZ jPZ
ÿ ÿ (2.15)
∆tϕnj gj´1{2
n
“ ∆tϕnj`1 gj`1{2
n
,
nPN nPN
jPZ jPZ

de sorte qu’une intégration par parties discrète donne:


ÿ ÿ ÿ
∆xujn`1 pϕnj ´ ϕjn`1 q ´ ∆xu0j ϕ0j ` n
∆tgj`1{2 pϕnj ´ ϕnj`1 q “ 0 .
nPN jPZ nPN
jPZ jPZ

En introduisant la fonction constante par morceaux gk px, tq “ gj`1{2 n


pour xj ă x ă xj`1 ,
t ă t ă t , l’égalité précédente se réécrit:
n n`1

ż ż
ϕk px, tq ´ ϕk px, t ` ∆tq
uk px, t ` ∆tq dxdt ´ uk px, 0qϕk px, 0qdx ,
RˆR` ∆t R
ż (2.16)
ϕk px ´ ∆x{2, tq ´ ϕk px ` ∆x{2, tq
` gk px, tq dxdt “ 0 .
RˆR` ∆x

39
Le premier terme de cette égalité se réécrit
ż
ϕk px, t ´ ∆tq ´ ϕk px, tq
uk px, tq dxdt
Rˆr∆t,`8s ∆t
ż ˆ ˙
ϕk px, t ´ ∆tq ´ ϕk px, tq
“ uk px, tq 1r∆t,`8s ` Bt ϕpx, tq dxdt
RˆR` ∆t
ż ż
´ puk px, tq ´ upx, tqq Bt ϕpx, tqdxdt ´ upx, tqBt ϕpx, tqdxdt .
RˆR` RˆR`

En utilisant u bornée et le point iiq du deuxième Lemme, on obtient la convergence du


premier terme du membre de droite vers 0. Concernant le second terme:
ˇż ˇ
ˇ ˇ
ˇ pu px, tq ´ upx, tqq B ϕpx, tqdxdtˇ ď }Bt ϕ}L8 pKq }uk ´ u}L1 pKq ÝÑ 0 ,
ˇ `
k t ˇ
RˆR

où l’onż a noté K “ Supppϕq. On a donc établi la convergence du premier terme de (2.16)


vers ´ upx, tqBt ϕpx, tqdxdt.
RˆR` ż ż
Par des arguments analogues, on montre que uk px, 0qϕk px, 0qdx ÝÑ u0 pxqϕpx, 0qdx.
R R

Concernant le troisième terme de (2.16), on remarque que gk px, tq “ gpuk px`∆x{2, tq, uk px´
∆x{2, tqq. En utilisant la consistance du flux, et en supposant la fonction de flux g Lipschitzi-
enne par rapport à chacune de ses variables, il existe une constante C ą 0 (indépendante de
u et de k) telle que:
` ´
|gk px, tq ´ gpu k px, tq, uk px, tqq | ď C|uk px, tq ´ uk px, tq| ` C|uk px, tq ´ uk px, tq| .
loooooooooomoooooooooon (2.17)
f puk px,tqq

On a donc:
ż
ϕk px ´ ∆x{2, tq ´ ϕk px ` ∆x{2, tq
gk px, tq dxdt
RˆR` ∆x
ż
ϕk px ´ ∆x{2, tq ´ ϕk px ` ∆x{2, tq
“ f puk px, tqq dxdt
RˆR` ∆x
ż
ϕk px ´ ∆x{2, tq ´ ϕk px ` ∆x{2, tq
` pgk px, tq ´ f puk px, tqqq dxdt .
RˆR` ∆x
ż
On montre alors que le premier terme converge vers ´ f pupx, tqqBx ϕpx, tqdxdt et le
RˆR`
second vers 0 en utilisant (2.17) avec la convergence uniforme de u˘ k vers u. Au final, (2.16)
donne:
ż ´ ¯ ż
upx, tqBt ϕpx, tq ` f pupx, tqqBx ϕpx, tq dxdt ` u0 pxqϕpx, 0qdx “ 0 ,
RˆR` R

qui n’est autre que la forme faible de l’équation pour u.


Exercice 2.3.1. Les schémas centré, de Lax-Friedrichs et de Lax-Wendroff se réécrivent
sous forme de schémas conservatifs. Trouver ces écritures.
40
2.4 Schéma de Godunov
Définition 2.4.1: Schéma de Godunov
Le schéma de Godunov est donné par
∆t ` n ˘
ujn`1 “ unj ´ n
gj`1{2 ´ gj´1{2 , (2.18)
∆x
avec les flux d’interface gj`1{2
n
“ gpunj , unj`1q définis par la fonction gpu, vq telle que
gpuL , uR q “ f pu˚q, où u˚ “ up0, tq avec upx, tq la solution entropique exacte du prob-
lème de Riemann entre uL et uR .

Remarque 2.4.1.
La solution du problème de Riemann est auto-similaire : upx, tq “ vpx{tq, donc up0, tq “ vp0q
est indépendant de t. Si uL “ uR , alors u˚ “ uL “ uR donc le schéma de Godunov est
consistant. La fonction ξ ÞÑ f pvpξqq est toujours continue en 0 (exercice).
Voir examen.

Ce schéma peut se définir d’une autre manière, qui se décompose en deux étapes:

i) On résout des problèmes de Riemann indépendants à chaque interface xj`1{2 , avec la


donnée initiale unj si x ă xj`1{2 et unj`1 si x ą xj`1{2 .

ii) Après un pas de temps ∆t (à préciser), on moyenne les solutions exactes sur chaque
maille pour obtenir une nouvelle donnée initiale constante par morceaux ujn`1 .

∆t

b b b

xj´1{2 xj xj`1{2 x

Figure 3: Construction du solveur de Godunov : problème de Riemann en xj´1{2 et xj`1{2 .

On choisit ∆t de telle manière que les problèmes de Riemann n’interagissent pas: si


maxjPZ |unj | ď M alors les vitesses d’onde σ sont plus petites que c :“ max|u|ďM |f 1 puq| où u
est solution exacte. On suppose donc la condition CFL
∆x
∆t ď . (2.19)
2c

41
Sous cette condition, la solution dans la maille (xj´1{2 , xj`1{2) est
"
uj´1{2 px, tq si x ă xj ,
upx, tq “
uj`1{2 px, tq si x ą xj ,
avec uj´1{2 px, tq solution du problème de Riemann en xj´1{2 entre unj´1 et unj et uj`1{2px, tq
est solution du problème de Riemann en xj`1{2 entre unj et unj`1. Ainsi, u étant solution
exacte, on a:
ż tn`1 ż xj`1{2 ż xj`1{2 ż xj`1{2
n`1
pBt u ` Bx f puqq dxdt “ 0 ô upx, t qdx ´ upx, tn qdx
tn xj´1{2 xj´1{2 xj´1{2
ż tn`1 ż tn`1
` f pupx, tqq|xj`1{2 dt ´ f pupx, tqq|xj´1{2 dt “ 0 .
tn tn
` ˘
Or, on a f pupx, tqq|xj`1{2 “ f uj`1{2pxj`1{2 , tq , et cette quantité ne dépend pas de t à
ż
1 xj`1{2
cause du caractère auto-similaire de la solution. Ainsi, en posant uj “ n
upx, tn qdx,
∆x xj´1{2
l’égalité précédente se réécrit:
∆t ` ˘
ujn`1 ´ unj ` f pu˚j`1{2q ´ f pu˚j´1{2q “ 0 ,
∆x
et on retrouve le schéma (2.3) avec gpunj , unj`1q “ f pu˚j`1{2q où u˚j`1{2 est la solution du prob-
lème de Riemann en xj`1{2 entre unj et unj`1.

Remarque 2.4.2.
‚ D’un point de vue pratique, on programme la première définition du schéma.
‚ Supposons f convexe. On rappelle que:
- Si uR ă uL , la solution exacte du problème de Riemann est une discontinuité se
f puRq ´ f puL q
propageant à la vitesse σ “ .
uR ´ uL
La solution du problème de Riemann en 0 est donc u˚p0, uL , uR q “ uL si σ ą 0
(c’est à dire si f puR q´f puL q est de signe opposé à uR ´uL , autrement dit f puL q ą
f puR q), et u˚p0, uL , uR q “ uR si σ ă 0 (c’est à dire f puL q ă f puR q).
- Si uR ą uL , la solution est une onde de détente définie par:
$
& uL si x ă σL t
upx, tq “ pf 1 q´1 px{tq si σL t ă x ă σR t
%
uR si σR t ă x
où σL “ f 1 puL q et σR “ f 1 puR q.
On en déduit:
$
& uL si 0 ă f 1 puL q
˚
u p0, uL , uR q “ pf 1 q´1 p0q si f 1 puL q ă 0 ă f 1 puR q
%
uR si f 1 puR q ă 0
42
Par conséquent:
$ #

’ uL ą uR (choc) et f puL q ą f puRq pi.e. σ ą 0q

’ f pu q si


L uL ă uR (détente) etf 1 puL q ą 0
& #
gpuL , uR q “ uL ą uR (choc) et f puL q ă f puR q pi.e. σ ă 0q

’ f pu q si


R uL ă uR (détente) etf 1 puR q ă 0


%
f ˝ pf 1 q´1 p0q si uL ă uR (détente) et f 1 puL q ă 0 ă f 1 puR q
"
˚ uL si f 1 ě 0
En particulier, si f est de plus monotone, on a u p0, uL , uR q “ , et le
uR si f 1 ď 0
schéma de Godunov se confond avec le schéma upwind.
Le schéma de Godunov est difficile à implémenter en pratique, car il requiert la connais-
sance de la solution exacte du problème de Riemann à chaque pas de temps et à chaque
interface du maillage. Ceci motive l’introduction de méthodes plus simples permettant
d’approcher la solution du problème de Riemann.

2.5 Variantes du schéma de Godunov


2.5.1 Schéma de Roe-Murman.
C’est une simplification du schéma de Godunov par linéarisation. On introduit:
f puq ´ f pvq
apu, vq “ si u ‰ v , apu, uq “ f 1 puq . (2.20)
u´v
On remplace le "vrai" problème de Riemann par sa version linéarisée, c’est à dire:
$
& Bt u ` apuR", uL qBx u “ 0 ,
uL si x ă 0 , (2.21)
% upx, 0q “
uR si x ą 0 .
Il s’agit d’une
" équation de transport dont la solution est très simple:
uL si x ă apuR , uL qt ,
upx, 0q “ On en déduit le schéma suivant:
uR si x ą apuR , uL qt .
∆t ` ˘
ujn`1 “ unj ´ gRoe punj , unj`1q ´ gRoe punj´1, unj q ,
∆x
"
f puL q si apuR , uL q ą 0 ,
avec gRoe puL , uR q “ Ainsi, si la solution du problème de Riemann
f puR q si apuR , uL q ă 0 .
est un choc, le schéma de Roe donne la solution exacte. Par contre, sans correction, le schéma
de Roe peut converger vers des solutions non physiques (c’est à dire non entropiques).
Par exemple, supposons f strictement convexe, et uL ă uR (la solution exacte est donc une
détente) et f puL q “ f puR q. Cependant apuL , uR q “ 0, donc gRoe “ f puL q “ f puR q, et par
suite les quantités gj`1{2
n n
´ gj´1{2 sont nulles pour tout n P N et tout j P Z: la solution
numérique est un choc stationnaire.

43
Remarque 2.5.1. A noter qu’avec (2.20), on a:
1 uL ´ uR
gRoe puL , uR q “ pf puR q ` f puL qq ` |apuL , uR q| ,
2
loooooooooomoooooooooon 2
looooooooooomooooooooooon
partie centrée viscosité numérique

et le schéma s’écrit:
∆t ” 1 ` ˘ 1 1 n ı
ujn`1 “ unj ´ f punj`1q ` f punj´1q ´ |anj`1{2 | pu n
´ u n
looooomooooonq ` |a | pu n
´ u n
q .
∆x 2 2 j`1 j
2 j´1{2 looooomooooon
j j´1

∆xBx uj`1{2 ∆xBx uj´1{2

On reconnaît une partie centrée d’ordre 2 en espace (instable) et une viscosité numérique
portée par Bx p|a|Bx uq de l’ordre de ∆x. Lorsque |a| “ 0, le schéma de Roe n’a pas de
viscosité numérique et il peut ne pas converger, ou bien converger vers une solution non
entropique. Pour corriger ce problème, on peut remplacer |a| par une version régularisée aε
qui ne s’annule pas.

2.5.2 Schéma de Enquist-Osher.


Ce schéma s’écrit sous la forme suivante:
ż n ż unj
∆t ´ ¯ ∆t ´ uj`1 1 ¯
n`1 n
uj “ uj ´ n n
f puj`1q ´ f puj´1q ` |f pξq|dξ ´ |f 1 pξq|dξ .
2∆x 2∆x unj un
j´1

Il s’agit d’un schéma conservatif dont le flux numérique est donné par:
ˆ żv ˙
1 1
gpu, vq “ f puq ` f pvq ´ |f pξq|dξ . (2.22)
2 u

Dans les domaines en u où le signe de f 1 est constant, on retrouve le schéma upwind.

2.5.3 Solveurs de Riemann approchés


On cherche à construire une solution approchée auto-similaire du problème de Riemann entre
deux états uL et uR , que l’on notera uppx, tq “ uppx{t, uL , uR q. On demande à ces solveurs
approchés de satisfaire certaines propriétés de consistance:
Définition 2.5.1: Consistance intégrale

Un solveur de Riemann approché u ppx{t, uL , uR q est dit consistant avec la forme


intégrale de la loi de conservation si la relation suivante est vérifiée:
ż ∆x{2
1 1 ∆t
ppx{∆t, uL , uR qdx “ puL ` uR q ´
u pf puR q ´ f puL qq , (2.23)
∆x ´∆x{2 2 ∆x

∆t 1
où ∆t est déterminé via la condition de non interaction max|u|ďM |f 1 puq| ď , avec
∆x 2
M “ maxjPZ |unj |.

44
Cette relation garantit que le schéma ainsi obtenu est conservatif et consistant (au sens
de la Définition 2.3). Plus précisément, on a le résultat suivant:
Proposition 2.5.1

On considère u ppx{t, uL , uRq un solveur de Riemann approché consistant avec la forme


intégrale de la loi de conservation. Alors le schéma:
ż ˆ ˙ ż ˆ ˙
n`1 1 xj x ´ xj´1{2 n n 1 xj`1{2 x ´ xj`1{2 n n
uj “ up , uj´1, uj dx ` p
u , uj , uj`1 dx
∆x xj´1{2 ∆t ∆x xj ∆t

est conservatif et consistant. Plus précisément, on a:


∆t ` n ˘
ujn`1 “ unj ´ n
gj`1{2 ´ gj´1{2 ,
∆x
avec
ż xj`1{2 ˆ ˙
n ∆x n 1 x ´ xj`1{2 n n
gj`1{2 “ f punj q ` uj ´ p
u , uj , uj`1 dx
2∆t xj ∆t ∆t
ż ˆ ˙
n ∆x n 1 xj`1 x ´ xj`1{2 n n
“ f puj`1q ´ u ` p
u , uj , uj`1 dx .
2∆t j`1 ∆t xj`1{2 ∆t

Exemple 2.5.1. Schéma de Godunov.


Le schéma de Godunov vérifie la consistance intégrale. En effet, que l’on soit en présence
d’un choc ou d’une détente, la solution exacte uR du problème de Riemann entre deux états
uL , uR vérifie: "
uL si x{t ă σ ´ ,
uR px{t, uL , uR q “
uR si x{t ą σ ` ,
avec σ ´ ď σ ` la vitesse des ondes. uR étant solution exacte, on a:
ż ∆t ż ∆x{2
rBt uR ` Bx f puRqs dtdx “ 0 .
0 ´∆x{2

Le calcul donne:

ż ∆t ż ∆x{2 ż ∆x{2
Bt uR px, t, uL , uR qdtdx “ puR px, ∆t, uL , uR q ´ u0 pxqq dx
0 ´∆x{2 ´∆x{2
ż ∆x{2
∆x
“ uRpx{∆t, uL , uR qdx ´ puL ` uR q .
´∆x{2 2
ż ∆t ż ∆x{2
Bx f puR px, ∆t, uL , uR qqdt
0 ´∆x{2
ż ∆t ” ı
“ f puRp∆x{2, t, uL , uR qq ´ f puR p´∆x{2, t, uL , uRqq dt .
0
45
σ´ t σ`
∆t

uL uR

∆tσ ´ ∆tσ ` x x
´∆x{2 ∆x{2

Figure 4: Intégration du solveur sur r´∆x{2, ∆x{2s.

∆t
Rappelons que max |f 1 puq| ď 1{2, et par conséquent ∆x{2 ě ∆tσ ` ě tσ ` pour tout
∆x
t P r0, ∆ts (voir Figure 4). Il vient uR p∆x{2, t, uL , uR q “ uR pour tout t P r0, ∆ts, et de même
uR p´∆x{2, t, uL , uR q “ uL pour tout t P r0, ∆ts. Ces quantités sont donc indépendantes de
t, et la quantité précédente s’écrit:
ż ∆t ż ∆x{2
Bx f puRpx, ∆t, uL , uR qqdt “ ∆t pf puR q ´ f puL qq .
0 ´∆x{2

On en tire:
ż ∆x{2
∆x
uR px, ∆t, uL , uR qdt “ puL ` uR q ´ ∆t pf puRq ´ f puL qq ,
´∆x{2 2

et on retrouve bien (2.23).


Exemple 2.5.2. Schéma de Rusanov. On pose λ “ maxminpuL ,uD qďvďmaxpuL ,uD q |f 1 pvq|, et on
introduit le solveur approché:
$
& uL si x{t ă ´λ ,
pR px{t, uL , uR q “
u u˚ si ´ λ ă x{t ă λ ,
%
uR si x{t ą λ .
On choisit l’état intermédiaire u˚ de sorte que u pR satisfasse la relation de consistance. En
se basant sur la Figure (4) avec σ “ ˘λ, et l’état intermédiaire u˚ dans le cône central
˘

séparant uL et uR , on a:
ż ∆x{2 ˆ ˙ ˆ ˙
∆x ˚ ∆x
pR px, ∆t, uL , uR qdx “
u ´ λ∆t uL ` 2λ∆tu ` ´ λ∆t uR
´∆x{2 2 2
∆x
“ puL ` uR q ` λ∆t p2u˚ ´ puL ` uR qq .
2
La définition (2.23) impose:
uL ` uR 1
λ p2u˚ ´ puL ` uR qq “ ´ pf puR q ´ f puL qq ô u˚ “ ´ pf puR q ´ f puL qq .
2 2λ
46
Ainsi, avec cette définition de u˚, le schéma est conservatif et consistant d’après la Proposition
(2.5.3), et s’écrit:
∆t ` n ˘
ujn`1 “ unj ´ n
gj`1{2 ´ gj´1{2 ,
∆x
ż xj`1{2 ˆ ˙
∆t 1 x ´ x j`1{2 n n
avec gj`1{2
n
“ f punj q ` un ´ pR
u , uj , uj`1 dx. On pourra vérifier,
2∆x j ∆t xj ∆t
à titre d’exercice, que:
ˆ n ˙
n 1` n n
˘ uj`1 ´ unj
gj`1{2 “ f puj`1q ` f puj q ´ λ .
2 2

Nous avons vu que la consistance seule ne suffit pas à garantir la convergence des schémas.
De plus, même en cas de convergence, certains schémas peuvent converger vers une solution
non physique (non entropique). Ceci nous incite à dégager un ensemble de critères permettant
de garantir de "bonnes" propriétés de convergence, que nous allons préciser maintenant. A
cet effet, les notions de monotonie et de schéma TVD jouent un rôle central dans l’analyse
des schémas numériques.

2.6 Schémas monotones et entropiques


2.6.1 Schémas monotones - schémas TVD
Définition 2.6.1
´ ¯
Un schéma explicite du type ujn`1 “ Hpunj1 q est dit monotone si la
j´Kďj 1 ďj`K
fonction H est croissante par rapport à chacun des ses arguments.

Proposition 2.6.1
´ ¯
n`1
Un schéma uj “ Hpuj q est monotone si et seulement si:
n

´ ¯ ´ ¯
u0j ě vj0 , @j P Z ùñ unj ě vjn , @j P Z , n P N . (2.24)

Proof. Le sens direct est immédiat puisque u1j “ Hpu0j q ě Hpvj0 q “ vj1 pour tout j P Z et on
conclut par récurrence.
Considérons un schéma explicite ujn`1 “ Hpunj´K , ¨ ¨ ¨ , unj`K q qui satisfait le principe de
comparaison (2.24). Fixons alors j P Z et 2K réels aj´K`1, ¨ ¨ ¨ , aj`K . Considérons deux
réels u, v tels que u ě v, et des données initiales discrètes pu0 q et pv 0 q telles que:

u0j 1 “ aj 1 si j ´ K ă j 1 ď j ` K , uj 1 “ u si j 1 “ j ´ K ,
(2.25)
vj01 “ aj 1 si j ´ K ă j 1 ď j ` K , vj 1 “ v si j 1 “ j ´ K .

47
Le principe de comparaison donne :

u1j “ Hpu, aj´K`1, ¨ ¨ ¨ , aj`K q ě Hpv, aj´K`1, ¨ ¨ ¨ , aj`K q “ vj1 , (2.26)

ce qui permet d’établir que H est croissante par rapport à sa première variable. Un raison-
nement analogue donne la croissance de H par rapport à toutes ses variables.

Proposition 2.6.2
On considère un schéma explicite conservatif à 3 points:
∆t ` n n ˘
ujn`1 “ unj ´ gpuj , uj`1q ´ gpunj´1, gjn q ,
∆x
avec g de classe C 1 . S’il est monotone alors u ÞÑ gpu, vq est croissante et v ÞÑ gpu, vq
est décroissante.
Réciproquement, si u ÞÑ gpu, vq est croissante et v ÞÑ gpu, vq est décroissante, alors
sous condition CFL le schéma est monotone.

∆t
Proof. On a Hpu, v, wq “ v ´ λ pgpv, wq ´ gpu, vqq avec λ “ . Notons tout d’abord que
∆x
Bu Hpu, v, wq “ λB1 gpu, vq et Bw Hpu, v, wq “ ´λB2 gpv, wq, de sorte que:
ˇ
ˇ u ÞÑ Hpu, v, wq varie comme u ÞÑ gpu, vq ,
ˇ
ˇ w ÞÑ Hpu, v, wq varie comme w ÞÑ ´gpv, wq ,

ce qui permet d’établir le sens direct.


Supposons à présent u ÞÑ gpu, vq croissante et v ÞÑ gpu, vq décroissante. D’après ce qui
précède, on a la croissance de H par rapport à u et w. D’autre part:

Bv Hpu, v, wq “ 1 ´ λ pB1 gpv, wq ´ B2 gpu, vqq . (2.27)

En supposant u, v, w P K Ă R avec K compact, on a |B1 gpv, wq´B2gpu, vq| ď M avec M ą 0.


∆t
Donc, sous la CFL ď 1{M on a le résultat.
∆t

Proposition 2.6.3
Si un schéma conservatif est monotone, alors l’opérateur H vérifie le principe du
maximum discret:

min punk q ď unj ď max punk q .


kPtj´1,j,j`1u kPtj´1,j,j`1u

Exemple 2.6.1.

48
‚ Schéma de Lax-Friedrichs:
` ˘
ujn`1 “ ujn`1 ´ λ gpunj , unj`1q ´ gpunj´1, gjn q ,
1 1
avec gpu, vq “ pf puq ` f pvqq ´ pv ´ uq. On a:
2 2λ
1 1 ∆x
Bu g “ f 1 puq ` λ ě 0 si ∆t ď avec K compact, (2.28)
2 2 maxuPK |f 1 puq|
1 1
et Bv g “ f 1 pvq ´ λ ď 0 sous la même condition, ce qui donne la monotonie sous
2 2
condition CFL.
‚ Schéma de Lax-Wendroff :
Ce schéma n’est pas monotone. Considérons ˇ par exemple l’équation Bt u ´ Bx u “ 0,
ˇ 0 si j ď 0,
avec la condition initiale discrète u0j “ ˇˇ Nous avons vu que le schéma est
1 si j ě 1 ,
donné par:
λ` n ˘ λ2 ` n ˘
ujn`1 “ unj ` uj`1 ´ unj´1 ` uj`1 ´ 2unj ` unj´1 ,
2 2
stable sous la condition 0 ď λ ď 1. Or, le calcul donne:
λ λ2 λ
u11 “ 1 ` ` p1 ´ 2 ` 0q “ 1 ` p1 ´ λq ą 1 si λ ă 1 . (2.29)
2 2 2
Le principe du maximum discret n’est pas respecté, donc le schéma n’est pas monotone.
‚ Schéma d’Enquist-Osher :
La fonction de flux numérique associée au schéma s’écrit:
ˆ żv ˙
1 1
gpu, vq “ f puq ` f pvq ´ |f pξq|dξ , (2.30)
2 u

1 1
d’où Bu g “ pf 1 puq ` |f 1 puq|q ě 0 et Bu g “ pf 1 pvq ` |f 1 pvq|q ě 0 et le schéma est
2 2
monotone.
‚ Schéma de Godunov :
On montre que ce schéma vérifie le principe de comparaison discret (2.24). Rappelons que
ce schéma peut être construit en deux étapes:
i) Résolution exacte de problèmes de Riemann à chaque interface du maillage.
ii) Moyennisation des solutions sur chaque maille.
Considérons deux données initiales pvj0 qjPZ et pu0j qjPZ telles que u0j ě vj0 pour tout j P Z. Rap-
pelons que la solution u associée à la condition initiale pu0j qjPZ dans la maille (xj´1{2 , xj`1{2 )
est "
uR p0, u0j´1, u0j q si x ă xj ,
uj px, tq “
uR p0, u0j , u0j`1q si x ą xj ,
49
ˆ ˙
x ´ xj´1{2 0
avec uR la solution exacte du problème de Riemann entre les états
, uj´1, u0j
t ˆ ˙
x ´ xj`1{2 0 0
uj´1 et uj au niveau de la maille xj´1{2 et uR
0 0
, uj , uj`1 la solution exacte du
t
problème de Riemann entre les états u0j et u0j`1 au niveau de la maille xj`1{2 . On définit
de même vj px, tq la solution associée à la condition initiale pvj0 qjPZ. D’après le principe de
comparaison des solution exactes, la condition u0j ě vj0 pour tout j P Z implique
ˆ ˙ ˆ ˙
x ´ xj´1{2 0 0 x ´ xj´1{2 0 0
uR , uj´1, uj ě uR , vj´1 , vj ,
t t

avec une inégalité similaire à l’interface xj`1{2 , et par suite uj px, tq ě vj px, tq. Par suite:
ż xj`1{2 ż xj`1{2
n`1
uj “ uj px, tqdx ě vj px, tqdx “ vjn`1 ,
xj´1{2 xj´1{2

de sorte que le principe de comparaison discret (2.24) est satisfait.


´ ¯
Introduisons à présent quelques notations. Pour une suite discrète unj nPN, on définit:
jPZ

u∆ px, tq “ unj si tn ď t ă tn`1 et xj´1{2 ď x ă xj`1{2 .

On introduit les normes et semi-normes discrètes suivantes:


ÿ
}u∆ p., tn q}L1 pRq “ ∆x|unj | ,
jPZ

}u∆ p., t q}L8 pRq “ max |unj | ,


n
jPZ
ÿ
T V pu∆ p., t qq “ |unj`1 ´ unj | .
n

jPZ

On notera parfois formellement un “ u∆ p., tn q.


Définition 2.6.2: Schéma TVD
Un schéma est dit TVD (Variation Totale Décroissante) si on a T V pun`1 q ď
T V pun q pour tout n P N. Un schéma est dit Stable L8 s’il existe c ą 0 indépendant
de n et de ∆t tel que }u∆ p., tn q}L8 pRq ď c.

Nous allons voir que de tels schémas sont bien adaptés au calcul de solutions entropiques.
Avant d’aller plus loin, nous aurons besoin du lemme suivant.

50
Lemme 2.6.1: Crandall-Tartar-Majda

Soit (Ω, µ) un espace mesurable et C un sous-espace de L1 pΩq tel que pf, g P Cq ùñ


maxpf, gq P C. Soit
T : C ÝÑ L1 pΩq
f ÞÝÑ T pf q
ż ż
un opérateur vérifiant T pf qdµ “ f dµ pour tout f P C. Alors les trois propriétés
Ω Ω
suivantes sont équivalentes:

i) @f, g P C , f ě g ùñ T pf q ě T pgq ,
ż ż
ii) pT pf q ´ T pgqq` dµ ď pf ´ gq` dµ ,
Ω Ω
ż ż
iii) |T pf q ´ T pgq| dµ ď pf ´ gq dµ .
Ω Ω

Proof. iq ñ iiq : On écrit maxpf, gq “ g ` pf ´ gq` ě g. Par suite, T pmaxpf, gqq ě T pgq,
et de même T pmaxpf, gqq ě T pf q et par conséquent T pmaxpf, gqq ě maxpT pf q, T pgqq. On
retranche alors T pgq:

T pmaxpf, gqq ´ T pgq ě max pT pf q, T pgqq ´ T pgq “ pT pf q ´ T pgqq` .

On intègre: ż ż
pT pf q ´ T pgqq` dµ ď rT pmaxpf, gqq ´ T pgqs dµ ,
Ω Ω
ż ż ż
et on conclut avec pT pmaxpf, gqq ´ T pgqq dµ “ maxpf, gq ´ gdµ car T est conservatif,
ż żΩ ż Ω Ω

et maxpf, gq ´ gdµ “ pf ´ gq` .


Ω Ω Ω
iiq ñ iiiq : On utilise |f ´ g| “ pf ´ gq` ` pg ´ f q` . D’où:
ż ż ż
|T pf q ´ T pgq|dµ “ pT pf q ´ T pgqq` dµ ` pT pgq ´ T pf qq` dµ
Ω żΩ ż Ω ż
ď pf ´ gq` dµ ` pg ´ f q` dµ “ |f ´ g|dµ .
Ω Ω Ω

1
iiiq ñ iiq : On utilise la relation pf ´ gq` “ p|f ´ g| ` pf ´ gqq et la conservativité de
2
l’opérateur T :
ż ˆż ż ˙
1
pT pf q ´ T pgqq` dµ “ |T pf q ´ T pgq|dµ ` pT pf q ´ T pgqq dµ
Ω 2 Ω Ω
ˆż ż ˙ ż
1
ď |f ´ g|dµ ` pf ´ gq dµ “ pf ´ gq` dµ .
2 Ω Ω Ω

51
ż
iiq ñ iq : Supposons f ě g, c’est à dire pg ´ f q` dµ “ 0. On a donc:

ż
0ď pT pgq ´ T pf qq` dµ ď 0 ùñ T pgq ď T pf q .

En particulier, pour un opérateur qui conserve la masse dans L1 pωq, la monotonie est
équivalente au caractère contractant L1 . Munis de ce Lemme, on établit un premier résultat
important.
Théorème 2.6.1
Soit un schéma explicite conservatif et monotone. Alors il est TVD et stable L8 , et il
vérifie:

i) }u∆ p., tq}L8 pRq ď }u∆ p., 0q}L8 pRq ,

ii) }u∆ p., tq}L1 pRq ď }u∆ p., 0q}L1pRq ,

iii) T V pu∆ p., tqq ď T V pu∆ p., 0qq ,

iv) }u∆ p., tq ´ v∆ p., tq}L1 pRq ď }u∆ p., 0q ´ v∆ p., 0q}L1 pRq ,

v) Si g est Lipschitz-continue, alors:

‚ }u∆ p., tn q ´ u∆ p., tm q}L1 pRq ď c|n ´ m|∆tT V pu∆ p., 0qq ,
‚ }u∆ p., tq ´ u∆ p., sq}L1 pRq ď c p|t ´ s| ` ∆tq T V pu∆ p., 0qq .

Proof. iq : Il s’agit d’établir que

min u0j ď unj ď max u0j , @j P Z , @n P N .


jPZ jPZ

On note formellement ujn`1 “ Hj pun q “ Hpuj´K , ¨ ¨ ¨ , uj`K q avec H monotone et conservatif.


Notons M “ maxjPZ u0j . Considérons la solution initiale définie par vj0 “ M pour tout j P Z
On a: ¨ ˛
` 0 ˘
Hj pMq “ M , @j P Z‚ ,
uj ď vj0 , @j P Z ñ ˝u1j “ Hj pun q ď loooooomoooooon (2.31)
Hconservatif

et on montre de même que u est minorée par m “ minjPZ u0j . On en déduit la stabilité L8 .
1

iiq : Ce point est une conséquence immédiate de (iv) avec v∆ “ 0.


iiiq : On montre que le schéma est T V D. Soit t tel que tn ď t ă tn`1 . On a:
ÿ 1
T V pu∆ p., tqq “ T V pu∆ p., tn qq “ |unj`1 ´ unj | “ }u∆ p. ` ∆x, tn q ´ u∆ p., tn q}L1 pRq .
jPZ
∆x

52
On applique alors (iv) en t “ tn :
1 1
}u∆ p. ` ∆x, tn q ´ u∆ p., tn q}L1 pRq ď }u∆ p. ` ∆x, 0q ´ u∆ p., 0q}L1pRq “ T V pu∆ p., 0qq
∆x ∆x
ivq : On utilise le Lemme précédent dans L1 pRq avec C “ tf P L1 pRq , constante par mailleu et
ÿ ÿ
l’opérateur T “ H n . En effet, cet opérateur est conservatif: la relation ujn`1 ∆x “ u0j ∆x
jPZ jPZ
pour toute donnée initiale garantie par la conservativité du schéma peut se réécrire
pu0j q
ż ż
n
H pu0 pxqqdx “ u0 pxqdx pour tout u0 P C .
R R

Ainsi, la monotonie de H n (qui se déduit immédiatement de celle de H) implique que H n


est un opérateur contractant dans L1 . Ainsi, pour tout t tel que tn ď t ă tn`1 :

}u∆ p., tq ´ v∆ p., tq}L1 pRq “ }u∆ p., tn q ´ v∆ p., tn q}L1 pRq
“ }H n pu∆ p., 0qq ´ H n pv∆ p., 0qq }L1 pRq
ď }u∆ p., 0q ´ v∆ p., 0q}L1pRq ,

ce qui donne le résultat.


vq : Soit m, n deux entiers avec n ą m. En utilisant (iv), on a:

ÿ
n´1
n m
}u∆ p., t q ´ u∆ p., t q}L1 pRq ď }u∆ p., tp`1 q ´ u∆ p., tp q}L1 pRq
p“m

ď |n ´ m|}u∆ p., ∆tq ´ u∆ p., 0q}L1 pRq .

D’autre part:
ÿ ÿ
0
}u∆ p., ∆tq ´ u∆ p., 0q}L1 pRq “ ∆x| H j pu q
loomoon ´u0j | “ 0
∆t|gj`1{2 0
´ gj´1{2 |.
jPZ jPZ
∆t
u0j ´ ∆x pgj`1{2
0 0
´gj´1{2 q
` ˘
On a gj`1{2
0
“ g u0j , u0j`1 , et en utilisant g Lipschitz-continue:
ÿ ÿ` ˘
0 0
∆t|gj`1{2 ´ gj´1{2 | ď κ∆t |u0j`1 ´ u0j | ` |u0j ´ u0j´1|
jPZ jPZ

ď 2κ∆tT V pu∆ p., 0qq ,

et on conclut le premier point avec c “ 2κ.


Pour s, t quelconques, on considère n, m entiers tels que tn ď s ă tn`1 et tm ď t ă tm`1 .
Il est alors aisé de vérifier que |n ´ m|∆t “ |tn ´ tm | ď |s ´ t| ` ∆t et u∆ p., tq ´ u∆ p., sq “
u∆ p., tn q ´ u∆ p., tm q, d’où le résultat.

53
2.6.2 Schémas entropiques
Définition 2.6.3: Entropie et consistance

Soit la condition d’entropie Bt Epuq ` Bx Gpuq ď 0, avec E une fonction convexe et G


le flux d’entropie associé.

‚ Un schéma explicite du type ujn`1 “ Hpunj´1, unj , unj`1q est dit consistant avec
la condition d’entropie s’il existe une fonction Gpu, vq:

i) consistante avec le flux d’entropie, c’est à dire Gpu, uq “ Gpuq pour tout
u P R.
ii) vérifiant l’inégalité d’entropie discrète:

∆t ` ` n n ˘ ` ˘˘
Epujn`1 q ´ Epunj q ` G uj , uj`1 ´ G unj´1, unj ď 0 .
∆x
La fonction G est appelée flu d’entropie.

‚ Si le schéma est consistant avec toutes les conditions d’entropie, on dit que ce
schéma est entropique.

Théorème 2.6.2
Un schéma explicite, conservatif, consistant et monotone est entropique.

Proof. Nous admettrons la démonstration. Notons qu’en pratique, il suffit de démontrer la


consistance avec les entropies de Kruzkov

Epuq “ |u ´ k| et Gpuq “ sgnpu ´ kq pf puq ´ gpuqq ,

avec k P R.
On déduit de ce théorème que les schémas de Godunov, Lax-Friedrichs, Enquist-Osher
sont entropiques (éventuellement sous CFL). Le schéma de Roe peut-être corrigé de manière
à être entropique.
La notion de consistance intégrale (2.23) introduite pour les solveurs de Riemann ap-
prochés s’étend aux schémas entropiques. Plus précisément:
Définition 2.6.4
On dit qu’un solveur de Riemann approché u ppx{t, ul , uR q vérifie la consistance inté-
grale avec l’inégalité d’entropie si la relation suivante est satisfaite:
ż ∆x{2
1 1 ∆t
E pp
upx{t, ul , uR qq dx ď pEpuL q ` EpuR qq ´ pGpuR q ´ GpuL qq . (2.32)
∆x ´∆x{2 2 ∆x

54
Théorème 2.6.3
Soit uppx{t, ul , uR q un solveur de Riemann approché vérifiant (2.23,2.32), c’est à dire
vérifiant la consistance intégrale de la loi de conservation et la consistance intégrale
avec l’inégalité d’entropie. Ce schéma vérifie l’inégalité d’entropie discrète:
∆t ` ` n n ˘ ` ˘˘
Epujn`1q ´ Epunj q ` G uj , uj`1 ´ G unj´1, unj ď 0 ,
∆x
avec
ˆ ˆ
ż xj`1{2 ˙˙
` ˘ ∆x 1 x ´ xj`1{2
G unj , unj`1 “ Gpunj q ` ηpuni q ´ E u p , uL , uR
2∆t xj∆t ∆t
ż ˆ ˆ ˙˙ (2.33)
n ∆x n 1 xj`1{2 x ´ xj`1{2
“ Gpuj q ´ ηpuj`1q ´ E up , uL , uR
2∆t ∆t xj ∆t

Théorème 2.6.4: Consistance avec les inégalités d’entropie


On considère un schéma explicite conservatif, consistant et entropique. On suppose
∆t
que lorsque ∆x tend vers 0, avec λ “ fixé, la suite de solutions discrètes u∆ vérifie:
∆x
1. }u∆ }L8 pRq ď c

2. u∆ ÝÑ u dans L1loc pR ˆ R` q et p.p. dans R ˆ R` .

Alors la limite u est solution faible entropique de l’équation

Proof. La preuve est calquée sur celle du théorème de Lax-Wendroff : on multiplie l’inégalité
d’entropie discrète par ϕnj , on effectue une intégration par parties discrète et on passe à la
limite.

Théorème 2.6.5: Convergence des schémas

On considère un schéma explicite conservatif, consistant et monotone (il est donc


1 şxj`1{2 0
entropique). Soit u0 pxq P pL1 X L8 X BV q pRq, et soit v 0 pxq “ vj0 “ u pxqdx
∆x xj´1{2
dans chaque maille (xj´1{2 , xj`1{2 q. On considère une suite de pas d’espace ∆x qui
∆t
tend vers 0 avec λ “ fixé et on note u∆ la suite de solutions discrètes associée.
∆x
Alors il existe une solution faible entropique u de donnée initiale u0 telle que, pour
une sous-suite, u∆ Ñ u p.p. et dans L8 pR` , L1loc pRqq

55
Chapter 3

Problème de Riemann pour les systèmes


hyperboliques

3.1 Systèmes hyperboliques unidimensionnels


Ce chapitre vise à établir certains résultats généraux sur les lois de conservation en 1d. On
considère dans cette section des systèmes d’équations aux dérivées partielles du premier ordre
s’écrivant sous la forme suivante :

Bt wpx, tq ` Bx pF pwpx, tqqq “ 0 , (3.1)

où w : RˆR` Ñ Rp est le vecteur d’état, et F : Rp Ñ Rp le flux correspondant. Le problème


de Cauchy associé à (3.1) s’écrit:
"
Bt wpx, tq ` Bx pF pwpx, tqqq “ 0 , x P R, t ą 0 , (3.2a)
wpx, 0q “ w0 px, 0q , x P R . (3.2b)

Dans la suite, on supposera que la solution vit dans un ensemble convexe d’état admissibles,
noté Ω.

3.1.1 Hyperbolicité, symétrisabilité et solutions fortes


Définition 3.1.1: Solutions fortes
On appelle solution forte (ou solution classique) du problème de Cauchy (3.2) sur
R ˆ r0, T s une fonction w P C 1 pR ˆ r0, T sq qui vérifie (3.2a) et (3.2b).

Remarque 3.1.1. Si w est une solution forte du système (3.1), alors c’est aussi une solution
du système suivant, appelé forme quasi-linéaire ou forme non conservative du système
du lois de conservation:
Bt w ` ApwqBx w “ 0, (3.3)
où Apwq “ ∇w F pwq est la matrice jacobienne de l’application w ÞÑ F pwq.

56
Une notion importante pour la construction de solutions fortes (comme de solutions
faibles) au problème de Cauchy (3.2) est la notion d’hyperbolicité.
Définition 3.1.2: Systèmes hyperboliques

Le système (3.1) est dit hyperbolique si la matrice jacobienne Apwq “ ∇w F pwq est
diagonalisable et à valeurs propres réelles pour tout w P Ω. Si les valeurs propres sont
toutes distinctes, le système est dit strictement hyperbolique

Dans la suite, nous noterons ces valeurs propres pλi pwqq1ďiďp , et les ordonnerons par
ordre croissant:
λ1 pwq ď λ2 pwq ď ¨ ¨ ¨ ď λp pwq . (3.4)

Remarque 3.1.2. L’hyperbolicité garantit la stabilité linéaire des constantes. Un état


constant w˚ est solution de (3.1). Cette solution est dite linéairement stable si toute per-
turbation w0 P L2 pRq de cette donnée initiale constante conduit à une solution bornée dans
L2 pRq, uniformément en temps, pour le système linéarisé autour de w˚ :
"
Bt w ` Apw˚ qBx w “ 0 , x P R, t ą 0 , (3.5a)
wpx, 0q “ w0 px, 0q , x P R. (3.5b)

Autrement dit, Apw˚ q est R-diagonalisable si et seulement si (3.5) est bien posé dans L2 ,
c’est-à-dire si et seulement si pour tout w0 P L2 pRq, (3.5) admet une unique solution w P
CpR` , L2 pRqq qui dépend continûment de w0 : il existe C ą 0 telle que

sup kwp., tqkL2 pRq ď C kw0 kL2 pRq .


tPR`

Une propriété pratique pour vérifier l’hyperbolicité d’un système est la suivante:
Proposition 3.1.1
La propriété d’hyperbolicité est invariante par changement de variable.

Proof. Soit w “ φpvq un changement de variable admissible. La matrice ∇v φpvqq est donc
inversible pour tout v. La fonction w est solution forte de (3.3), si et seulement si, v est
solution forte de
∇v φpvqBtv ` Apφpvqq∇v φpvqBx v “ 0
ðñ Bt v ` BpvqBx v “ 0,
avec Bpvq “ p∇v φpvqq´1Apφpvqq∇v φpvq. Les matrices Bpvq et Apφpvqq sont semblables.
Ainsi, Apwq est R-diagonalisable pour tout w si et seulement si Bpvq est R-diagonalisable
pour tout v.
En pratique, on s’efforce de privilégier un changement de variable rendant la diagonali-
sation de la matrice Bpvq la plus simple possible.

57
Exemple 3.1.1. En 1d, en regroupant les variables d’écoulement dans la quantité w “
ph, qqT avec q “ hu, les équations de Saint-Venant sur fond plat s’écrivent sous la forme:

Bt w ` Bx pF pwqq “ 0 , (3.6)
¨ ˛
q
avec F pwq “ ˝ q 2 h2 ‚. L’ensemble des états admissibles est le sous-ensemble de R2
`g
h 2
caractérisé par Ω “ tph, qq , h ą 0 , q P Ru. On vérifie que la jacobienne associée au flux est
donnée par: ˆ ˙
0 1
Apwq “ 2 . (3.7)
´ hq 2 ` gh 2 hq
? ?
Les valeurs propres associées sont λ1 pwq “ hq ´ gh, λ2 pwq “ hq ` gh. On en déduit que le
système est hyperbolique, et même strictement hyperbolique lorsque h ą 0.

En fait, l’hyperbolicité d’un système ne suffit pas à garantir l’existence de solutions fortes
même localement en temps. Pour ce faire, on introduit la notion de symétrisabilité:
Définition 3.1.3: Système symétrisable

Le système (3.1) est dit symétrisable si pour tout w P Ω, il existe une matrice Spwq,
symétrique, définie, positive telle que la matrice SpwqApwq soit symétrique.

3.1.2 Solutions faibles


Rappelons qu’une solution w de (3.2) est appelée solution classique (ou solution solution
forte) si w P C 1 pR ˆ R` q et satisfait l’équation (3.2) en tout point.
Comme dans le cas scalaire, les solutions des systèmes hyperboliques peuvent développer des
discontinuités en temps fini, d’où la nécessité donner du sens à la notion de solution pour
des fonctions non régulières. Dans un premier temps, si l’on suppose w P C 1 pR ˆ R` q, on
peut multiplier (3.2) par une fonction test φ P C08 pR ˆ R` q et intégrer sur R ˆ R` . Après
intégration par parties:
ij ż
pwBt φ ` F pwqBx φq dxdt ` w0 φpx, 0qdx “ 0 .
R
RˆR`

58
Notons qu’alors la dérivabilité de w n’est plus requise pour donner sens à la formule précé-
dente. Ceci motive la définition suivante:
Définition 3.1.4: Solution faible
On dit que w P L1loc pR ˆ R` q est solution faible de (3.2, 3.2b) si w vérifie:
ż ż ż
pwBtφ ` F pwqBx φq dxdt ` w0 φpx, 0qdx “ 0 , @φ P C08 pR, R˚` q . (3.8)
RˆR` R

Nous venons de voir que toute solution classique de (3.2) était une solution faible au sens
de (3.8). Réciproquement, si l’on considère une solution faible w (3.8) de classe C 1 pR ˆ R` q,
alors on vérifie aisément que w est une solution classique de (3.2).

Dans le cas d’une discontinuité, nous allons voir à présent que les valeurs de w et F pwq
sont reliées par une relation, appelée relation de Rankine-Hugoniot. Plaçons nous dans
un ouvert Σ Ă RˆR` et supposons que la solution w soit séparée en deux solutions régulières
wG et wD par une courbe de discontinuité paramétrée par une fonction Γ : t ÝÑ pγptq, tq
de classe C 1 pR` q (Figure 1). L’ouvert Σ se partitionne alors en deux ouverts Σ´ et Σ` sur
lesquels évoluent respectivement wG et wD . En tout point px, tq de la courbe Γ, on note:

wG px, tq “ lim wG py, sq et wD px, tq “ lim wD py, sq . (3.9)


py,sqÑpx,tq py,sqÑpx,tq

En premier lieu, il convient de rappeler que w est solution classique en tout point où elle
est C 1 . En d’autres termes, wG est solution classique sur Σ´ et wD est solution classique
sur Σ` . On s’intéresse à présent au comportement de la solution au niveau de la courbe de
discontinuité. w étant une solution faible, on a:
ż
pwBt φ ` F pwqBxφq dxdt “ 0 , @φ P C08 pΣq ,
Σ

ce qui se réécrit:
ż ż
pwBt φ ` F pwqBx φq dxdt ` pwBt φ ` F pwqBxφq dxdt “ 0 , @φ P C08 pΣq . (3.10)
Σ´ Σ`

Après application de la formule de Green 1


ż ż
pwBt φ ` F pwqBxφq dxdt “ ´ pBt w ` Bx F pwqq φ dxdt
Σ´ Σż
´
` ˘
` wG n´ ´
t ` F pwG qnx φ dxdt ,
ż ż Σ´ XΓ (3.11)
pwBt φ ` F pwqBxφq dxdt “ ´ pBt w ` Bx F pwqq φ dxdt
Σ` Σż
`
` ˘
` wD n` `
t ` F pwD qnx φ dxdt ,
Σ` XΓ
1
ş ş ş
upx, tqdivpvpx, tqqdν “ ´
Σ Σ
upx, tqdivpvpx, tqqdν ` BΣ
upx, tqvpx, [Link]σ, appliquée ici avec upx, tq “ φ
T
et vpx, tq “ pw, F pwqq .

59
˘ T
où ~n˘ “ pn˘ ˘
x , nt q désignent les vecteurs normaux extérieurs à Γ relativement à Σ (Figure
˘
1). La solution étant régulière sur Σ , les intégrales de surface disparaissent dans l’équation
précédente, pour donner:
ż ż
` ´ ´
˘ ` ˘
wD nt ` F pwD qnx φ dxdt ` wD n`t ` F pw D qn`
x φ dxdt “ 0 , (3.12)
Σ´ XΓ Σ` XΓ

ou encore ż
` ˘
pwG ´ wD q n´ ´
t ` pF pwG q ´ F pwD qq nx φ dxdt “ 0 . (3.13)
ΣXΓ
Finalement, en utilisant le fait que la direction normale à la courbe Γ est portée par le vecteur
p1, ´γ 1 ptqqT , on obtient les relations de Rankine-Hugoniot:

γ 1 ptq pwD ´ wG q “ F pwD q ´ F pwG q , (3.14)

ou de manière plus compacte:


γ 1 ptqrrwss “ rrF pwqss , (3.15)
La quantité γ 1 est appelée vitesse de propagation de la discontinuité, et usuellement notée
σ dans la littérature. Nous venons d’établir le résultat suivant:
t

Γ
Σ´
wG

~n´
Σ`
wD

Figure 1: Relation de Rankine-Hugoniot : séparation des états le long de la courbe de dis-


continuité γ.

Théorème 3.1.1: Théorème de Rankine Hugoniot

w est une solution faible de (3.1) si et seulement si w est solution classique là où elle
est de classe C 1 , et vérifie la relation

σrrwss “ rrF pwqss , (3.16)

le long de la discontinuité Γ.

60
Remarque 3.1.3. Cas d’un choc à vitesse constante On suppose la courbe de discon-
tinuité caractérisée par γptq “ σt. Dans ce contexte, la solution s’écrit

w “ wG ` pwD ´ wG q Hpx ´ σtq , (3.17)

où H désigne la fonction de Heaviside :


"
0 si x ă 0 ,
Hpxq “ (3.18)
1 si x ą 0 .

Une démonstration rapide consiste à dériver (3.17) par rapport au temps (au sens des dis-
tributions)2 , pour obtenir:

Bt w “ Bt wG ` pBt wD ´ Bt wG q Hpx ´ σtq ´ σ pwD ´ wG q δx´σt . (3.19)

D’autre part :
F pwq “ F pwG q ` pF pwD q ´ F pwG qq Hpx ´ σtq ,
et par un raisonnement analogue:

Bx F pwq “ Bx F pwG q ` pBx F pwD q ´ Bx F pwG qq Hpx ´ σtq ` pF pwD q ´ F pwG qq δx´σt . (3.20)

En sommant (3.19) et (3.20), et en utilisant (3.2), on obtient immédiatement:

pF pwD q ´ F pwG qq δx´σt “ σ pwD ´ wG q δx´σt ,

d’où l’on tire les relations de saut.

3.1.3 Inégalités d’entropie


Comme dans le cas scalaire, la problème de Cauchy peut admettre une infinité de solutions
faibles. Afin de sélectionner les solutions physiquement admissibles (en d’autres termes, les
"bonnes" solutions), on adjoint le problème d’une nouvelle condition, communément appelée
inégalité d’entropie. Nous commençons par la définition générale suivante:
Définition 3.1.5: Entropie

Une fonction convexe η P C 1 pΩq est une entropie du système (3.1) s’il existe une
fonction ψ P C 1 pΩq appelée flux d’entropie telle que

Bt ηpwq ` Bx pψpwqq “ 0 (3.21)

pour toute solution classique w de (3.1) .

2
On rappelle que la dérivée de la fonction H est la distribution de Dirac δ et que si τx,σ désigne le
changement de variables t ÞÑ x ´ σt alors on a pT ˝ τx,σ q1 “ ´σ.T 1 ˝ τx,σ pour toute distribution T .

61
Remarque 3.1.4. Dans le cas des équations de Saint-Venant 1d, l’équation d’énergie (6.26)
est donnée par: ˆ 2 ˙ ˆ ˙
h 1 2 2 1 3
Bt g ` hu ` Bx gh u ` hu “ 0 . (3.22)
2 2 2
h2 1
La fonction définie par ηpwq “ g ` hu2 est une entropie du système, associée au flux
2 2
1 3
ψpwq “ gh u ` hu . La définition précédente s’étend de manière naturelle aux systèmes
2
2
2d. Par exemple, nous avons vu que pour toute fonction régulière w, l’énergie totale ηpwq “
1 2 1
gh ` h}u}2 du système de Saint-Venant 2d satisfaisait l’équation Bt ηpwq`divpψpwqq “ 0,
2 2 ˆ ˙
1
où ψpwq “ gh ` h}u} u. Dans les deux cas évoqués précédemment il reste bien sûr à
2 2
2
s’assurer que η est convexe, ce qui revient à montrer que sa Hessienne est définie positive.

Cette notion s’étend au cas des solutions faibles à travers la définition suivante:
Définition 3.1.6: Solution faible entropique

Une solution faible de (3.1) est dite entropique si pour tout couple entropie/flux pη, ψq
on a l’inégalité d’entropie suivante:
ij ż
pηpwqBt φ ` ψpwqBxφq dxdt ` ηpw0 pxqqφpx, 0qdx ě 0 , @φ P C08 pR, R˚` q , φ ě 0 .
R
RˆR`

Remarque 3.1.5. Cette définition se réécrit au sens faible:

Bt ηpwq ` Bx pψpwqq ď 0 . (3.23)

Ce critère peut se justifier rigoureusement d’un point de vue physique, par le biais
d’approches avec viscosité evanescente.3 Considérant la remarque précédente (3.1.4), ce
critère admet une interprétation simple d’un point de vue physique: il garantit que l’énergie
du système est dissipée à travers les chocs. Dans le cas scalaire, on peut montrer que ce
critère garantit l’unicité des solutions, mais à ce jour ce résultat est seulement conjecturé dans
le cadre général des systèmes hyperboliques. Nous verrons toutefois qu’il garantit l’existence
et l’unicité de la solution entropique du problème de Riemann sous certaines conditions. En
pratique, on utilise le critère de Lax pour caractériser les solutions faibles entropiques:
3
L’idée est de de montrer (3.23) en voyant les solutions faibles de l’EDP comme limites de solutions d’un
système visqueux Bt pwǫ q ` Bx pF pwǫ qq “ ´ǫ∆w.

62
Définition 3.1.7: Conditions d’admissibilité de Lax
Une discontinuité séparant des états wG et wD se propageant à une vitesse σ est une
choc physiquement admissible si et seulement s’il existe un indice k associé à un champ
vraiment non-linéaire tel que:

λk pwG q ě σ ě λk pwD q ,
λj pwG q ď σ et λj pwD q ď σ , @j ă k , (3.24)
λj pwG q ě σ et λj pwD q ě σ , @j ą k .

La démonstration de ce résultat général est assez technique et dépasse le cadre ce ce


cours. Nous verrons toutefois dans quelle mesure ce critère assure le caractère entropique
des solutions faibles pour le problème de Riemann dans le cas des équations de Saint-Venant.

3.2 Problème de Riemann


Dans la perspective de la construction de schémas numériques, il est essentiel de caractériser
les solutions du problème de Riemann, qui correspond à un problème de Cauchy (3.2),
(3.2b) avec donnée initiale discontinue:
$

&Bt w ` Bx pF pwqq “ 0" , x P R, t ą 0 , (3.25a)
wG si x ă 0

% wpx, 0q “ (3.25b)
wD si x ą 0

La construction de la solution du problème de Riemann se base sur le résultat général (3.2.4),


que nous verrons à la fin de ce chapitre. Il est d’abord nécessaire d’introduire les notions
suivantes.

3.2.1 Champs caractéristiques


Dans la suite, on appellera k ième champ caractéristique la donnée du triplet
pλk pwq, rk pwq, lk pwqq, où lk pwq, rk pwq sont les vecteurs propres à gauche et à droite associés
à la valeur propre λk pwq. Nous verrons ultérieurement que la donnée des ces champs est
cruciale dans la construction de la solution du problème de Riemann, à la base des méthodes
numériques que nous allons considérer.
Définition 3.2.1: Nature des champs caractéristiques

On dit que le k ième champ caractéristique est vraiment non linéaire s’il vérifie:

∇λk [Link] pwq ‰ 0 , @w P Ω . (3.26)

On dit que le k ième champ caractéristique est linéairement dégénéré s’il vérifie:

∇λk [Link] pwq “ 0 , @w P Ω . (3.27)

63
Remarque 3.2.1. Revenons à l’exemple des équations ? de Saint-Venant. Les
? vecteurs pro-
pres à droite associés aux valeurs propres λ1 pwq “ u ´ gh et λ2 pwq “ u ` gh sont donnés
par: ˆ ˙ ˆ ˙
1? 1?
r1 pwq “ et r2 pwq “ . (3.28)
u ´ gh u ` gh
Il est alors aisé de vérifier que pour k “ 1, 2, ∇λk [Link] pwq ‰ 0 et ∇λk [Link] pwq ‰ 0 pour
tout w P Ω, de sorte que les deux champs caractéristiques sont vraiment non linéaires. A
titre d’exercice, on pourra vérifier que les équations de Saint-Venant avec traceur passif:
Bt h ` Bx phuq
ˆ “ 0, ˙
h2
2
Bt hu ` Bx hu ` g “ 0, (3.29)
2
Bt hϕ ` Bx phuϕq “ 0 ,
donnent lieu à un troisième champ caractéristique linéairement dégénéré. Nous verrons plus
tard que dans le cas des équations de Saint-Venant 2d, il existe aussi un troisième champ
caractéristique linéairement dégénéré.

3.2.2 Courbes intégrales - Ondes de détente


Une courbe intégrale associée à rk est une courbe c : ξ P R ÞÝÑ cpξq P Ω de classe C 8
vérifiant:
c1 pξq “ αpξqrk pcpξqq . (3.30)
En d’autres termes, la tangente à la courbe c en tout point est colinéaire à rk . On définit
les invariants de Riemann associés au champ caractéristique comme étant les fonctions
différentiables Ik : Rp Ñ R telles que ∇Ik [Link] pwq “ 0 pour tout w P Ω, que l’on appellera
k-invariants dans la suite. La terminologie “k-invariants” de Riemann tient du fait que ces
fonctions sont constantes le long des courbes intégrales de rk . En effet:
d
pIk pcpξqqq “ c1 pζq.∇Ik pcpξqq “ αpξqrk pcpζqq.∇Ik pcpξqq “ 0 . (3.31)

Il est alors important de préciser que dans le cas d’un champ linéairement dégénéré k (3.27),
la valeur propre λk est nécessairement un k´invariant de Riemann.

De manière générale, lorsqu’un choc n’est pas admissible, on peut tenter de relier deux
états wG et wD avec une solution régulière de la forme:
$ x

’ wG si ď σG ,

& x t x
wpt, xq “ νp q si σG ď ď σD (3.32)

’ t x t

% wD si σD ď ,
t
En revenant à la forme quasi-linéaire Bt w ` ApwqBx w “ 0, on a alors:
´ x ¯ ´x¯ ´ ´ x ¯¯ ˆ 1 ˙ ´ x ¯
1
´ 2 ν `A ν ν1 , (3.33)
t t t t t
64
x
d’où l’on tire, en posant ξ “ :
t
pApνpξqq ´ ξIq ν 1 pξq “ 0 .

Si l’on s’intéresse aux solutions non triviales, ceci implique que ξ est une valeur propre du
système associée au vecteur propre ν 1 pξq. En d’autres termes, il existe k tel que:

ν 1 pξq “ αpξqrk pνpξqq et ξ “ λk pνpξqq , (3.34)

et ν est alors une courbe intégrale associée à rk . On remarque alors qu’en dérivant la
deuxième équation: ν 1 pξq.∇λk pνpξqq “ 1, et donc:

αpξq∇λk pνpξ[Link] pνpξqq “ 1 ,

ce qui ne peut se produire si le champ est linéairement dégénéré. L’équation précédente se


réécrit:
1
αpξq “ ,
∇λk pνpξ[Link] pνpξqq
et en revenant à la première équation de (3.34):

rk pνpξqq
ν 1 pξq “ , (3.35)
∇λk pνpξ[Link] pνpξqq

système qu’on peut alors résoudre en s’aidant des conditions initiales. Notons aussi que la
deuxième relation de (3.34) impose σG “ λk pνpσG qq “ λk pwG q et de même σD “ λk pwD q.
Ainsi, nous venons de voir que deux états wG et wD peuvent être reliés par une solution
régulière si et seulement s’ils sont situés sur la courbe intégrale associée à un champ vraiment
non linéaire, avec la condition λk pwG q ă λk pwD q.

3.2.3 Discontinuités de contact


Supposons à présent le k ième champ caractéristique linéairement dégénéré, et cherchons à
connecter un état wD à un état wG donné. Nous avons vu que dans ce cas il n’est pas
possible de connecter ces deux états par une courbe régulière de la forme (3.32). Considérons
une une courbe intégrale c : ξ ÞÝÑ cpξq passant par wG associée à rk :
" 1
c pξq “ αpξqrk pcpξqq ,
(3.36)
c1 p0q “ wG .

Notons que λk étant un k´invariant de Riemann, on a λk pwG q “ λk pcpξqq. Le calcul donne:

d` ˘
F pcpξqq ´ F pwq ´ λk pwG q pcpξq ´ wq “ Apcpξqqc1pξq ´ λk pwG qc1 pξq
dξ ` ˘ (3.37)
“ αpξq Apcpξqq ´ λk pcpξqqI rk pcpξqq
“ 0,
65
ce qui signifie que la quantité F pcpξqq ´ F pwG q ´ λk pwG q pcpξq ´ wG q est constante. Con-
sidérant qu’elle s’annule en 0, on a donc:

F pcpξqq ´ F pwG q “ λk pwG q pcpξq ´ wq , (3.38)


ce qui signifie que les relations de Rankine-Hugoniot sont vérifiées le long de cette courbe
intégrale, avec une vitesse de propagation égale à λk pwG q. On en déduit que l’ensemble des
états pouvant être connectés à wG sont situés sur cette courbe intégrale. Par ailleurs, pour
tout état wD situé sur cette courbe, on aura σ “ λk pwD q “ λk pwG q, ce qui proscrit d’une
part une construction de la forme (3.32) (comme nous l’avons déjà évoqué). Il résulte de
cette analyse que la solution définie par
$ x
& wG si ă λk pwG q ,
wpx, tq “ t
% wD si x ą λk pwG q,
t
est une solution faible du problème, qui peut être vue comme une version dégénérée de (3.32)
avec σG “ σD “ λk pwG q. Il convient enfin de préciser que les k´invariants de Riemann sont
préservés à travers cette discontinuité, car ils sont constants sur les courbes intégrales.

3.2.4 Solution du problème de Riemann


Théorème 3.2.1: Théorème local d’existence et d’unicité (Lax)

On considère le problème de Riemann (3.25a,3.25b). Alors, si wG est suffisamment


proche de wD , et si le système n’admet que des champs vraiment non linéaires ou
linéairement dégénérées, il existe une unique solution faible entropique composée de p
ondes simples (choc, détente ou discontinuité de contact), séparant p ` 1 états con-
stants.

66
Chapter 4

Résolution pratique du problème de


Riemann

4.1 Equations de Saint-Venant 1d


On rappelle que le système de Saint-Venant sur fond plat s’écrit:

Bt w ` Bx pF pwqq “ 0 , (4.1)
¨ ˛
ˆ ˙ hu
h
avec w “ et F pwq “ ˝ 2 h2 ‚.
hu hu ` g
2
D’après le théorème (3.2.4), on sait que la solution du problème de Riemann consistera
en trois états constants séparés par deux ondes simples correspondant à des chocs ou des
détentes (il n’y a pas de discontinuité de contact), comme illustré dans la Figure 1.

ˆ ˙ ˆ ˙
hG ˆ ˙ hD
wG h ˚
wD
huG w˚ huD
b
hu˚ b
b

1-onde 2-onde
(Détente ou choc) (Détente ou choc)

Figure 1: Problème de Riemann - Connexion des états wG et wD .

4.1.1 Caractérisation des chocs


Pour fixer les idées, on s’intéresse aux chocs associés au premier champ caractéristique (on
renvoie à l’exemple (3.2.1) pour l’expression et la nature des champs). En d’autres termes,
on s’intéresse à caractériser les états w˚ pouvant être reliés à wG par un 1-choc admissible.

67
Par les relations de Rankine-Hugoniot, on a:
σph˚ ´ hG q “ h˚ u˚ ´ hG uG ,
1 1 (4.2)
σph˚ u˚ ´ hG uG q “ h˚ pu˚q2 ` gph˚q2 ´ hG u2G ´ gphG q2 .
2 2
˚ ˚
h u ´ hG uG
La première équation donne σ “ , et en injectant dans la deuxième:
h˚ ´ hG
ˆ ˆ ˙ ˙
˚ 2 ˚ 2 1 hG h˚ ˚
pu q ´ 2u uG ` uG ´ g ´ phG ´ h q , (4.3)
2 h˚ h
ce qui permet d’exprimer u˚ en fonction de h˚ :
d ˆ ˙
˚ g hG h˚
u “ uG ˘ ´ phG ´ h˚ q . (4.4)
2 h˚ hG
La théorie linéaire1 permet alors d’identifier la solution affectée du signe “-” comme celle
correspondant au premier champ caractéristique. Le calcul donne:
c
hG ` h˚ h˚
σG “ uG ´ sgnph˚ ´ hG q g , (4.5)
2 hG
Par ailleurs, les conditions d’admissibilité de Lax (3.24) imposent λ1 pwG q ě σ ě λ1 pw˚ q,
c’est à dire: a a
uG ´ ghG ě σ ě u˚ ´ gh˚ ,
d’où l’on déduit δ “ 1 et la condition nécessaire et suffisante h˚ ě hG . La relation (4.4) se
réécrit alors de la façon suivante:
d ˆ ˙
˚ ˚ ˚ g 1 1
u “ C1 ph q :“ uG ´ ph ´ hG q ` . (4.6)
2 h˚ hG
Par un raisonnement analogue, on montre qu’un 2-choc connectant un état w˚ à un état wD
se caractérise par :
d ˆ ˙
˚ ˚ ˚ g 1 1
u “ C2 ph q :“ uD ` ph ´ hD q ` , (4.7)
2 h˚ hD
à travers un choc se propageant à la vitesse
c
hD ` h˚ h˚
σD “ uD ` g , (4.8)
2 hD
sous la condition d’admissibilité h˚ ě hD .
Les relations (4.6) et (4.7) définissent des courbes C1 , C2 dans le plan ph, uq caractérisant
l’ensemble des états pouvant être respectivement connectés à wG et wD . Ces courbes sont
communément appelées courbes de Hugoniot et représentées sur la Figure 2. Lorsque ces
courbes s’intersectent, ceci signifie que les états wG et wD peuvent être connectés via un
état intermédiaire w˚ par le biais de deux chocs admissibles (Fig. 2 gauche), mais ce n’est
pas toujours le cas (Fig. 2 (droite)).
1
Il s’agit en fait d’étudier le comportement de la solution autour de petites perturbations hG “ h˚ ` α
et choisir un signe approprié pour observer la convergence vers w˚ dans une direction tangente au vecteur
propre r1 pw˚ q lorsque α tend vers 0.

68
u u
2 ´ choc pC2 q C2
1 ´ choc pC1 q w˚
b
b
wG
wD b

b
wG
˚
b
w
h h
wD b C1

Figure 2: Courbes de Hugoniot pour les 1-choc et 2-choc. Cas de deux chocs admissibles
(gauche) et d’un choc non admissible (droite).

4.1.2 Ondes de détente


Lorsqu’un choc n’est pas admissible, on cherche à connecter un état w˚ à wG par une 1-onde
de détente. Rappelons que pour tout vecteur d’état w “ ph, huqT on a ∇λ1 pwq.r1 pνpwqq “
3a
´ g{h de sorte que l’équation (3.35) s’écrit:
2
ˆ ˙
1 2a 1?
ν pξq “ ´ h{g . (4.9)
3 u ´ gh
2a
La première équation donne h1 pξq “ ´ hpξq{g, dont la solution générale est donnée par
3
1
hpξq “ pA ´ ξq2 , avec A une constante réelle2 . On utilise alors la condition hpξq “ hG
9g ?
pour ξ “ σG . En exploitant (3.34) et plus particulièrement σG “ λ1 pνpσG qq “ uG ´ ghG ,
on en déduit a
A “ uG ` 2 ghG . (4.10)
Ceci permet de déterminer h en fonction de ξ. Il suffit alors d’exploiter l’expression des
invariants de Riemann afin de déterminer la vitesse:
a a
uG ` 2 ghG “ u˚ ` 2 gh˚ ,

ce qui peut se réécrire:


´a a ¯
u˚ “ D1 ph˚ q :“ uG ` 2 ghG ´ gh˚ . (4.11)

Notons enfin que la condition λ1 pwG q ď λ1 pw˚ q doit être satisfaite pour pouvoir connecter
ces deux états, c’est à dire: a a
uG ´ ghG ă u˚ ´ gh˚ ,
2
Il s’agit d’un cas particulier d’une équation de Bernouilli : y 1 ` aptqy ` bptqy α “ 0, où a, b sont deux
fonctions continues et α P Rzt0, 1u, qui peut se résoudre en introduisant le changement de variable z “ y 1´α .

69
ou encore, avec (4.11):
a ´a a ¯ a
uG ´ ghG ă uG ` 2 ghG ´ gh˚ ´ gh˚ , (4.12)

ce qui revient à avoir h˚ ď hG . De manière similaire, les états w˚ pouvant être connectés à
wD à travers une 2-onde sont caractérisés par le biais des invariants de Riemann:
a a
uD ´ 2 ghD “ u˚ ´ 2 gh˚ , (4.13)

c’est à dire: ´a a ¯
˚ ˚
u “ D2 ph q :“ uD ´ 2 ˚
ghD ´ 2 gh , (4.14)
sous la condition h˚ ď hD . Une nouvelle fois, l’ensemble des états intermédiaires pouvant
être connectés à wG et wD par ces ondes de détentes peuvent être représentés dans le plan
ph, uq par les deux courbes de Hugoniot D1 , D2 définies par les relations (4.11) et (4.14).

u t
C2 σG σD
b
C1
wG w˚
wD
˚
b
w
wG
h
wD b
x

Figure 3: Cas choc - choc. Représentation dans le plan ph, uq (gauche) et dans le plan px, tq
(droite). Les vitesses des chocs σG , σD sont données par (4.5) et (4.8).

4.1.3 Bilan
Au final, il existe quatre cas possibles pour connecter wG et wD , selon si les 1-ondes et
2-ondes sont des détentes ou des chocs. Ces différents cas sont représentés dans les Figures
3 à 6 dans le plan ph, uq (gauche) ainsi que dans le plan px, tq (droite)

4.1.4 Un exemple : rupture de barrage


Il s’agit d’un cas d’étude très classique, correspondant à un cas particulier de problème de
Riemann pour lequel uG “ uD “ 0. Pour fixer les idées, considérons le cas hG ą hD ą
0. Nous savons que la solution du problème de Riemann va être constituée de trois états
constants reliés par deux ondes simples (choc ou ondes de détente). En d’autres termes, il
s’agit de déterminer un état intermédiaire w˚ permettant de relier wG à wD , soit par des

70
u wD b
D2 t
λ2 pw˚ q λ2 pwD q
σG
wG b C1

b w˚ wD
wG
h
x

Figure 4: Cas choc - détente. Représentation dans le plan ph, uq (gauche) et dans le plan
px, tq (droite). La vitesse du 1-choc σG est donnée par (4.5). La deuxième onde est une
2-détente se présentant sous la forme (3.32) avec les vitesses λ2 pw˚ q et λ2 pwD q.

chocs, soit par des détentes, portées par le premier et le deuxième champ caractéristique.
Concernant le premier champ caractéristique, l’étude générale nous a permis de dégager des
dépendances entre h˚ et u˚ , plus précisément:
$ `? ? ˚˘

’ D 1 ph ˚
q “ u G ` 2 ghG ´ gh si h˚ ă hG (onde de détente)
& d ˆ ˙
u˚ ph˚ q “ g 1 1
’ ˚ ˚
si h˚ ą hG (choc)
% C1 ph q “ uG ´ ph ´ hG q 2 h˚ ` hG

(4.15)
De manière similaire, les états pouvant être connectés à wD sont caractérisés par:
$ `? ? ˘

’ D2 ph˚ q “ uD ´ 2 ghD ´ gh˚ si h˚ ă hD (onde de détente)
& d ˆ ˙
u˚ ph˚ q “ g 1 1
’ ˚ ˚
si h˚ ą hD (choc)
% C2 ph q “ uD ` ph ´ hD q 2 h˚ ` hD

(4.16)
˚
Il s’agit donc de déterminer h de manière à rendre compatible les relations (4.15) et (4.16).
Après une étude de cas basique, ceci n’est possible que si hG ą h˚ ą hD , ce qui signifie que
la solution sera composée d’une 1-onde de détente et d’un 2-choc (cas de la Figure 5). On
est donc ramenés chercher la solution h˚ de l’équation non linéaire suivante :
d ˆ ˙
´a a ¯ g 1 1
˚
2 ˚
ghG ´ gh “ ph ´ hD q ` . (4.17)
2 h˚ hD

71
u t
C2 λ1 pwG q λ1 pw˚q σD
D1

w˚ wD
b

wG
h
wD b b
x
wG

Figure 5: Cas détente - choc. Représentation dans le plan ph, uq (gauche) et dans le plan
px, tq (droite). La première onde est une 1-détente se présente sous la forme (3.32) avec les
vitesses λ2 pwG q et λ2 pw˚q. La vitesse du 2-choc σD est donnée par (4.8).

4.2 Le cas 2d
4.2.1 Généralités
En regroupant les variables d’écoulement dans la quantité w “ ph, hu, hvqT , supposée ap-
partenir à l’ensemble (convexe) des états admissibles:
(
Ω :“ ph, hu, hvqT , h ą 0, hu, hv P R Ă R3 ,

les équations (6.21) se reformulent sur fond plat sous la forme compacte:

Bt w ` div pGpwqq “ 0 , (4.18)

où la fonctions de flux G : Ω ÞÝÑ R3 ˆ R3 est définie par:


¨ ˛
hu hv
˚ 2 1 2 ‹
˚hu ` gh huv ‹
où Gpwq “ ˚ 2 ‹. (4.19)
˝ 1 ‚
huv hv 2 ` gh2
2
Comme dans le cas 1d ces équations peuvent se réécrire sous forme quasi-linéaire:

Bt w ` Ax pwqBx w ` Ay pwqBy w “ 0 , (4.20)

où les matrices Ax pwq, Ay pwq sont données par:


¨ ˛ ¨ ˛
0 1 ´uv 0 0 1
Ax pwq “ ˝gh ´ u2 2u v ‚ , Ay pwq “ ˝ ´uv v u ‚. (4.21)
2
0 0 u gh ´ u 0 2v

72
u wD b
D2 t
λ1 pwG q λ1 pw˚q λ2 pw˚ q λ2 pwD q


D1
b

wD
wG
h
b
x
wG

Figure 6: Cas choc - choc. Représentation dans le plan ph, uq (gauche) et dans le plan px, tq
(droite). Les vitesses des chocs σG , σD sont données par (4.5) et (4.8).

Le système de Saint-Venant 2d (4.18) est strictement hyperbolique. Plus précisément, on


peut montrer que pour tout vecteur ~n “ pnx , ny q, la matrice A~n pwq “ nx Ax pwq ` ny Ay pwq
est diagonalisable et admet trois valeurs propres distinctes si w P Ω. Cette analyse peut se
poursuivre en dimension quelconque et permet d’étendre la notion d’hyperbolicité au cas des
systèmes multi-dimentionnels.

De façon similaire, notons aussi que les notions de solutions fortes et solutions faibles, ainsi
que les définitions de solutions entropiques s’étendent au cas de 2d de manière naturelle.

Dans la perspective de l’étude du problème de Riemann, nous nous placerons dans un con-
texte "pseudo-1d ", en supposant l’écoulement constant selon l’axe y, et noterons
ˆ ˙T
2 1 2
F pwq “ hu, hu ` gh , huv (4.22)
2
la composante horizontale du flux G, de sorte que le problème de Chauchy correspondant
s’écrit:
Bt w ` Bx pF pwqq “ 0 , x P R, t ą 0 . (4.23)
Remarque 4.2.1. Il est important de préciser ici que cette cette stratégie d’analyse trouve
ses motivations d’un point de vue numérique, car même en deux dimensions, la construction
de schémas s’effectue en se ramenant à des problèmes de Riemann 1d dans une direction
spatiale privilégiée (portée par la direction normale à l’interface séparant les états "gauche"
wG et "droit" wD ).

4.2.2 Champs caractéristiques


La matrice jacobienne associée au problème (4.22), (4.23) est donnée par:
¨ ˛
0 1 0
˝
Ax pwq “ gh ´ u2 2u 0 ‚
´uv v u
73
Une étude similaire au cas 1d permet de dégager la nature des champs caractéristiques, ainsi
que les invariants de Riemann associés. Le calcul donne:
¨ ˛
? 1
?
λ1 pwq “ u ´ gh , r1 pwq “ ˝u ´ gh‚
¨ ˛ v
0
λ2 pwq “ u , r2 pwq “ 0‚ ˝ (4.24)
¨1 ˛
? 1
?
λ3 pwq “ u ` gh , r3 pwq “ ˝u ` gh‚
v

Le premier et le troisième champ caractéristique


? sont vraiment?non linéaires. Les invariants
de Riemann associés sont respectivement u ` 2 gh, v et u ´ 2 gh, v. Il est facile de vérifier
que le deuxième champ caractéristique est linéairement dégénéré, et admet u et h pour
invariants de Riemann.

4.2.3 Problème de Riemann


Considérant la remarque (4.2.1), le problème de Riemann associé est donné par (3.25a),
(3.25b), avec la fonction de flux donnée par (4.22). Sur la base du Théorème (3.2.4), la
solution du problème de Riemann sera constituée de quatre états constants séparés par 3
ondes simples. Considérant la nature des champs caractéristiques, la solution sera constituée
d’un choc ou d’une détente, d’une discontinuité de contact, et enfin d’un choc ou d’une
détente, comme illustré dans la Figure 7.

¨ ˛ ¨ ˛ ¨ ˛
hG h˚ ¨ ˛ hD
˚
wG ˝huG ‚ ˚ ˝
wG hu˚ ‚ h
wD ˝huD ‚
hvG h˚ vG
˚
wD ˝ hu˚ ‚
b
h˚ vD hvD
b b

1-onde b

(Détente ou choc) 2-onde 3-onde


(Discontinuité de contact) (Détente ou choc)

Figure 7: Problème de Riemann - Connection des états wG et wD .

4.2.4 Caractérisation des chocs et des détentes


L’étude du premier et du troisième champ caractéristique se mène de manière analogue au
cas 1d. Concernant les chocs, les conditions (4.2) sont strictement identiques, et on vérifie
par le calcul que la troisième condition:

σph˚ v ˚ ´ hG vG q “ h˚ u˚ v ˚ ´ hG uG vG
74
implique v ˚ “ vG . Pour les détentes, les calculs sont inchangés; il suffit d’utiliser le fait que v
est un invariant de Riemann à travers le premier champ caractéristique pour avoir v ˚ “ vG .
De manière similaire, v ˚ “ vD à travers un 3-choc ou une 3-onde de détente.

4.2.5 Discontinuité de contact


˚ ˚
Examinons le cas de la deuxième onde. On notera respectivement wG et wD les états
connectés à wG à travers la première onde et à wD à travers la troisième. Nous venons
de voir que l’on ne peut pas relier ces deux états intermédiaires par une solution continue à
travers cette onde, ni par un choc admissible. Revenant à l’étude effectuée en Section (3.2.3),
la solution sera sous la forme:
$ x
˚ ˚
& wG si ă λk pwG q,
wpx, tq “ t (4.25)
% w˚ si x ą λk pw˚ q ,
D G
t
˚ ˚
avec wG et wD tous deux situés sur une courbe intégrale associée au deuxième champ car-
actéristique. Notons que h et u sont des quantités constantes à travers cette onde car ce
˚
sont tous deux des 2-invariants de Riemann. Par ailleurs, nous avons vu que vG “ vG et
˚ ˚ ˚ ˚ ˚ ˚ ˚
vD “ vD . On note alors h “ hG “ hD et u “ uG “ uD , de sorte que la solution consiste
simplement en un saut de vitesse normale, comme illustré sur le schéma présenté en Figure. 7.

Ceci permet de caractériser les solutions du problème de Riemann d’une manière analogue
au cas 1d, à la différence près que le saut de vitesse normale est pris en compte par la
discontinuité de contact.

75
Chapter 5

Résolution numérique des équations de


Saint-Venant

On s’intéresse dans cette partie au système 1d sans gradient de fond:


$
& Bt h ` Bx huˆ“ 0 , ˙
2 1 2
% Bt hu ` Bx hu ` gh “ 0 ,
2
que l’on peut écrire sous forme compacte Bt w ` Bx f pwq “ 0.

5.1 Exemples de schémas numériques


On considère des schémas Volumes Finis conservatifs, sous la forme:
∆t ` n ˘
wjn`1 “ wjn ´ n
gj`1{2 ´ gj´1{2 , (5.1)
∆x
où gj`1{2
n
“ gpwjn , wj`1
n
q. A noter que wjn est une quantité vectorielle et g : R2 ˆ R2 Ñ R2 .
La notion de consistance au sens faible s’étend directement au cas système: gpw, wq “
f pwq , @w P R2 .

5.1.1 Schéma de Godunov


Le schéma de Godunov s’étend de manière naturelle au cas système, ce qui exige la résolution
préalable du problème de Riemann associé. Le schéma de Godunov s’écrit donc sous la forme
(5.1), avec gj`1{2
n
“ gpwjn, wj`1
n
q “ f pw ˚q, où w ˚ “ wR p0, tq est la solution exacte du problème
de Riemann entre wjn et wj`1 n
évaluée en x “ 0.
Remarque 5.1.1. Le pas de temps est alors contraint par les valeurs propres λk pwjn q,
n
λk pwj`1 q, où k “ 1, 2, qui représentent les vitesses de propagation des ondes du système.
Plus précisément:
∆t 1
max |λk pwjn q| ď , k “ 1, 2 ,
jPZ ∆x 2
? ?
où λ1 pwq “ u ´ gh et λ2 pwq “ u ` gh.
76
Comme dans le cas scalaire (et c’est d’autant plus valable dans le cas systèmes), la réso-
lution exacte du problème de Riemann peut s’avérer complexe et coûteuse numériquement.
On considère alors des solveurs approchés, dont le principe de construction est calqué sur le
cas scalaire.

5.1.2 Schéma de Rusanov


Les notions de consistance avec la forme intégrale et l’inégalité d’entropie vues dans le cas
scalaire se généralisent de manière directe. La consistance avec la forme intégrale s’écrit:
ż ∆x
1 2 1 ∆t
ŵR px{∆t, wL , wR qdx “ pwL ` wR q ´ pf pwR q ´ f pwL qq , (5.2)
∆x ´ ∆x
2
2 ∆x

à ceci près qu’il s’agit là encore d’égalités vectorielles. Sur cette base, il est possible de
reproduire l’analyse effectuée dans le cas scalaire pour obtenir le schéma de Rusanov:

n 1` ˘ λnj`1{2 ` n ˘
gj`1{2 “ gpwjn , wj`1
n
q “ n n
f pwj q ` f pwj`1q ´ wj`1 ´ wjn , (5.3)
2 2
avec λnj`1{2 “ max |λk pwq| , k “ 1, 2 .
wPtwjn ,wj`1
n
u

5.1.3 Schéma HLL (Harten - Lax - Van Leer)


Si on a une estimation de la plus grande et de la plus petite valeur propre impliquées dans
le problème de Riemann (notées respectivement λ` et λ´ ), on peut envisager de construire
un solveur approché moins diffusif, sous la forme suivante:
$
& wL si x{t ă λ´
ŵR px{t, wL , wR q “ w ˚ si λ´ ă x{t ă λ`
%
wR si λ` ă x{t

En revenant à (5.2), on peut déterminer w ˚ de manière à obtenir la consistance avec la forme


intégrale: t
´
λ`
λ

∆t w˚
wL wR

´∆x{2 λ´ ∆t λ` ∆t ∆x{2 x

Figure 1: Cône de dépendance associé au solveur approché.

77
ż „ ¯ 
1 ´ ∆x ¯ ´ ¯ ´ ∆x
∆x
1 2
´ ` ´ ˚ `
ŵR px{∆t, wL , wR qdx “ ` λ ∆t wL ` ∆t λ ´ λ w ` ´ λ ∆t wR
∆x ´ ∆x ∆x 2 2
∆t ´ ` ¯
2
wL ` wR ´ ˚ ´ `
“ ` pλ ´ λ q w ` λ wL ´ λ wR .
2 ∆x
´ ¯
On en tire pλ` ´ λ´ q w ˚ ` λ´ wL ´ λ` wR “ ´ f pwR q ´ f pwL q , c’est à dire:
´ ¯
λ` wR ´ λ´ wL ´ f pwR q ´ f pwL q
w˚ “ . (5.4)
λ` ´ λ´
Par analogie avec le cas scalaire, le schéma s’écrit sous la forme (5.1), avec
ż
n n ∆x n 1 xj`1{2
gj`1{2 “ f pwj q ` w ´ ŵR ppx ´ xj`1{2 q{∆t, wjn , wj`1
n
qdx
2∆t j ∆t xj
ż
n ∆x n 1 xj`1
“ f pwj`1q ´ w ` ŵR ppx ´ xj`1{2 q{∆t, wjn , wj`1
n
qdx .
2∆t j`1 ∆t xj`1{2

Si λ´ ą 0, alors ŵR ppx ´ xj`1{2 q{∆t, wjn , wj`1


n
@x P rxj , xj`1{2 s, et donc en utilisant
q “ wjn
∆x n ∆x n
la première ligne dans l’égalité précédente: gj`1{2
n
“ f pwjn q ` wj ´ wj “ f pwjn q. On
2∆t 2∆t
montre de même en utilisant la deuxième définition que si λ` ă 0, alors gj`1{2 n n
“ f pwj`1 q.
´ `
Enfin, si λ ă 0 ă λ , on peut utiliser la première expression par exemple, pour obtenir:
n
gj`1{2 “ f pwjn q ´ λ´ wjn ` λ´ w ˚ ,

ce qui donne, avec (5.4):

n
λ` f pwjn q ´ λ´ f pwj`1
n
q ´ λ´ λ` pwj`1
n
´ wjn q
gj`1{2 “ .
λ` ´ λ´
Au final:
$

’ f pw n q si λ´ ą 0
& ` j n
n λ f pwj q ´ λ´ f pwj`1
n
q ´ λ´ λ` pwj`1
n
´ wjn q
gj`1{2 “ si λ´ ă 0 ă λ`

’ λ ` ´ λ´
% f pw n q si λ` ă 0
j`1

A noter que λ˘ est une estimation locale de la plus grande et de la plus petite valeur propre
à l’interface xj`1{2 . Dans le cas des équations de Saint-Venant : λ˘ “ λ˘
j`1{2 , avec:
´ a ¯ ´ a ¯
λ´
j`1{2 “ min u´ gh , λ`
j`1{2 “ max u ` gh .
wPtwjn ,wj`1
n
u wPtwjn ,wj`1
n
u

78
5.2 Stabilité
Les notions de stabilité L8 , L1 et de schéma TVD s’étendent naturellement au cas système,
de même que la notion de stabilité entropique. En particulier, il existe un théorème de
type Lax-Wendroff assurant la convergence de schémas entropiques vers la solution faible
entropique du problème, sous certaines hypothèses de régularité sur la solution discrète (et
sous réserve de la convergence du schéma bien sûr). On s’intéresse particulièrement ici
à la stabilité L2 , ainsi qu’à d’autres critères de stabilité non-linéaire, plus spécifiques aux
équations de Saint-Venant.

5.2.1 Positivité
Il s’agit ici d’étudier la capacité des schémas à maintenir la solution
( discrète dans l’ensemble
convexe des états admissibles Ω “ ph, huq , h ą 0, hu P R Ă R . On dit qu’un schéma
T 2

préserve la positivité (de la hauteur d’eau) si:

hnj ą 0 , @j P Z ñ hjn`1 ą 0 , @j P Z .

Remarque 5.2.1. Propriété indispensable en pratique, car elle assure la stabilité des sché-
mas au voisinage des fronts secs. Configurations couramment rencontrées par exemple en
hydraulique fluviale ou en océanographie côtière pour étudier les phénomènes de submersion,
crues, “run-up” et “run-down” associés au va-et-vient des vagues.

Exemples : Les schémas de Godunov, Rusanov, Lax-Friedrichs et HLL préservent la


positivité sous CFL. Pour le schéma de Godunov, c’est évident en revenant à la définition
moyennée du schéma: la solution exacte wR du problème de Riemann associé à deux états
wL , wR tels que hL ą 0 et hR ą 0 vérifie hR ą 0. Ce raisonnement s’applique aussi aux
solveurs approchés tels que Rusanov, Lax-Friedrichs et HLL, sous réserve de la positivité
de l’état intermédiaire w ˚ . On peut toutefois mettre ceci en évidence de manière directe en
revenant à l’expression classique des schémas, par exemple, pour le schéma de Rusanov:
∆t ` n ˘
hjn`1 “ hnj ´ n
gj`1{2 ´ gj´1{2 ,
∆x
1` ˘ λj`1{2 ` n ˘
avec gj`1{2
n
“ n
f pwj`1 q ` f pwjn q ´ wj`1 ´ wjn . Le calcul donne:
2 2
ˆ ˙
∆t ` ˘ ∆t ` ˘ ∆t ` ˘
hjn`1 “ hnj 1 ´ λj`1{2 ` λj´1{2 `hnj`1 λj`1{2 ´ unj`1 `hnj´1 λj´1{2 ` unj´1 .
2∆x 2∆x 2∆x
?
Rappelons alors que λj`1{2 “ max |u˘ gh| ą maxpuj`1, uj q, ce qui assure la positivité
wPtwjn ,wj`1
n
u
∆t ` ˘
des deux derniers coefficients. Par ailleurs, 1 ´ λj`1{2 ` λj´1{2 ě 0 sous la contrainte
2∆x
∆t n a n
max |uj ˘ ghj | ď 1.
jPZ ∆x

79
5.2.2 États stationnaires
On s’intéresse à présent aux solutions stationnaires du système avec topographie:
$
& Bt h ` Bx huˆ“ 0 , ˙
2 1 2
% Bt hu ` Bx hu ` gh “ ´ghBx z .
2
Ces solutions sont caractérisées par les équations suivantes, parfois appelées équations de
Bernoulli dans la littérature:
$
& Bx hu
ˆ “ 0, ˙
2 1 2
% Bx hu ` gh “ ´ghBxz .
2
ˆ 2˙
u
Remarquons alors que l’on a Bx ` gBx η “ 0, où η “ h ` z. Ces états d’équilibre
2
traduisent un équilibre entre la force de pression (énergie potentielle) et l’énergie cinétique.
Un autre critère de stabilité numérique consiste à demander aux schémas de préserver ces
équilibres au niveau discret. Ces équilibres sont difficile à résoudre, même au niveau continu.
C’est pourquoi on se place généralement dans le cadre simplifié u “ 0, qui constitue déjà un
critère de stabilité essentiel. Cet équilibre est caractérisé par:

u“0 , η “ cte .

L’interprétation est la suivante : on considère un niveau d’eau total constant et une vitesse
nulle sur tout le domaine, avec des conditions aux limites de type paroi. En l’absence de
forçage extérieur (vent, pluie, débit entrant ou sortant, ...), on attend à ce que le schéma
maintienne cet état d’équilibre, parfois appelé lac au repos dans la littérature. Cette propriété

η “h`z
hpx, tq u“0

bpxq
x
Figure 2: État d’équilibre type lac au repos.

se formule donc ainsi:


´ ˆ n˙ ¯ ´ ¯
n hj
wj “ , avec hnj ` zj “ cte @j P Z ñ wjn`1 “ wjn @j P Z . (5.5)
0
80
Approche naïve
Pour mettre le problème en évidence, on considère le schéma de Lax-Friedrichs avec une
discrétisation centrée du terme source:
˜ ¸
∆t ` ˘ ´ 0 ¯
wjn`1 “ wjn ´ gn n
´ gj´1{2 ´ ∆tSjn , avec Sjn “ zj`1 ´ zj´1 .
∆x j`1{2 ghnj
2∆x
On a donc :
$ ˆ ˙
’ ∆t hunj`1 ´ hunj´1 1` n ˘

’ hn`1
“ hn
´ ` h ´ 2hn
` h n


j j
∆x ˜ 2 2 j`1 j j´1

& ¸
hu n hu n
∆t f pw j`1 q ´ f pw j´1 q
’ hujn`1 “ hunj ´

’ 2∆x 2



% 1` ˘ ∆t n ´ zj`1 ´ zj´1 ¯
` hunj`1 ´ 2hunj ` hunj´1 ´ g h
2 ∆x j 2
Considérons alors un état d’équilibre unj “ 0 , hnj ` zj “ cte, pour tout j P Z. Le schéma
donne:
$
’ 1` n ˘
& hjn`1 “ hnj ` hj`1 ´ 2hnj ` hnj´1
2 „ ˆ ˙ 
’ n`1 n g ∆t ` n n
˘ hnj`1 ` hnj´1 n
% huj “ huj ´ hj`1 ´ hj´1 ` hj pzj`1 ´ zj´1 q
2 ∆x 2
On remarque que la condition hnj ` zj “ cte @j P Z permet de réécrire la deuxième équation
sous la forme:
g ∆t ` n ˘
hujn`1 “ hunj ´ hj`1 ´ 2hnj ` hnj´1 .
4 ∆x
On voit alors clairement que l’équilibre est perturbé par la quantité hnj`1 ´2hnj `hnj´1 , qui n’est
pas nulle en général. Cette perturbation provient du terme de diffusion sur h ainsi que d’une
discrétisation inappropriée des flux et du terme source dans la quantité de mouvement: la
solution numérique est alors non physique. Depuis plusieurs décennies, de nombreux travaux
ont été menés afin d’apporter des solutions au problème.

Reconstruction hydrostatique
Il s’agit d’une méthode générale s’appliquant à tout schéma associé à des flux consistants.
L’idée est d’injecter des états reconstruits dans les flux numériques (afin d’en exploiter la
consistance), tout en proposant une discrétisation appropriée du terme source. Étant donné
un flux numérique consistant gpu, vq, ce schéma s’écrit:
∆t ` n ˘
wjn`1 “ wjn ´ n
g̃j`1{2 ´ g̃j´1{2 ´ ∆tS̃jn ,
∆x
´ ` ˘
avec g̃j`1{2
n
“ gpwj`1{2 , wj`1{2 q, où wj`1{2 désignent des états intermédiaires définis à l’interface
j ` 1{2. On définit ces états intermédiaires comme suit:
h` n
j`1{2 “ hj`1 ` zj`1 ´ zj`1{2 , hu` ` n
j`1{2 “ hj`1{2 uj`1 ,
h´ n
j`1{2 “ hj ` zj ´ zj`1{2 , hu´ ´ n
j`1{2 “ hj`1{2 uj ,

81
et zj`1{2 » zpxj`1{2 q correspond à une approximation de z à l’interface. Considérons alors
un état d’équilibre du type lac au repos: unj “ 0 , hnj ` zj “ η, pour tout j P Z. On a alors:
ˆ ` ˙ ˆ n ˙ ˆ n ˙ ˆ ´ ˙
` hj`1{2 hj`1 ` zj`1 ´ zj`1{2 hj ` zj ´ zj`1{2 hj`1{2 ´
wj`1{2 “ “ “ “ “ wj`1{2
0 0 0 0
` ´
On note alors wj`1{2 “ wj`1{2 “ wj`1{2 . Les flux numériques étant supposés consistants, on
ˆ ˙
´ ` 0
a g̃j`1{2 “ gpwj`1{2, wj`1{2q “ f pwj`1{2q “
n
, et par suite:
gphj`1{2q2 {2
$ n`1
& hj “ hnj ˆ ˙
n`1 n ∆t phj`1{2q2 ´ phj´1{2 q2
% huj “ huj ´ g ´ ∆tS̃jn .
∆x 2
Il reste à définir le terme source S̃jn de manière à compenser les termes provenant du flux
numérique. On introduit alors la définition générale:
˜ ´ 2 ` 2
¸
g phj`1{2 q ´ phj´1{2 q
S̃jn “ ´ ,
2 ∆x

g ∆t ` ˘
de sorte que ∆tS̃jn “ ´ phj`1{2q2 ´ phj´1{2 q2 dans un contexte de lac au repos, assur-
2 ∆x
ant ainsi un équilibre exact.
zj`1 ` zj
Remarque 5.2.2. En posant zj`1{2 “ , le calcul donne:
2
˜ ´ ¸˜ ´ ¸
hj`1{2 ` h` j´1{2 hj`1{2 ´ h`
j´1{2
S̃jn “ ´g
2 ∆x
˜ ´ ¸
`
hj`1{2 ` hj´1{2 ´ zj`1 ´ zj´1 ¯
“g
2 2∆x
n
» g hpxj , t qBx zpxj q ,
de sorte que la discrétisation du terme source est consistante avec le terme source continu.
En réalité, le choix centré ne permet pas de récupérer un schéma stable et on prend zj`1{2 “
maxpzj , zj`1 q afin d’assurer des inégalités d’entropie semi-discrètes.
Remarque 5.2.3. Depuis peu, il existe des schémas capables de préserver des états d’équilibre
plus complexes (u ‰ 0), voire même tous les états stationnaires, mais leur analyse et leur
mise en oeuvre est bien plus fastidieuse.

Équations pre-balanced
L’idée est de remplacer directement h par η “ h ` z dans les équations. Ce changement
de variables est une simple translation et n’altère pas la structure du système (il reste hy-
perbolique, mêmes valeurs propres, même nature des champs caractéristiques et solution
problème de Riemann, ...). Le système se réécrit:
Bt v ` Bx gpvq “ T pvq ,
82
¨ ˛
ˆ ˙ hu ˆ ˙
η ˝ 2 2 ‚ 0
avec v “ , gpvq “ pη ´ zq phuq , et T pvq “ .
hu g ` ´gpη ´ zqBx z
´g 2 ¯ 2pη ´ zq ´g ¯
En utilisant la relation Bx 2
pη ´ zq `gpη ´zqBx z “ Bx pη 2 {2 ´ ηzq `gηBx z, on aboutit
2 2
à une version simplifiée du modèle:

Bt v ` Bx fˆpvq “ Ŝpvq ,
¨ ˛
hu ˆ ˙
0
avec f pvq “ ˝
ˆ phuq ‚, et Ŝpvq “
2
. On peut aisément établir
g pη 2 {2 ´ ηzq ` ´gηBxz
2pη ´ zq
que ce système est équivalent au sens faible au système d’origine si z est de classe C 1 .
Considérons à présent un schéma consistant pour ce modèle:
∆t ` n ˘
vjn`1 “ vjn ´ n
ĝj`1{2 ´ ĝj´1{2 ´ ∆tŜjn ,
∆x

avec ĝj`1{2
n n
“ ĝpvj`1 , vjn q, et ĝ consistant (par rapport à fˆ, i.e. , ĝpv, vq “ fˆpvqq. On injecte
ˆ n˙
ηj
alors une solution à l’équilibre vjn “ , avec ηjn “ η pour tout j P Z. Tous les états au
0
temps n sont donc identiques, et la consistance des flux donne:
ˆ ˙
n n n ˆ n ` 0 ˘
ĝj`1{2 “ ĝpvj`1, vj q “ f pvj q “ ,
g η 2 {2 ´ ηzj`1{2

zj`1{2 désignant une approximation de z à l’interface j ` 1{2. On a donc:


#
ηjn`1 “ ηjn
∆t ` ˘
hujn`1 “ hunj ` gη zj`1{2 ´ zj´1{2 ´ ∆tŜjn .
∆x
ˆ ˙
zj`1{2 ´ zj´1{2
Il suffit alors de définir Ŝj “ gηj
n n
pour garantir un équilibre avec le terme
∆x
source, et assurer ainsi hujn`1 “ hunj . Cette méthode permet d’éviter d’avoir recours à une
étape de reconstruction au niveau des flux.

5.2.3 Stabilité L2
Même dans le cas de systèmes non-linéaires, il est possible d’analyser la stabilité L2 des
schémas par une approche similaire à celle vue dans le cas scalaire linéaire. Ceci fournit une
condition nécessaire de stabilité qui constitue en général un bon indicateur de la contrainte
à imposer sur le pas de temps. A titre d’exemple, on étudie ici le schéma de Lax-Friedrichs:
∆t ´ n ¯
wjn`1 “ wjn ´ n
gj`1{2 ´ gj´1{2 ,
∆x

83
1` ˘ ∆x ` n ˘
avec gj`1{2
n
“ f pwjn q ` f pwj`1
n
q ´ wj`1 ´ wjn . Le schéma se réécrit :
2 2∆t
$ ˆ n n ˙
’ ∆t qj`1 ´ qj´1 1` n ˘

’ h n`1
“ h n
´ ` h ´ 2h n
` hn


j j
∆x ˜ 2 2 j`1 j j´1

& n 2 n 2
` n ˘2 ` n ˘2 ¸
n`1 n ∆t pqj`1 q pqj´1 q g hj`1 g hj´1
’ q j “ q j ´ ´ ` ´

’ 2∆x hnj`1 hnj´1 2 2

’ ` ˘

% 1 n
` qj`1 ´ 2qjn ` qj´1
n
2
On linéarise alors ce schéma autour d’un état d’équilibre ph̄, q̄q; on pose:

hnj “ h̄ ` h̃nj , qjn “ q̄ ` q̃jn ,

le symbole “˜” désignant une petite perturbation autour de l’équilibre. Notons alors que le
schéma sur h est déjà sous forme linéaire: en utilisant l’expression précédente, on obtient
directement:
ˆ n n ˙
n`1 n ∆t q̃j`1 ´ q̃j´1 1´ n ¯
h̃j “ h̃j ´ ` h̃j`1 ´ 2h̃nj ` h̃nj´1 . (5.6)
∆x 2 2

Concernant le schéma sur q, on a des contributions non linéaires provenant de la partie


centrée des flux numériques:
n
` n
˘2 n
pqj´1 q2 q̄ ` q̃j´1 q̄ 2 ` 2q̄pq̃j´1 q
n
“ “ ` Opw̃ 2q
hj´1 n
h̄ ` h̃j´1 n
h̄˜` h̃j´1 ¸
ˆ 2 n ˙
q̄ ` 2q̄pq̃j´1 q h̄
“ n
` Opw̃ 2 q
h̄ h̄ ` h̃
˜ ¸
j´1
` ˘ n
n
h̃j´1
“ q̄ū ` 2ūpq̃j´1 q 1´ ` Opw̃ 2 q

n
“ q̄ū ` 2ūpq̃j´1 q ´ h̃nj´1 ū2 ` Opw̃ 2 q ,

où l’on a noté ū “ q̄{h̄. D’autre part:

phnj´1q2 ph̃nj´1 ` h̄q2 h̄2


g “g “ g ` g h̄h̃nj´1 ` Opw̃ 2q .
2 2 2
On en déduit le schéma linéarisé:

n`1 n ∆t ´ ` n n
˘ ` 2
˘´ n n
¯¯
q̃j “ q̃j ´ 2ū q̃j`1 ´ q̃j´1 ` g h̄ ´ ū h̃j`1 ´ h̃j´1
2∆x (5.7)
`
1 n ˘
` q̃j`1 ´ 2q̃jn ` q̃j´1
n
2
Rappelons que l’analyse de stabilité L2 garantit un contrôle uniforme sur la norme L2 de la
solution. Cette notion peut ainsi être utilisée pour garantir des estimations discrètes sur de

84
quantités conservées par le modèle, et en particulier l’énergie. En introduisant le changement
de variables:
a
h̃nj Ñ h̃nj “ h̃nj g{2
a
n n n
q̃j Ñ q̃j “ q̃j { 2h̄ ,

la norme L2 de la solution s’écrit:


¨ ´ ¯2 ˛
` ˘
ÿ ´` ˘2 ` ˘2 ¯ ÿ ˚ h̃j
n
1 q̃jn ‹
2

}un }L2 “ hjn`1 ` qjn`1 “ ˝g ` ‚,


jPZ jPZ
2 h̄ 2

et peut être associée à une énergie discrète. L’équation de la masse (5.6) et l’équation de
quantité de mouvement (5.7) se réécrivent:
ˆ ˙
n`1 n ∆t qnj`1 ´ qnj´1 1´ n n n
¯
h̃j “ h̃j ´ c ` h̃ ´ 2h̃j ` h̃j´1 , (5.8)
∆x 2 2 j`1
∆t ´ ` n ˘ ´ ¯¯
q̃jn`1 “ q̃nj ´ n 2 n n
2ū q̃j`1 ´ q̃j´1 ` pc ´ ū {cq h̃j`1 ´ h̃j´1
2∆x (5.9)
1` n ˘
` q̃j`1 ´ 2q̃nj ` q̃j´1
n
,
2
a
où c “
ˆ n ˙ g h̄.ˆ On ˙ considère alors la fonction constante par morceaux définie par W n pxq “
h pxq h̃nj
n “ si x Psxj´1{2 , xj`1{2 r, de sorte que (5.8) se réécrit :
q pxq q̃nj
ˆ ˙
n`1 n ∆t q n px ` ∆xq ´ q n px ´ ∆xq 1
h pxq “ h pxq´c ` phn px ` ∆xq ´ 2hn pxq ` hn px ´ ∆xqq ,
∆x 2 2
ou encore:
ˆ ˙
n`1 ∆t
n τ∆x q n ´ τ´∆x q n 1
h “h ´c ` pτ∆x hn ´ 2hn ` τ´∆x hn q ,
∆x 2 2

où l’on a introduit l’opérateur τδ défini par τδ f : x ÞÝÑ f px ` δq. De même, l’équation sur la
quantité de mouvement (5.9) se réécrit:
∆t ` ` ˘ ˘
q n`1 “ q n ´ 2ū pτ∆x q n ´ τ´∆x q n q ` c ´ ū2 {c pτ∆x hn ´ τ´∆x hn q
2∆x
1
` pτ∆x q n ´ 2q n ` τ´∆x q n q .
2

On considère alors la transformée de Fourier de W n (on rappelle que τz ∆x f pξq “ e


i∆xξ p
f pξq).
On obtient:
¨ ˛
1 iω ´iω 1 iω ´iω
1 ` pe ´ 2 ` e q ´ cα pe ´ e q
x n`1 pξq “ ˚
W ˝1 2 2 ‹ xn
‚W pξq ,
1
α pū2 {c ´ cq peiω ´ e´iω q 1 ´ αū peiω ´ e´iω q ` peiω ´ 2 ` e´iω q
2 2
85
∆t
où l’on a posé ω “ ∆xξ et α “ . On tombe ainsi sur une expression de la forme
∆x
x n`1 “ AW
W x n , d’où l’on tire W
x n “ An W
x 0 , et par suite }W x 0}L2 , où ~.~
x n }L2 “ ~An ~}W
désigne la norme subordonnée L2 . On a donc stabilité L2 si et seulement s’il existe une
constante K ą 0 telle que ~An ~ ă K pour tout ξ P R et tout n P N. En utilisant le
résultat classique ρpAqn ď ~An ~ ď ~A~n , on obtient ρpAq ď 1 pour condition suffisante, et
~A~ ď 1 pour condition nécessaire. En pratique, on étudie le spectre de A afin d’obtenir
des conditions nécessaires. Dans notre cas, une analyse basique donne SppAq “ tX ˘ u “
! ´ a ¯) 1 1
c ´ iαs ū ˘ g h̄ , avec c “ peiω ` e´iω q “ cospωq et s “ peiω ´ e´iω q “ sinpωq. On
2 2i
∆t ˘ a
˘
vérifie alors que la condition |X | ď 1 se réécrit |λ | ď 1, où λ˘ “ ū ˘ g h̄. On retrouve
∆x
ainsi la contrainte standard basée sur la plus grande valeur propre du système.

86
Chapter 6

Complément - Aspects physiques -


Dérivation des équations

6.1 Dérivation des équations SW


On part des équations de Navier Stokes incompressibles sans viscosité (en d’autres termes,
on considère le cas d’un fluide parfait), communément appelées équations d’Euler incom-
pressibles. On considère l’écoulement irrotationnel. On notera dans la suite
¨ ˛
Upx, y, z, tq
V “ Vpx, y, z, tq “ ˝ V px, y, z, tq ‚ (6.1)
W px, y, z, tq

la variable d’écoulement, fonction du temps et des trois variables d’espace, à valeurs dans
R3 . Dans ce contexte, rappelons que la condition d’incompressibilité s’écrit:

divpVq “ 0 , (6.2)

et que l’évolution de la vitesse est donnée par:


BV 1 ~
` divpV b Vq “ ~g ´ grad p, (6.3)
Bt ρ

où ~g “ p0, 0, ´gqT , g étant la constante de pesanteur. Le système précédent s’écrit en


formulation "composante" comme suit:
$

’ BU BU 2 BUV BUW 1 Bp

’ ` ` ` “´ ,

’ Bt Bx By Bz ρ Bx
&
BV BUV BV 2 BV W 1 Bp
` ` ` “´ , (6.4)

’ Bt Bx By Bz ρ By

’ BW BUW BV W BW 2 1 Bp


% ` ` ` “ ´g ´ .
Bt Bx By Bz ρ Bz
D’un point de vue pratique, la résolution numérique de ces équations est extrêmement coû-
teuse. D’une part, il s’agit d’un système à trois inconnues, fonctions des trois dimensions
87
d’espace et du temps. Par ailleurs, du fait que le système décrit “seulement” l’évolution
du champ de vitesse dans le fluide, on ne dispose pas d’information directe sur le domaine
de calcul, qui peut varier avec le temps en fonction des déformations de la surface libre, et
devient de facto une inconnue supplémentaire du problème. Depuis plusieurs décennies, de
nombreux travaux ont été menés afin de se ramener à des équations plus faciles à étudier et à
traiter numériquement. Une approche classique consiste à intégrer ces équations sur la pro-
fondeur, afin de réduire la dimension du problème en éliminant la variable verticale. Cette
technique doit s’adjoindre de conditions aux limites appropriées et de certaines hypothèses
simplificatrices sur les régimes d’écoulement considérés. L’un des modèles les plus classiques
reposant sur cette approche est le modèle de Saint-Venant (Shallow Water en anglais),
très couramment utilisé pour décrire les écoulements de faible profondeur relative.

6.1.1 Conditions aux limites


Commençons par préciser les conditions aux limites qui vont intervenir dans la dérivation
des équations. Dans la suite, nous noterons Zpx, y, tq l’élévation totale de la surface libre
et hpx, y, tq la hauteur d’eau. On se donne une représentation des déformations du fond
b : px, yq ÞÝÑ bpx, yq supposée de classe C 1 , de sorte que Z ˆ
“ h ` b (Figure ˙ 1). A noter que
Upx, y, z, tq
b ne dépend pas de t. On notera aussi U “ Upx, y, z, tq “ le vecteur vitesse
V px, y, z, tq
horizontal, constitué des deux premières composantes de V (6.1).

Zpx, y, tq
hpx, y, tq Upx, y, z, tq

bpx, yq
yb x
Figure 1: Représentation des différentes variables d’écoulement.

Condition de non pénétration


Il s’agit de traduire le fait que le fluide ne peut pas traverser la surface de fond. La vitesse
d’une particule au fond est donc nulle dans la direction perpendiculaire à la surface. Consid-
érons un point Mb de coordonnées px0 , y0 , bpx0 , y0 qq de cette surface. Pour toute trajectoire
d’une particule à la surface du fond, donnée par φ : ζ ÞÝÑ pxpζq, ypζq, bpxpζq, ypζqqT telle
ˆ ˙T
Bb Bb
que φp0q “ Mb , on vérifie que le vecteur ~n “ ´ px0 , y0 q , ´ px0 , y0q , 1 satisfait
Bx By

88
~n.φ1 p0q “ 0, ce qui signifie que ~n est un vecteur normal à surface au point Mb . La condition
V.~n “ 0 au point Mb s’écrit alors:
~
Wb “ Ub .gradpbq , (6.5)

où l’on a formellement noté Wb “ W px0 , y0 , bpx0 , y0 q, tq et Ub “ Upx0 , y0 , bpx0 , y0 q, tq.

Condition cinématique à la surface libre


Cette condition exprime le fait qu’une particule située à la surface libre reste à la surface
libre. Considérons donc la trajectoire d’une particule φ : t ÞÝÑ pxptq, yptq, zptqqT à la surface,
caractérisée par:
zptq “ Zpxptq, yptq, tq . (6.6)
ˆ ˙T
dxptq dyptq
La vitesse horizontale UZ de cette particule est donnée par UZ “ , . On
dt dt
dzptq
déduit alors de (6.6) l’expression de la composante verticale de la vitesse WZ “ à la
dt
surface:
BZ ~
WZ “ ` UZ .gradpZq . (6.7)
Bt

6.1.2 Adimensionnement des équations


Contexte
Nous nous intéressons ici à des régimes de faible profondeur relative, c’est à dire pour lesquels
la longueur d’onde λ est grande par rapport à la profondeur caractéristique h0 , ou encore
lorsque le paramètre µ “ h0 {λ est très petit (µ ăă 1) (voir Figure 2). A ce titre, le modèle
de Saint-Venant, qui repose sur cette approximation, fournit une bonne représentation des
mécanismes à proximité des côtes (λ est de l’ordre de la dizaine de mètres et la profondeur
de l’ordre du mètre). Toutefois, ces équations peuvent aussi être utilisées pour étudier
des phénomènes océanographiques à grande échelle, comme les tsunamis par exemple. En
effet, ces derniers se propagent sur des profondeurs de quelques kilomètres, alors que les
longueurs d’onde en jeu peuvent atteindre plusieurs centaines de kilomètres. Ces régimes de
propagation sont caractérisés par une faible variation des grandeurs de l’écoulement selon
la coordonnée verticale, ce qui, comme nous allons le voir, va se traduire par certaines
hypothèses d’échelle dans l’adimensionnement du problème.

Adimensionnement
On introduit les quantités sans dimension suivantes:
h0 x y z a
µ“ x̃ “ ỹ “ z̃ “ t̃ “ µt g{h0 ,
L L L h0
U V W p
Ũ “ ? Ṽ “ ? W̃ “ ? p̃ “ .
gh0 gh0 µ gh0 ρgh0

89
λ

h0

Figure 2: Longueur d’onde et profondeur caractéristique.

Notons que l’on retrouve ici les facteurs d’échelle précédemment mentionnés sur la variable
verticale z (adimensionnée
? par h0 “ µL) et la composante verticale de la vitesse w (adimen-
sionnée par µ gh0 ). Rappelons ici que le paramètre µ est supposé petit: dans la suite, nous
négligerons les termes d’ordre µ2 dans la dérivation des équations. On vérifie aisément que
la condition d’incompressibilité sur les variables adimensionnées s’écrit:

divpṼq “ 0 , (6.8)

et que le système (6.4) donne:


$

’ B Ũ B Ũ 2 B Ũ Ṽ B Ũ W̃ Bp

’ ` ` ` “´ ,

’ B t̃ B x̃ B ỹ B z̃ B x̃

& B Ṽ 2
B Ũ Ṽ B Ṽ B Ṽ W̃ B p̃
` ` ` “´ , (6.9)

’ B t̃˜ B x̃ B ỹ B z̃ B ỹ¸

’ 2

’ 2 B W̃ B Ũ W̃ B Ṽ W̃ B W̃ B p̃

% µ ` ` ` “ ´1 ´ .
B t̃ B x̃ B ỹ B z̃ B z̃

On omettra les “tildes” dans la suite pour raisons de simplicité d’écriture. Les conditions aux
limites adimensionnées ont alors exactement la même expression que celles vues précédem-
ment (6.5), (6.7).

6.1.3 Équation de la masse


On intègre (6.8) sur la profondeur de la colonne d’eau 1

ż z“Z ż z“Z ż z“Z


BW
divpVqdz “ divpUqdz ` dz
z“b z“b z“b Bz
ˆż Z ˙
“ div ~
U dz ´ UZ .gradpZq ~
` Ub .gradpbq ` WZ ´ Wb ,
b
ż ż vpxq
d vpxq B
1
Ce résultat provient de la formule de Leibniz f px, tqdt “ f px, tqdt ` v 1 pxqf pvpxqq ´
dx upxq upxq Bx
u1 pxqf pupxqq pour toute fonction f continue et toutes fonctions u, v dérivables, qui se généralise en dimension
supérieure.

90
où l’on a gardé les notations formelles de la Section 6.1.1 pour définir les valeurs de la vitesse
à la surface (UZ , WZ ) et au fond (Ub , Wb ). Nous noterons dans la suite
ˆ ˙ ż
u 1 z“Z
u“ “ U dz
v h z“b

la valeur moyenne de U sur la verticale. En utilisant les conditions aux limites (6.5) et
BZ Bh
(6.7), et le fait que “ (en effet Z “ h ` b et b est indépendant du temps) on obtient
Bt Bt
l’équation suivante, communément appelée équation de la masse:
Bh
` divphuq “ 0 . (6.10)
Bt

6.1.4 Equation de quantité de mouvement


On intègre à présent (6.9) sur la profondeur de la colonne d’eau. La première équation donne:
żZ żZ żZ żZ żZ
BU BU 2 BUV BUW Bp
dz ` dz ` dz ` dz ` dz “ 0 . (6.11)
b Bt b Bx b By b Bz b Bx
Le premier terme s’écrit:
żZ ˆż Z ˙
BU B BZ
dz “ Udz ´ UZ
b Bt Bt b Bt (6.12)
Bhu BZ
“ ´ UZ .
Bt Bt
Les termes faisant intervenir les composantes horizontales de la vitesse donnent:
żZ 2 ˆż Z ˙
BU B BZ Bb
dz “ U dz ´ UZ2
2
` Ub2 ,
b Bx Bx b Bx Bx
żZ ˆż Z ˙
BUV B BZ Bb
dz “ UV dz ´ UZ VZ ` Ub Vb ,
b By By b By By

de sorte que
żZ żZ ˆż Z ˙ ˆż Z ˙
BU 2 BUV B 2 B
dz ` dz “ U dz ` UV dz
b Bx b By Bx b By b (6.13)
~
´ UZ pUZ .gradpZqq ~
` Ub pUb .gradpbqq .

D’autre part :
żZ
BUW
dz “ WZ UZ ´ Wb Ub
b Bz (6.14)
BZ ~ ~
“ UZ ` UZ pUZ .gradpZqq ´ Ub pUb .gradpbqq ,
Bt

91
en vertu des conditions aux limites (6.5) et (6.7). Au final, en utilisant (6.12), (6.13) et
(6.14), l’équation (6.11) se réécrit:
ˆż Z ˙ ˆż Z ˙ żZ
Bhu B 2 B Bp
` U dz ` UV dz ` dz “ 0 . (6.15)
Bt Bx b By b b Bx

Il reste alors à traiter les termes advectifs et les termes de pression.

Termes advectifs
Première approximation: Dans la littérature, on considère parfois simplement
ż Z une vitesse
constante selon la verticale, de sorte que u “ U. On a alors immédiatement U 2 dz “ hu2
żZ b

et UV dz “ huv. Cette hypothèse peut paraître brutale, mais en réalité, elle exprime le
b
fait que l’on peut négliger les variations de vitesse sur l’axe vertical. Ce point est détaillé
plus précisément ci dessous:
Deuxième approximation: De manière plus générale, on formule une hypothèse d’écoulement
faiblement cisaillés. En d’autres termes, on autorise une déviation de la vitesse selon l’axe
vertical, mais suffisamment faible pour être négligée dans le régime asymptotique considéré
(Figure 3):
Upx, y, z, tq “ upx, y, tq ` µu1 px, y, z, tq . (6.16)
żZ
Notons que par définition de upx, y, tq, on a u1 dz “ 0. Ainsi:
b
żZ żZ
U 2 dz “ pu ` µu1 q2 dz
b b
żZ żZ żZ
(6.17)
“ 2
u dz ` µ 2 1 2
pu q dz ` 2µu u1 dz
b b b
“ hu2 ` Opµ2 q .
żZ
On vérifie de même que UV dz “ huv ` Opµ2 q. Les termes résiduels sont de l’ordre Opµ2 q
b
et peuvent ainsi être négligés, conformément à ce qui a été énoncé au début de la Section
6.1.2.

Terme de pression
Bp
La troisième équation de (6.9) donne “ ´1 ` Opµ2 q. Après intégration :
Bz
ppZq ´ ppzq “ z ´ Z ` Opµ2 q .

Remarque 6.1.1. En revenant aux équations dimensionnées, et en négligeant le terme


d’ordre Opµ2 q, cette équation s’écrit
Bp
“ ´ρg .
Bz
92
z

upx, tq
u1 px, z, tq

x
Figure 3: Illustration 1d de la décomposition (6.16) en une vitesse moyenne u et une dévia-
tion de vitesse horizontale u1 .

C’est ce que l’on appelle hypothèse de pression hydrostatique, à savoir que l’on suppose
une distribution linéaire de la pression selon la verticale. Elle traduit le fait que la pression
en un point px, yq dépend uniquement du poids de la colonne d’eau au dessus de ce point.

En considérant nulle la pression atmosphérique à la surface ppZq “ 0, on a immédiatement


ppbq “ h ` Opµ2 q, et le calcul donne:
żZ
1
ppzqdz “ h2 ` Opµ2 q .
b 2
On écrit alors:
żZ ż
Bp B Z BZ Bb
dz “ ppzqdz ´ ppZq ` ppbq
b Bx Bx b Bx Bx
ˆ ˙ (6.18)
B 1 2 Bb
“ h ` h ` Opµ2 q .
Bx 2 Bx

En considérant dans (6.15) le travail effectué sur les termes advectifs (6.17) et de pression
(6.18), on aboutit à l’équation suivante:
ˆ ˙
Bhu B 2 1 2 B Bb
` hu ` h ` huv “ ´h , (6.19)
Bt Bx 2 By Bx

et par des considérations similaires:


ˆ ˙
Bhv B B 1 2 Bb
` phuvq ` huv ` h “ ´h . (6.20)
Bt Bx By 2 By

Il s’agit des équations de quantité de mouvement. Le système final se compose des


équations (6.10) et (6.19),(6.20), et peut se reformuler (en version dimensionnée) sous la
forme compacte:
Bh
` divphuq “ 0 ,
Bt ˆ 2˙ (6.21)
Bhu ~ h ~
` divphu b uq ` grad g “ ´ghgradpbq .
Bt 2

93
6.1.5 Équation d’énergie
1
On définit l’énergie cinétique du système par K “ h}u}2 et l’énergie potentielle par E “
2
1 2
gh ` ghz. L’énergie totale est définie par E “ K ` E, et vérifie elle aussi une équation de
2
conservation. Pour la déterminer, on s’intéresse d’abord à l’évolution du champ de vitesse.
Notons tout d’abord qu’en utilisant l’équation de la masse, on a:
Bhu Bu Bh Bu
“h `u “h ´ udivphuq . (6.22)
Bt Bt Bt Bt
~
Par suite, en utilisant l’égalité divphu b uq ´ udivphuq “ [Link], les équations (6.21)
et (6.22) donnent l’équation sur la vitesse:
Bu ~ ~ phq “ ´g gradpbq
~
` [Link] ` g grad . (6.23)
Bt
L’évolution de l’énergie cinétique étant donnée par:
ˆ ˙
BK B 1 1 Bh Bu
“ h}u} “ }u}2
2
` hu. ,
Bt Bt 2 2 Bt Bt
l’invocation combinée de (6.23) et de l’équation de la masse permet de déduire l’équation
suivante: ˆ ˙
BK h}u}2 ~ ~
` div u “ ´[Link] ´ [Link] . (6.24)
Bt 2
Par suite, en revenant à la définition de E, et en utilisant une nouvelle fois l’équation de la
masse, on a:
BE Bh
“ gph ` bq “ ´gph ` bqdivphuq ,
Bt Bt
ce qui peut se réécrire
BE ~
` divpghuph ` bqq “ [Link] ` bq . (6.25)
Bt
Finalement, en additionnant les contributions (6.24) et (6.25):
ˆˆ ˙ ˙
BE h2
` div E `g `K u “ 0. (6.26)
Bt 2
Nous verrons ultérieurement que l’énergie totale E joue le rôle d’une entropie du système
(6.21).

6.2 Analyse linéaire


6.2.1 Relation de dispersion
D’un point de vue général, une grande partie des dynamiques dominantes des modèles
d’écoulement sont contenues dans les équations linéarisées. A ce titre, l’analyse linéaire
94
est un outil puissant permettant d’extraire des informations importantes sur les propriétés
et caractéristiques physiques des modèles en considération. Dans le cadre des équations de
Saint-Venant, l’étude du système linéaire est notamment utilisée pour les dynamiques grande
échelle, pour lesquelles les effets de la rotation terrestre ont une importance capitale. Au
prix de légères adaptations, il est tout à fait possible de reproduire la dérivation effectuée
ci-dessus en considérant l’effet dû à la rotation dans (6.3):

BV 1 ~
` divpV b Vq “ ~g ´ grad p ´ f~e3 ˆ V , (6.27)
Bt ρ

où f désigne la force de Coriolis. On aboutit alors au système de Saint-Venant (6.21) avec


rotation:
Bh
` divphuq “ 0 ,
Bt ˆ 2˙ (6.28)
Bhu ~ h ~
` divphu b uq ` grad g “ ´ghgradpbq ´ f huK ,
Bt 2
ˆ ˙
K ´v
où u “ . Afin d’étudier les dynamiques linéaires du système, nous allons considérer
u
des solutions perturbées autour d’un état d’équilibre au repos (vitesse nulle):

hpx, y, tq “ h0 ` h1 px, y, tq ,
(6.29)
upx, y, tq “ u1 px, y, tq .

En injectant dans (6.28), et en négligeant les termes d’ordre 2 (correspondant au produit de


perturbations), on aboutit au système linéarisé:

Bh1
` h0 divpu1 q “ 0 ,
Bt (6.30)
Bu1 ~ 1 “ ´f h0 u1K .
` g gradh
Bt
On cherche alors des solutions du système sous forme d’ondes planes, c’est à dire (on omettra
les “primes”):

hpx, y, tq “ h̄eipk.~x´ωtq ,
upx, y, tq “ ūeipk.~x´ωtq , (6.31)
vpx, y, tq “ v̄eipk.~x´ωtq ,

où h̄, ū, v̄ sont des réels. Dans l’écriture précédente, k “ pkx , ky qT est le vecteur d’onde
et ω la pulsation. Le vecteur d’onde quantifie le nombre d’ondes par unités d’espace, ou
encore le nombre d’oscillations contenues dans une longueur d’onde λ. Le vecteur d’onde est
relié à la longueur d’onde via la relation }k} “ 2π{λ. La pulsation ω, ou encore fréquence
angulaire, quantifie le nombre d’ondes par unités de temps. Cette quantité est reliée à la

95
période T par la relation classique ω “ 2π{λ. En injectant (6.31) dans (6.30), on aboutit à
un système matriciel: ¨ ˛¨ ˛ ¨ ˛
´ω kx ky h̄ 0
˝gh0 kx ´ω if ‚˝ū‚ “ ˝0‚ . (6.32)
gh0 ky if ´ω v̄ 0
Ce système admet des solutions non triviales si et seulement si le déterminant est nul, ce qui
donne l’équation: ` ` ˘˘
ω ω 2 ´ f 2 ` }k}2 gh0 “ 0 , (6.33)
qui admet pour solutions:
a
ω“0 , ω“˘ }k}2 c2 ` f 2 , (6.34)
?
où l’on a posé c “ gh. Cette relation, qui relie la pulsation au nombre d’onde, est com-
munément appelée relation de dispersion.

6.2.2 Ondes de gravité


On se place dans le cas f “ 0. Dans ce cas, les solutions d’intérêt sont données par:

ω “ ˘}k}c . (6.35)

Les ondes vérifiant cette relation sont généralement appelées ondes de gravité, du fait que
leur évolution est essentiellement pilotée par la constante de gravité g. Une conséquence
importante est que la vitesse de phase est indépendante de la longueur d’onde, et égale
à la vitesse de groupe2 . En d’autres termes, toutes les ondes se propagent à la même
vitesse dans le milieu, indépendamment de leur longueur d’onde. On dit alors que le milieu
de propagation est non-dispersif.

Remarque 6.2.1. Les ondes de gravité régissent essentiellement la propagation des tsunamis
et sont à ce titre souvent utilisées pour estimer leur vitesse en première approximation. Par
exemple, un tsunami
? généré à une profondeur de h0 “ 3000m se déplacera à une vitesse ap-
proximative de gh0 « 172m.s´1 , soit à peu près 620km{h, indépendamment de sa longueur
d’onde.

6.2.3 Ondes inertie-gravité


On s’intéresse à présent au cas f ‰ 0. La première solution de (6.34) ω “ 0, correspond
à un mode stationnaire caractérisé par l’équilibre géostrophique (obtenu à partir de la
deuxième équation (6.32) en considérant les dérivées en temps nulles).
ˆ 2˙
~ h
grad g ` f huK “ 0 . (6.36)
2
ω
2
En 1d, la vitesse de phase d’une onde monochromatique eipkx´ωtq est définie par vφ “ . La vitesse de
k
dωpkq
groupe d’un paquet d’onde est définie par la relation vg “ . Dans ce contexte, on a vϕ “ vg “ c.
dk
96
Cet équilibre entre les forces de pression et l’effet de rotation terrestre a une importance
prépondérante en météorologie et en océanographie grande échelle. Les solutions restantes
a
ω “ ˘ }k}2 c2 ` f 2 , (6.37)

correspondent à des modes rapides, et sont couramment appelées ondes de Poincaré ou


ondes inertie-gravité. En considérant
c la solution positive dans le cas 1d, on voit que la
2
f
vitesse de phase vϕ “ ω{k “ ` c2 n’est alors plus indépendante du nombre d’onde :
k2
On parle alors de milieu dispersif. A noter que la vitesse de groupe est cette fois donnée
Bω kc2
par vg “ “a . Deux cas limites sont à considérer:
Bk f 2 ` c2 k 2
Limite ondes courtes
Considérons le cas des ondes courtes, caractérisées par λ ăă c{f . En considérant la relation
λ “ 2π{}k}, on a f ăă c}k}, de sorte que la relation de dispersion peut être estimée par
ω 2 “ }k}2 c2 . En d’autres termes, les ondes courtes ne voient pas l’effet dû à la rotation
terrestre et leur comportement s’aligne sur les ondes de gravité.

Limite ondes longues Si maintenant on suppose λ ąą c{f , ou de manière équivalente


f ąą }k}c, la relation de dispersion peut-être approchée par ω 2 “ f 2 . Ces ondes de très
grande échelle sont purement inertielles et ne voient pas les effets de la gravité.

Le ratio c{f , qui détermine la compétition entre les effets inertiels et gravitaires, est ap-
pelée rayon de déformation de Rossby.

97

Vous aimerez peut-être aussi