- Email: [email protected]
Convergence d'un schéma volumes finis explicite en temps pour les systèmes hyperboliques linéaires symétriques en domaines bornés
C. R. Acad. Sci. Paris, t. 331, Série I, p. 95–100, 2000 Analyse numérique/Numerical Analysis
Convergence d’un schéma volumes finis explicite en temps pour les systèmes hyperboliques linéaires symétriques en domaines bornés Yves COUDIÈRE a , Jean-Paul VILA b , Philippe VILLEDIEU c a b c
Université de Provence, 39, rue F.-Joliot-Curie, 13453 Marseille cedex 13, France MIP, UMR CNRS 5640, Insa, 135, avenue de Rangeuil, 31077 Toulouse cedex 4, France Onera, centre de Toulouse, 2, avenue Edouard-Belin, 31055 Toulouse cedex, France
(Reçu le 23 mars 2000, accepté le 27 avril 2000)
Résumé.
En utilisant la technique développée dans [7], nous donnons une estimation d’erreur, d’ordre 1/2 par rapport à la taille du maillage, pour une approximation volumes finis explicite en temps, pour la résolution des systèmes hyperboliques linéaires symétriques. Le résultat s’applique à des maillages non structurés quelconques. Une attention particulière est portée sur la discrétisation des conditions aux limites, afin d’assurer la stabilité du schéma. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
Convergence of a finite volume time-explicit scheme for symmetric linear hyperbolic systems on bounded domains Abstract.
Using the technique of [7], we prove an error estimate, of order 1/2 with respect to the size of the meshes, for a finite volume time explicit approximation of symmetric linear hyperbolic systems. Our result apply to any unstructured mesh. Special attention is paid to to the discretization of the boundary condition necessary for the scheme to be stable. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
Abridged English version In this Note, we consider the approximate solution of linear symmetric hyperbolic systems: d P (An − M )u = 0 (]0, T [×∂Ω), Ai ∂i u = 0 ]0, T [×∂Ω , ∂t u + u(0, ·) = u0 (Ω). i=1
(1)
For sake of clarity, we suppose that Ω is a bounded polyhedral subset of Rd , the Ai are m × m symmetric Pd matrices with constant coefficients, and An = i=1 ni Ai , where n = (n1 , . . . , nd ) is the unit normal to ∂Ω outward to Ω. We suppose that M is a m × m matrix valued function, constant along the faces of ∂Ω, such that: M + t M > 0, ker(An − M ) + ker(An + M ) = Rm . Note présentée par Olivier P IRONNEAU. S0764-4442(00)00502-4/FLA 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés.
95
Y. Coudière et al.
Under appropriate smoothness conditions [6], (1) has a unique solution u ∈ C0 ([0, T ], H1 (Ω)) ∩ C1 ([0, T ], L2(Ω)). In that case, we prove that the piecewise constant finite volume approximation vh of u, for an explicit discretization in time, verifies: Z t (Ph vh )(M − An)Ph vh dσ(x) < c h1/2 , (2) ku − vh kL2 (]0,T [×Ω) < c h1/2 , ]0,T [×∂Ω
under a CFL condition on the time step (Ph is a projection on ker(An + M )). Finite element approximations of (1), using the discontinuous Galerkin method introduced in [4] on space-time triangulations, are well known to converge with an error estimate similar to (2) (for a piecewise constant approximation) [3]. These approximations are not explicit in time. A recent adaptation of the method has been proposed [2], that yields an explicit scheme. It is achieved on space-time meshes, under restrictive assumptions on the shape of the elements. We use a different approach. The space discretization is of finite volumes type, and the time integration is performed by using the classical explicit Euler scheme. The result of convergence (2) is established without particular assumptions on the mesh. Our proof is based on [7]. After proving some stability results (Lemma 4.1), for a small enough time step (assumption (7)), we compare vh and u in terms of weak errors of consistency, from which we deduce (2). The key ingredient of the stability and convergence results is that the dissipation of energy along the boundary ∂Ω can be controlled. This is not possible in general, unless a suitable term is added to the numerical fluxes on ∂Ω.
1. Introduction On considère ici l’approximation des systèmes hyperboliques symétriques linéaires suivant : d P (An − M )u = 0 (]0, T [×∂Ω), Ai ∂i u = 0 ]0, T [×∂Ω , ∂t u + u(0, ·) = u0 (Ω). i=1 Pour clarifier l’exposé nous supposons que Ω est un ouvert polyédral borné de Rd , les Ai sont des matrices P m × m, symétriques, à coefficients constant, An = di=1 ni Ai , où n = (n1 , . . . , nd ) est la normale unitaire sortante le long de ∂Ω. On suppose que M est une matrice m × m constante par morceaux sur les faces de ∂Ω, telle que : M + t M > 0, ker(An − M ) + ker(An + M ) = Rm . Sous des conditions de régularité sur ∂Ω et sur u0 [6], on sait que ce problème admet une solution unique u ∈ C0 ([0, T ], H1 (Ω)) ∩ C1 ([0, T ], L2 (Ω)). Dans ce cas, on démontre, sous une condition de CFL, l’estimation d’erreur (2), d’ordre 1/2 par rapport à la taille du maillage, entre la solution approchée vh et la solution exacte u. L’approximation de (1) par la méthode des éléments finis discontinus de Galerkin a été largement étudiée, mais elle conduit à des schémas implicites en temps [4,3]. À notre connaissance, il n’existe aucun résultat de convergence d’un schéma explicite, mis à part un résultat récent [2] ; mais obtenu à partir de maillages espace temps, avec des restrictions sur la forme des éléments. Notre approche est différente. Nous utilisons une discrétisation volumes finis en espace, puis un schéma d’Euler explicite en temps. Nous démontrons le résultat de convergence (2) sans hypothèse particulière sur le maillage. La preuve est basée sur des idées de [7]. La difficulté supplémentaire est d’obtenir un contrôle suffisant de la dissipation d’énergie au bord ∂Ω. Pour cela, on doit introduire un terme de pénalisation dans l’expression du flux numérique sur les faces frontières.
96
Convergence d’un schéma volumes finis explicite en temps ...
2. Définition du schéma numérique Nous désignons par Th un maillage polygonal quelconque de Ω (⊂ Rd ), de taille h ; c’est-à-dire une partition de Ω en polygones tel qu’il existe des constantes a, b, c, telles que : ∀ K ∈ Th ,
m(K) > a hd ,
m(∂K) 6 b hd−1 ,
diam(K) 6 h.
(3)
m(K) et m(∂K) désignent les mesures de Lebesgue d’un polygone K ∈ Th et de son bord ∂K. diam(K) désigne le diamètre de K. Nous appelons ∆t ∈ R+ le pas de temps, et notons tn = n∆t les instants de discrétisation. Nous approchons la solution exacte de (1) par la fonction vh constante par morceaux telle que, pour n n . Les valeurs vK (K ∈ Th , et tn+1 6 T ) sont déterminées par le tout (t, x) ∈ [tn , tn+1 [ ×K, vh (t, x) = vK schéma explicite : Z ∆t P n 1 n+1 n 0 = vK − gKe m(e), vK = u0 (x) dx. (4) vK m(K) e∈∂K m(K) K n est le flux numérique sur la face e de K, défini ci-dessous. gKe On note Sh l’ensemble des faces intérieures du maillage Th , et Γh l’ensemble de ses faces frontières. Pour un volume de contrôle K ∈ Th , on note ∂K l’ensemble des faces e de sa frontière, et Ke son voisin le long de la face e (si e ∈ Sh ), n n n gKe = DKe vK − CKe vK si e ∈ Sh , e (5) ∀ K ∈ Th , ∀ e ∈ ∂K, n n si e ∈ Γh . gKe = DKe vK
Pour nKe = (n1 , . . . , nd ), la normale sortante à K le long de e (on note simplement ne = nKe la normale sortante à Ω le long de e si e ∈ Γh ), DKe = AnKe + et CKe = −AnKe − si e ∈ Sh , (6) ∀ K ∈ Th , ∀ e ∈ ∂K, si e ∈ Γh . DKe = 12 (Me + Ane ) + Qe Me est la valeur de M le long de e ⊂ ∂Ω. Le choix de DKe et CKe sur les faces internes correspond au choix d’un schéma décentré amont. Sur les faces frontières, on a Ane u = 12 (Me + Ane )u si u ∈ ker(Me − Ane ), d’où le choix de DKe . Qe est un terme de stabilisation, qui ne doit pas modifier la consistance du flux, d’où la condition ker(Ane − Me ) ⊂ ker(Qe ). Pour choisir Qe , soit la décomposition suivante de Rm : Ee ⊕ Ge = ker(Ane − Me ), Fe ⊕ Ge = ker(Ane + Me ), Rm = Ee ⊕ Ge ⊕ Fe , Ge = ker(Me ) ∩ ker(Ane ). On appelle Pe la projection sur Fe , parallèlement à Ee ⊕ Ge . Alors on prend : Qe =
1t 2 (Me
− Ane )Pe .
3. Résultat principal T HÉORÈME 3.1 (Convergence du schéma). – Soit T > 0, ε > 0 et (Th )h>0 une famille de maillages vérifiant les hypothèses (3). Si ∆t > 0 est tel que ( kAnKe − k ∆tm(∂K) 61−ε si e ∈ Sh ,
1 M +t M m(K) ∆tm(∂K) (7) ∀ K ∈ Th , ∀ e ∈ ∂K, e e
− Ane 6 1 − ε si e ∈ Γh , 2
2
m(K)
alors la solution discrète vh , définie par (4)–(6), converge vers u, et l’on a : Z t 1/2 (Ph vh )(M − An)Ph vh dσ 6 c h1/2 , ku − vh kL2 (]0,T [×Ω) 6 c h , ]0,T [×∂Ω
où ∀ e ∈ Γh , ∀ x ∈ K, Ph (x) = Pe . De plus, les dérivées discrètes sont contrôlées au sens suivant P P t n n n n 1/2 (N ∆t < T ) : Nn=0 . e∈∂K∩Sh (vK − vKe )CKe (vK − vKe ) m(e)∆t 6 c h K∈Th
97
Y. Coudière et al.
4. Estimation d’énergie pour la solution discrète On note que dans la version implicite du schéma (4), il n’est pas nécessaire de rajouter un terme de stabilisation. Dans ce cas le schéma s’interprète comme un cas particulier de la méthode des éléments finis discontinus de Galerkin sur le maillage espace temps constitué des éléments ]tn , tn+1 [ ×K. La stabilité du schéma résulte alors de la coercivité de la forme bilinéaire associée [3]. Il n’est plus possible d’interpréter le schéma explicite (4) comme un schéma d’éléments finis discontinus. De ce fait la stabilité est difficile à obtenir. Elle fait apparaître, d’une part, une condition sur le pas de temps, et d’autre part la nécessité d’un terme de stabilisation qui permette de contrôler la dissipation du schéma numérique sur les faces frontière. C’est le terme Qe , qui permet d’obtenir le résultat suivant : L EMME 4.1 (Stabilité). – Soit ε > 0. Si ∆t vérifie la condition (7), alors la solution discrète vh vérifie :
vh (t, ·)
L2 (Ω)
N X
6 ku0 kL2 (Ω) ,
t n n vK Qe vK
m(e)∆t 6
n=0 e∈Γh N X
X
t
n=0 e∈∂K∩Sh K∈Th
1 ku0 kL2 (Ω) , ε
1 n n n n − vK − vK vK CKe vK m(e)∆t 6 ku0 kL2 (Ω) . e e ε
Ce résultat est démontré précisément dans [1]. La difficulté principale par rapport au cas étudié dans [7] est liée à la présence des conditions aux limites. Un résultat analogue pour les équations de Maxwell et des conditions aux limites particulières est également prouvé dans [5]. On peut vérifier que le flux numérique sur la frontière métallique donné dans [5] correspond à un choix particulier de la matrice Qe . 5. Démonstration du théorème de convergence La démonstration du théorème 3.1 se fait en trois étapes. Étape 1 : on commence par démontrer un résultat de comparaison général. Si v est une fonction constante par morceaux telle que : Z
Z −
|v|2 ∂t φ dt dx − [0,T ]×Ω
Z
d Z X
v∂t ψ dt dx − [0,T ]×Ω i=1
−
t
An + t M ψ dt dσ(x) − vAi ∂i ψ dt dx + v [0,T ]×Ω [0,T ]×∂Ω 2 t
Z
d Z X i=1
Z
t
t
vAi v∂i φ dt dx +
[0,T ]×Ω
t
t
M +M vφ dt dσ(x) − [0,T ]×∂Ω 2 t
v
u0 ψ(0, ·) dx = hµ, ψi ,
Ω
(8)
Z
|u0 |2 φ(0, ·) dx 6 hν, φi , Ω
(9) m ¯ telles que ψ(T, ·) = 0 (presque et φ ∈ C1 ([0, T ] × Ω), pour toutes fonctions ψ ∈ H1 (]0, T [ × Ω) m partout), φ(T, ·) = 0, et φ > 0, et où µ et ν sont des formes linéaires continues sur H1 (]0, T [ × Ω) ¯ alors on a et C0 ([0, T ] × Ω),
ku − vk2L2 (]0,T [×Ω) 6 hν, φi − 2 hµ, φui ,
avec φ(t, x) = T − t.
(10)
Étape 2 : on montre alors que vh vérifie (8) et (9), et on exhibe les formes linéaires correspondantes, µh et νh . On a νh = νh0 − ενh1 , et νh0 , νh1 , µh sont donnés par :
98
Convergence d’un schéma volumes finis explicite en temps ...
0 νh , φ =
Z
N X t n n vh (0, ·) 2 − |u0 |2 φ(0, ·) dx + 2 vK (Ce − Qe )vK φne − φK tn+1 m(e)∆t
Ω
n=0 e∈Γh N X
+
X
t
n n n n n − vKe + vKe vK CKe vK φe − φK tn+1 m(e)∆t,
n=0 e∈∂K∩Sh K∈Th N X
1 νh , φ =
X
t
n n n n − vKe − vKe vK CKe vK φK tn+1 m(e)∆t
n=0 e∈∂K∩Sh K∈Th
+2
N X
t n n vK Qe vK φK
tn+1 m(e)∆t,
n=0 e∈Γh
hµh , ψi =
N X
n n − vKe vK CKe ψen − ψK tn+1 m(e)∆t +
t
Z t Ω
n=0 K∈Th N X
+
vh (0, ·) − u0 ψ(0, ·) dx
t nt vK Ce
N X t nt vK Qe ψen m(e)∆t. ψen − ψK tn+1 m(e)∆t −
n=0 e∈Γh
n=0 e∈Γh
Pour une fonction ϕ (= φ ou ψ), on note Z 1 ϕ(t, ·) dx, ϕK (t) = m(K) K
ϕne =
1 m(e)∆t
Z
tn+1
Z ϕ dt dσ.
tn
e
On remarque que νh1 est une forme linéaire positive liée à la dissipation numérique du schéma. L’analogue de l’estimation BV faible obtenue dans le cadre des lois de conservations non linéaires s’écrit ici hνh , 1li 6 c. Les termes µh et νh0 sont les erreurs de consistance du schéma. Étape 3 : on estime précisément les termes de reste µh et νh0 . Démonstration. – Étape 1. Pour démontrer le résultat de comparaison, le point clé est la régularité de u : d’une part on peut prendre ψ = uφ comme fonction test dans (8), et d’autre part, u est solution de (1) au sens classique (fort). L’inégalité (10) résulte de manipulations algébriques simples, en remarquant de plus que u vérifie (8) et (9) pour µ = 0 et ν = 0 [1]. Étape 2. Le calcul de µh est immédiat (intégration par partie discrète), et le calcul de νh résulte des inégalités de stabilité discrètes. n n Qe vK (φne − φK (tn+1 )) > Étape 3. Sachant que φ(t, x) = T − t, on a φne − φK (tn+1 ) = ∆t/2, d’où t vK t n 0 ; et u ∈ ker(Ane − Me ) le long de e, et ψ = φu, donc Qe ψe = 0 ; on peut donc éliminer les termes correspondant de hνh0 , φi − 2hµk , ψi. Notons n n CKe t vK − vKe si e ∈ Sh , n si e ∈ Sh , = ∆vKe CKe = Me + Me n − An , 2 si e ∈ Γ si e ∈ Γh . Pe vK e h 2 En utilisant (10), après un calcul immédiat, on trouve :
ku − vh k2L2 (]0,T [×Ω) + ε νh1 , φ 6 c T ku0k2H1 (Ω) h2 +
N X X n=0
K∈Th
e∈∂K
t
n ∆vKe CKe εnK + εnKe
N X X ∆t t n m(e)∆t + ∆vKe CKe RnKe m(e)∆t, 2 n=0 K∈Th
e∈∂K
99
Y. Coudière et al. n où εK = vK − unK (et en convenant que εnKe = 0 si e ∈ Γh ), et ( (unK + unKe )(φne − φK (tn+1 )) − 2(ψen − ψK (tn+1 )) n RKe = unK (φne − φK (tn+1 )) − (ψen − ψK (tn+1 ))
si e ∈ Sh , si e ∈ Γh .
Grâce à l’inégalité de Cauchy–Schwartz, sous la condition (7), on a N X X n=0
K∈Th
t
n ∆vKe CKe εnK + εnKe
e∈∂K
1/2 1/2 ∆t m(e)∆t 6 c νh1 , 1l ∆t ku − vh kL2 (]0,T [×Ω) . 2
n n φK (tn+1 ) + δRKe , avec On vérifie ensuite que RnKe = RKe ( n+1 R n t un 2 K +uKe − ue (s) (tn+1 − s) ds 2(uK (tn+1 ) − une ), n n n ∆t 2 t δRKe = RKe = R tn+1 n 1 (uK (tn+1 ) − une ), u − u (s) (tn+1 − s) ds K
∆t tn
e
si e ∈ Sh , si e ∈ Γh .
Grâce à des résultats de régularité classiques pour v ∈ C ([0, T ], H (Ω)) ∩ C ([0, T ], L (Ω)), on a Z tn+1 X
2
n 2
u(t, ·) 2 1 RKe m(e) 6 c h + ∂t u(t, ·) L2 (K) dt, H (K) 0
e∈∂K
1
2
tn
Z X δRn 2 m(e) 6 c h3 Ke
e∈∂K
1
tn+1
tn
2
u(t, ·) 2 1 + ∂t u(t, ·) L2 (K) dt. H (K)
On en déduit N X X n=0
K∈Th
e∈∂K
t
1/2
ε 1 1 n n ν ,φ + kuk2V h, ∆vKe CKe RKe φK tn+1 m(e)∆t 6 νh1 , φ kukV h1/2 6 2 h 2ε
N X X n=0
K∈Th
t
1/2 n n ∆vKe CKe δRKe m(e)∆t 6 νh1 , 1l kukV h3/2 ,
e∈∂K
1/2 . En regroupant les différentes contributions, on où kukV = sup[0,T ] ku(t, ·)k2H1 (Ω) + k∂t u(t, ·)k2L2 (Ω)
a finalement : ku − vh k2L2 (]0,T [×Ω) + 2ε νh1 , φ 6 c h (1 + h1/2 + h) + c h1/2 ku − vh kL2 (]0,T [×Ω) , ce qu’il fallait démontrer. Références bibliographiques [1] Coudière Y., Analyse de schémas volumes finis sur maillages non structurés pour des problèmes linéaires hyperboliques et elliptiques, Thèse de doctorat, Université Paul-Sabatier, 1999. [2] Falk R.S., Richter G.R., Explicit finite element methods for symmetric hyperbolic equations, SIAM J. Numer. Anal. 36 (3) (1999) 935–952. [3] Johnson C., Pitkaranka J., An analysis of the discontinuous galerkin method for a scalar hyperbolic equation, Math. of Comp. 47 (1986) 285–312. [4] Lesaint P., Raviart P.-A., On a finite element method for solving the neutron transport equation, in: Mathematical Aspects of Finite Element in Partial Differential Equations, De Boor C. (Ed.), Academic Press, New York, 1974, pp. 89–123. [5] Piperno S., L2 -stability of finite volumes schemes for the Maxwell system in two and three dimensions on arbitrary unstructured meshes, Modél. Math. Anal. Numér. 34 (1) (2000) 139–158. [6] Rauch J., Symmetric positive systems with boundary characteristic of constant multiplicity, Trans. Amer. Math. Soc. (1985). [7] Vila J.-P., Villedieu P., Convergence de la méthode des volumes finis pour les systèmes de Friedrichs, C. R. Acad. Sci. Paris, Série I 325 (1997) 671–676.
100
Convergence d’un schéma volumes finis explicite en temps pour les systèmes hyperboliques linéaires symétriques en domaines bornés Yves COUDIÈRE a , Jean-Paul VILA b , Philippe VILLEDIEU c a b c
Université de Provence, 39, rue F.-Joliot-Curie, 13453 Marseille cedex 13, France MIP, UMR CNRS 5640, Insa, 135, avenue de Rangeuil, 31077 Toulouse cedex 4, France Onera, centre de Toulouse, 2, avenue Edouard-Belin, 31055 Toulouse cedex, France
(Reçu le 23 mars 2000, accepté le 27 avril 2000)
Résumé.
En utilisant la technique développée dans [7], nous donnons une estimation d’erreur, d’ordre 1/2 par rapport à la taille du maillage, pour une approximation volumes finis explicite en temps, pour la résolution des systèmes hyperboliques linéaires symétriques. Le résultat s’applique à des maillages non structurés quelconques. Une attention particulière est portée sur la discrétisation des conditions aux limites, afin d’assurer la stabilité du schéma. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
Convergence of a finite volume time-explicit scheme for symmetric linear hyperbolic systems on bounded domains Abstract.
Using the technique of [7], we prove an error estimate, of order 1/2 with respect to the size of the meshes, for a finite volume time explicit approximation of symmetric linear hyperbolic systems. Our result apply to any unstructured mesh. Special attention is paid to to the discretization of the boundary condition necessary for the scheme to be stable. 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS
Abridged English version In this Note, we consider the approximate solution of linear symmetric hyperbolic systems: d P (An − M )u = 0 (]0, T [×∂Ω), Ai ∂i u = 0 ]0, T [×∂Ω , ∂t u + u(0, ·) = u0 (Ω). i=1
(1)
For sake of clarity, we suppose that Ω is a bounded polyhedral subset of Rd , the Ai are m × m symmetric Pd matrices with constant coefficients, and An = i=1 ni Ai , where n = (n1 , . . . , nd ) is the unit normal to ∂Ω outward to Ω. We suppose that M is a m × m matrix valued function, constant along the faces of ∂Ω, such that: M + t M > 0, ker(An − M ) + ker(An + M ) = Rm . Note présentée par Olivier P IRONNEAU. S0764-4442(00)00502-4/FLA 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés.
95
Y. Coudière et al.
Under appropriate smoothness conditions [6], (1) has a unique solution u ∈ C0 ([0, T ], H1 (Ω)) ∩ C1 ([0, T ], L2(Ω)). In that case, we prove that the piecewise constant finite volume approximation vh of u, for an explicit discretization in time, verifies: Z t (Ph vh )(M − An)Ph vh dσ(x) < c h1/2 , (2) ku − vh kL2 (]0,T [×Ω) < c h1/2 , ]0,T [×∂Ω
under a CFL condition on the time step (Ph is a projection on ker(An + M )). Finite element approximations of (1), using the discontinuous Galerkin method introduced in [4] on space-time triangulations, are well known to converge with an error estimate similar to (2) (for a piecewise constant approximation) [3]. These approximations are not explicit in time. A recent adaptation of the method has been proposed [2], that yields an explicit scheme. It is achieved on space-time meshes, under restrictive assumptions on the shape of the elements. We use a different approach. The space discretization is of finite volumes type, and the time integration is performed by using the classical explicit Euler scheme. The result of convergence (2) is established without particular assumptions on the mesh. Our proof is based on [7]. After proving some stability results (Lemma 4.1), for a small enough time step (assumption (7)), we compare vh and u in terms of weak errors of consistency, from which we deduce (2). The key ingredient of the stability and convergence results is that the dissipation of energy along the boundary ∂Ω can be controlled. This is not possible in general, unless a suitable term is added to the numerical fluxes on ∂Ω.
1. Introduction On considère ici l’approximation des systèmes hyperboliques symétriques linéaires suivant : d P (An − M )u = 0 (]0, T [×∂Ω), Ai ∂i u = 0 ]0, T [×∂Ω , ∂t u + u(0, ·) = u0 (Ω). i=1 Pour clarifier l’exposé nous supposons que Ω est un ouvert polyédral borné de Rd , les Ai sont des matrices P m × m, symétriques, à coefficients constant, An = di=1 ni Ai , où n = (n1 , . . . , nd ) est la normale unitaire sortante le long de ∂Ω. On suppose que M est une matrice m × m constante par morceaux sur les faces de ∂Ω, telle que : M + t M > 0, ker(An − M ) + ker(An + M ) = Rm . Sous des conditions de régularité sur ∂Ω et sur u0 [6], on sait que ce problème admet une solution unique u ∈ C0 ([0, T ], H1 (Ω)) ∩ C1 ([0, T ], L2 (Ω)). Dans ce cas, on démontre, sous une condition de CFL, l’estimation d’erreur (2), d’ordre 1/2 par rapport à la taille du maillage, entre la solution approchée vh et la solution exacte u. L’approximation de (1) par la méthode des éléments finis discontinus de Galerkin a été largement étudiée, mais elle conduit à des schémas implicites en temps [4,3]. À notre connaissance, il n’existe aucun résultat de convergence d’un schéma explicite, mis à part un résultat récent [2] ; mais obtenu à partir de maillages espace temps, avec des restrictions sur la forme des éléments. Notre approche est différente. Nous utilisons une discrétisation volumes finis en espace, puis un schéma d’Euler explicite en temps. Nous démontrons le résultat de convergence (2) sans hypothèse particulière sur le maillage. La preuve est basée sur des idées de [7]. La difficulté supplémentaire est d’obtenir un contrôle suffisant de la dissipation d’énergie au bord ∂Ω. Pour cela, on doit introduire un terme de pénalisation dans l’expression du flux numérique sur les faces frontières.
96
Convergence d’un schéma volumes finis explicite en temps ...
2. Définition du schéma numérique Nous désignons par Th un maillage polygonal quelconque de Ω (⊂ Rd ), de taille h ; c’est-à-dire une partition de Ω en polygones tel qu’il existe des constantes a, b, c, telles que : ∀ K ∈ Th ,
m(K) > a hd ,
m(∂K) 6 b hd−1 ,
diam(K) 6 h.
(3)
m(K) et m(∂K) désignent les mesures de Lebesgue d’un polygone K ∈ Th et de son bord ∂K. diam(K) désigne le diamètre de K. Nous appelons ∆t ∈ R+ le pas de temps, et notons tn = n∆t les instants de discrétisation. Nous approchons la solution exacte de (1) par la fonction vh constante par morceaux telle que, pour n n . Les valeurs vK (K ∈ Th , et tn+1 6 T ) sont déterminées par le tout (t, x) ∈ [tn , tn+1 [ ×K, vh (t, x) = vK schéma explicite : Z ∆t P n 1 n+1 n 0 = vK − gKe m(e), vK = u0 (x) dx. (4) vK m(K) e∈∂K m(K) K n est le flux numérique sur la face e de K, défini ci-dessous. gKe On note Sh l’ensemble des faces intérieures du maillage Th , et Γh l’ensemble de ses faces frontières. Pour un volume de contrôle K ∈ Th , on note ∂K l’ensemble des faces e de sa frontière, et Ke son voisin le long de la face e (si e ∈ Sh ), n n n gKe = DKe vK − CKe vK si e ∈ Sh , e (5) ∀ K ∈ Th , ∀ e ∈ ∂K, n n si e ∈ Γh . gKe = DKe vK
Pour nKe = (n1 , . . . , nd ), la normale sortante à K le long de e (on note simplement ne = nKe la normale sortante à Ω le long de e si e ∈ Γh ), DKe = AnKe + et CKe = −AnKe − si e ∈ Sh , (6) ∀ K ∈ Th , ∀ e ∈ ∂K, si e ∈ Γh . DKe = 12 (Me + Ane ) + Qe Me est la valeur de M le long de e ⊂ ∂Ω. Le choix de DKe et CKe sur les faces internes correspond au choix d’un schéma décentré amont. Sur les faces frontières, on a Ane u = 12 (Me + Ane )u si u ∈ ker(Me − Ane ), d’où le choix de DKe . Qe est un terme de stabilisation, qui ne doit pas modifier la consistance du flux, d’où la condition ker(Ane − Me ) ⊂ ker(Qe ). Pour choisir Qe , soit la décomposition suivante de Rm : Ee ⊕ Ge = ker(Ane − Me ), Fe ⊕ Ge = ker(Ane + Me ), Rm = Ee ⊕ Ge ⊕ Fe , Ge = ker(Me ) ∩ ker(Ane ). On appelle Pe la projection sur Fe , parallèlement à Ee ⊕ Ge . Alors on prend : Qe =
1t 2 (Me
− Ane )Pe .
3. Résultat principal T HÉORÈME 3.1 (Convergence du schéma). – Soit T > 0, ε > 0 et (Th )h>0 une famille de maillages vérifiant les hypothèses (3). Si ∆t > 0 est tel que ( kAnKe − k ∆tm(∂K) 61−ε si e ∈ Sh ,
1 M +t M m(K) ∆tm(∂K) (7) ∀ K ∈ Th , ∀ e ∈ ∂K, e e
− Ane 6 1 − ε si e ∈ Γh , 2
2
m(K)
alors la solution discrète vh , définie par (4)–(6), converge vers u, et l’on a : Z t 1/2 (Ph vh )(M − An)Ph vh dσ 6 c h1/2 , ku − vh kL2 (]0,T [×Ω) 6 c h , ]0,T [×∂Ω
où ∀ e ∈ Γh , ∀ x ∈ K, Ph (x) = Pe . De plus, les dérivées discrètes sont contrôlées au sens suivant P P t n n n n 1/2 (N ∆t < T ) : Nn=0 . e∈∂K∩Sh (vK − vKe )CKe (vK − vKe ) m(e)∆t 6 c h K∈Th
97
Y. Coudière et al.
4. Estimation d’énergie pour la solution discrète On note que dans la version implicite du schéma (4), il n’est pas nécessaire de rajouter un terme de stabilisation. Dans ce cas le schéma s’interprète comme un cas particulier de la méthode des éléments finis discontinus de Galerkin sur le maillage espace temps constitué des éléments ]tn , tn+1 [ ×K. La stabilité du schéma résulte alors de la coercivité de la forme bilinéaire associée [3]. Il n’est plus possible d’interpréter le schéma explicite (4) comme un schéma d’éléments finis discontinus. De ce fait la stabilité est difficile à obtenir. Elle fait apparaître, d’une part, une condition sur le pas de temps, et d’autre part la nécessité d’un terme de stabilisation qui permette de contrôler la dissipation du schéma numérique sur les faces frontière. C’est le terme Qe , qui permet d’obtenir le résultat suivant : L EMME 4.1 (Stabilité). – Soit ε > 0. Si ∆t vérifie la condition (7), alors la solution discrète vh vérifie :
vh (t, ·)
L2 (Ω)
N X
6 ku0 kL2 (Ω) ,
t n n vK Qe vK
m(e)∆t 6
n=0 e∈Γh N X
X
t
n=0 e∈∂K∩Sh K∈Th
1 ku0 kL2 (Ω) , ε
1 n n n n − vK − vK vK CKe vK m(e)∆t 6 ku0 kL2 (Ω) . e e ε
Ce résultat est démontré précisément dans [1]. La difficulté principale par rapport au cas étudié dans [7] est liée à la présence des conditions aux limites. Un résultat analogue pour les équations de Maxwell et des conditions aux limites particulières est également prouvé dans [5]. On peut vérifier que le flux numérique sur la frontière métallique donné dans [5] correspond à un choix particulier de la matrice Qe . 5. Démonstration du théorème de convergence La démonstration du théorème 3.1 se fait en trois étapes. Étape 1 : on commence par démontrer un résultat de comparaison général. Si v est une fonction constante par morceaux telle que : Z
Z −
|v|2 ∂t φ dt dx − [0,T ]×Ω
Z
d Z X
v∂t ψ dt dx − [0,T ]×Ω i=1
−
t
An + t M ψ dt dσ(x) − vAi ∂i ψ dt dx + v [0,T ]×Ω [0,T ]×∂Ω 2 t
Z
d Z X i=1
Z
t
t
vAi v∂i φ dt dx +
[0,T ]×Ω
t
t
M +M vφ dt dσ(x) − [0,T ]×∂Ω 2 t
v
u0 ψ(0, ·) dx = hµ, ψi ,
Ω
(8)
Z
|u0 |2 φ(0, ·) dx 6 hν, φi , Ω
(9) m ¯ telles que ψ(T, ·) = 0 (presque et φ ∈ C1 ([0, T ] × Ω), pour toutes fonctions ψ ∈ H1 (]0, T [ × Ω) m partout), φ(T, ·) = 0, et φ > 0, et où µ et ν sont des formes linéaires continues sur H1 (]0, T [ × Ω) ¯ alors on a et C0 ([0, T ] × Ω),
ku − vk2L2 (]0,T [×Ω) 6 hν, φi − 2 hµ, φui ,
avec φ(t, x) = T − t.
(10)
Étape 2 : on montre alors que vh vérifie (8) et (9), et on exhibe les formes linéaires correspondantes, µh et νh . On a νh = νh0 − ενh1 , et νh0 , νh1 , µh sont donnés par :
98
Convergence d’un schéma volumes finis explicite en temps ...
0 νh , φ =
Z
N X t n n vh (0, ·) 2 − |u0 |2 φ(0, ·) dx + 2 vK (Ce − Qe )vK φne − φK tn+1 m(e)∆t
Ω
n=0 e∈Γh N X
+
X
t
n n n n n − vKe + vKe vK CKe vK φe − φK tn+1 m(e)∆t,
n=0 e∈∂K∩Sh K∈Th N X
1 νh , φ =
X
t
n n n n − vKe − vKe vK CKe vK φK tn+1 m(e)∆t
n=0 e∈∂K∩Sh K∈Th
+2
N X
t n n vK Qe vK φK
tn+1 m(e)∆t,
n=0 e∈Γh
hµh , ψi =
N X
n n − vKe vK CKe ψen − ψK tn+1 m(e)∆t +
t
Z t Ω
n=0 K∈Th N X
+
vh (0, ·) − u0 ψ(0, ·) dx
t nt vK Ce
N X t nt vK Qe ψen m(e)∆t. ψen − ψK tn+1 m(e)∆t −
n=0 e∈Γh
n=0 e∈Γh
Pour une fonction ϕ (= φ ou ψ), on note Z 1 ϕ(t, ·) dx, ϕK (t) = m(K) K
ϕne =
1 m(e)∆t
Z
tn+1
Z ϕ dt dσ.
tn
e
On remarque que νh1 est une forme linéaire positive liée à la dissipation numérique du schéma. L’analogue de l’estimation BV faible obtenue dans le cadre des lois de conservations non linéaires s’écrit ici hνh , 1li 6 c. Les termes µh et νh0 sont les erreurs de consistance du schéma. Étape 3 : on estime précisément les termes de reste µh et νh0 . Démonstration. – Étape 1. Pour démontrer le résultat de comparaison, le point clé est la régularité de u : d’une part on peut prendre ψ = uφ comme fonction test dans (8), et d’autre part, u est solution de (1) au sens classique (fort). L’inégalité (10) résulte de manipulations algébriques simples, en remarquant de plus que u vérifie (8) et (9) pour µ = 0 et ν = 0 [1]. Étape 2. Le calcul de µh est immédiat (intégration par partie discrète), et le calcul de νh résulte des inégalités de stabilité discrètes. n n Qe vK (φne − φK (tn+1 )) > Étape 3. Sachant que φ(t, x) = T − t, on a φne − φK (tn+1 ) = ∆t/2, d’où t vK t n 0 ; et u ∈ ker(Ane − Me ) le long de e, et ψ = φu, donc Qe ψe = 0 ; on peut donc éliminer les termes correspondant de hνh0 , φi − 2hµk , ψi. Notons n n CKe t vK − vKe si e ∈ Sh , n si e ∈ Sh , = ∆vKe CKe = Me + Me n − An , 2 si e ∈ Γ si e ∈ Γh . Pe vK e h 2 En utilisant (10), après un calcul immédiat, on trouve :
ku − vh k2L2 (]0,T [×Ω) + ε νh1 , φ 6 c T ku0k2H1 (Ω) h2 +
N X X n=0
K∈Th
e∈∂K
t
n ∆vKe CKe εnK + εnKe
N X X ∆t t n m(e)∆t + ∆vKe CKe RnKe m(e)∆t, 2 n=0 K∈Th
e∈∂K
99
Y. Coudière et al. n où εK = vK − unK (et en convenant que εnKe = 0 si e ∈ Γh ), et ( (unK + unKe )(φne − φK (tn+1 )) − 2(ψen − ψK (tn+1 )) n RKe = unK (φne − φK (tn+1 )) − (ψen − ψK (tn+1 ))
si e ∈ Sh , si e ∈ Γh .
Grâce à l’inégalité de Cauchy–Schwartz, sous la condition (7), on a N X X n=0
K∈Th
t
n ∆vKe CKe εnK + εnKe
e∈∂K
1/2 1/2 ∆t m(e)∆t 6 c νh1 , 1l ∆t ku − vh kL2 (]0,T [×Ω) . 2
n n φK (tn+1 ) + δRKe , avec On vérifie ensuite que RnKe = RKe ( n+1 R n t un 2 K +uKe − ue (s) (tn+1 − s) ds 2(uK (tn+1 ) − une ), n n n ∆t 2 t δRKe = RKe = R tn+1 n 1 (uK (tn+1 ) − une ), u − u (s) (tn+1 − s) ds K
∆t tn
e
si e ∈ Sh , si e ∈ Γh .
Grâce à des résultats de régularité classiques pour v ∈ C ([0, T ], H (Ω)) ∩ C ([0, T ], L (Ω)), on a Z tn+1 X
2
n 2
u(t, ·) 2 1 RKe m(e) 6 c h + ∂t u(t, ·) L2 (K) dt, H (K) 0
e∈∂K
1
2
tn
Z X δRn 2 m(e) 6 c h3 Ke
e∈∂K
1
tn+1
tn
2
u(t, ·) 2 1 + ∂t u(t, ·) L2 (K) dt. H (K)
On en déduit N X X n=0
K∈Th
e∈∂K
t
1/2
ε 1 1 n n ν ,φ + kuk2V h, ∆vKe CKe RKe φK tn+1 m(e)∆t 6 νh1 , φ kukV h1/2 6 2 h 2ε
N X X n=0
K∈Th
t
1/2 n n ∆vKe CKe δRKe m(e)∆t 6 νh1 , 1l kukV h3/2 ,
e∈∂K
1/2 . En regroupant les différentes contributions, on où kukV = sup[0,T ] ku(t, ·)k2H1 (Ω) + k∂t u(t, ·)k2L2 (Ω)
a finalement : ku − vh k2L2 (]0,T [×Ω) + 2ε νh1 , φ 6 c h (1 + h1/2 + h) + c h1/2 ku − vh kL2 (]0,T [×Ω) , ce qu’il fallait démontrer. Références bibliographiques [1] Coudière Y., Analyse de schémas volumes finis sur maillages non structurés pour des problèmes linéaires hyperboliques et elliptiques, Thèse de doctorat, Université Paul-Sabatier, 1999. [2] Falk R.S., Richter G.R., Explicit finite element methods for symmetric hyperbolic equations, SIAM J. Numer. Anal. 36 (3) (1999) 935–952. [3] Johnson C., Pitkaranka J., An analysis of the discontinuous galerkin method for a scalar hyperbolic equation, Math. of Comp. 47 (1986) 285–312. [4] Lesaint P., Raviart P.-A., On a finite element method for solving the neutron transport equation, in: Mathematical Aspects of Finite Element in Partial Differential Equations, De Boor C. (Ed.), Academic Press, New York, 1974, pp. 89–123. [5] Piperno S., L2 -stability of finite volumes schemes for the Maxwell system in two and three dimensions on arbitrary unstructured meshes, Modél. Math. Anal. Numér. 34 (1) (2000) 139–158. [6] Rauch J., Symmetric positive systems with boundary characteristic of constant multiplicity, Trans. Amer. Math. Soc. (1985). [7] Vila J.-P., Villedieu P., Convergence de la méthode des volumes finis pour les systèmes de Friedrichs, C. R. Acad. Sci. Paris, Série I 325 (1997) 671–676.
100