Cours Lahm
Cours Lahm
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
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
2
Chapter 1
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 .
Ωptq
x1 ptq x2 ptq
Figure 1: Volume matériel Ω se déplaçant à la vitesse du fluide.
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
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
Remarque 1.1.2.
‚ 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:
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
pente
apu0 pyqq
y x
px, tq
b
y “ x ´ at x
u0 pxq u0 px ´ ctq
ct
‚ Si y ď 0, alors xptq “ y ` t.
‚ 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
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.
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
Par suite, en définissant u par upx, tq “ u0 pypx, tqq pour tout px, tq P R ˆ R` , on obtient:
“ 0.
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
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:
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 .
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.
Remarque 1.3.1.
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
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
‚ 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.
13
t
Ω´ Γ
~n
nt
Ω`
nx
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
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 ,
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.
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.
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 .
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
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
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
‚ 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.
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).
Soit u P L8loc pRˆs0, T rq de classe C en dehors d’une courbe régulière Γ. Alors u est
1
‚ Le long de Γ on a l’inégalité:
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 .
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)
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:
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
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 .
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:
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 .
rf puqs
´ p2a ´ 1q |u` ´ u´ | ` sgnpu` ´ u´ q pf pu` q ` f pu´ q ´ 2f pkqq ď 0 ,
rus
ou encore:
22
f puq non admissible f puq
admissible
u´ u´ u` u u` u´ u´ u
Proof.
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.
fc “ sup g.
g convexe sur I
gďf sur I
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).
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).
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
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.
˚ ˚ 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 ,
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
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é:
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
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
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:
31
L’erreur de troncature est donnée par:
∆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:
pour tout ∆t, ∆x (avec ∆t{∆x fixé), tout n P N, et toute donnée initiale pu0j q.
Proposition 2.2.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
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
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
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.
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.
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).
∆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.
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:
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:
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
ż ż
ϕ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`
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
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:
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
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.
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
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.
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
∆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
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
∆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
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.
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 :
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 ,
Proposition 2.6.3
Si un schéma conservatif est monotone, alors l’opérateur H vérifie le principe du
maximum discret:
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
jPZ
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
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:
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,
ż żΩ ż Ω Ω
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:
iv) }u∆ p., tq ´ v∆ p., tq}L1 pRq ď }u∆ p., 0q ´ v∆ p., 0q}L1 pRq ,
‚ }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 .
et on montre de même que u est minorée par m “ minjPZ u0j . On en déduit la stabilité L8 .
1
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
}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 ,
ÿ
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
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
53
2.6.2 Schémas entropiques
Définition 2.6.3: Entropie et consistance
‚ 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.
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
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.
55
Chapter 3
Dans la suite, on supposera que la solution vit dans un ensemble convexe d’état admissibles,
noté Ω.
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)
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
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.
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:
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)
Σ´ Σ`
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:
Γ
Σ´
wG
~n´
Σ`
wD
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
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
Une démonstration rapide consiste à dériver (3.17) par rapport au temps (au sens des dis-
tributions)2 , pour obtenir:
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)
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
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`
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 .
On dit que le k ième champ caractéristique est vraiment non linéaire s’il vérifie:
On dit que le k ième champ caractéristique est linéairement dégénéré s’il vérifie:
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é.
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:
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:
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.
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:
66
Chapter 4
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)
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).
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)
70
u wD b
D2 t
λ2 pw˚ q λ2 pwD q
σG
wG b C1
w˚
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˚
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:
72
u wD b
D2 t
λ1 pwG q λ1 pw˚q λ2 pw˚ q λ2 pwD q
w˚
D1
b
w˚
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).
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 ).
¨ ˛ ¨ ˛ ¨ ˛
hG h˚ ¨ ˛ hD
˚
wG ˝huG ‚ ˚ ˝
wG hu˚ ‚ h
wD ˝huD ‚
hvG h˚ vG
˚
wD ˝ hu˚ ‚
b
h˚ vD hvD
b b
1-onde b
σ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.
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
à 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
∆t w˚
wL wR
´∆x{2 λ´ ∆t λ` ∆t ∆x{2 x
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
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
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.
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.
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
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:
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
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̄ ,
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
86
Chapter 6
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)
Zpx, y, tq
hpx, y, tq Upx, y, z, tq
bpx, yq
yb x
Figure 1: Représentation des différentes variables d’écoulement.
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)
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
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)
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).
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
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
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 .
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 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
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).
BV 1 ~
` divpV b Vq “ ~g ´ grad p ´ f~e3 ˆ V , (6.27)
Bt ρ
hpx, y, tq “ h0 ` h1 px, y, tq ,
(6.29)
upx, y, tq “ u1 px, y, tq .
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.
ω “ ˘}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.
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