arXiv:0704.0783v2  [math.NA]  20 Nov 2007
Mathematical Modelling and Numerical Analysis Will be set by the publisher
Modélisation Mathématique et Analyse Numérique
STABILITY OF A FINITE VOLUME SCHEME
FOR THE INCOMPRESSIBLE FLUIDS
Sébastien Zimmermann1
Abstract. We introduce a finite volume scheme for the two-dimensional incompressible Navier-Stokes
equations. We use a triangular mesh. The unknowns for the velocity and pressure are respectively
piecewise constant and affine. We use a projection method to deal with the incompressibility constraint.
We show that the differential operators in the Navier-Stokes equations and their discrete counterparts
share similar properties. In particular we state an inf-sup (Babuška-Brezzi) condition. Using these
properties we infer the stability of the scheme.
Résumé. Nous introduisons ici un schéma volumes finis pour les équations de Navier-Stokes in-
compressibles en deux dimensions. Les maillages considérés sont formés de triangles. Les inconnues
associées à la vitesse et la pression sont respectivement constantes et affines par morceaux. Nous
utilisons une méthode de projection pour traiter la contrainte d’incompressibilité. Nous vérifions que
les opérateurs différentiels apparaissant dans les équations de Navier-Stokes et leurs analogues dis-
crets vérifient des propriétés similaires. Nous prouvons en particulier une condition inf-sup. Nous en
déduisons la stabilité du schéma.
1991 Mathematics Subject Classification. 76D05, 74S10, 65M12.
Received: 13 august 2007.
1. Introduction
We consider the flow of an incompressible fluid in a polyhedral set Ω ⊂ R2 during the time interval [0, T ].
The velocity field u : Ω× [0, T ] → R2 and the pressure field p : Ω× [0, T ] → R satisfy the Navier-Stokes equations
∆u+ (u · ∇)u+∇p = f , (1.1)
div u = 0 , (1.2)
with the boundary and initial conditions
u|∂Ω = 0 , u|t=0 = u0.
The terms ∆u and (u·∇)u are associated with the physical phenomena of diffusion and convection, respectively.
The Reynolds number Re measures the influence of convection in the flow. For equations (1.1)–(1.2), finite
element and finite difference methods are well known and mathematical studies are available (see [9] for example).
Keywords and phrases: Incompressible fluids, Navier-Stokes equations, projection methods, finite volume.
1 17 rue Barrème - 69006 LYON. e-mail: Sebastien.Zimmermann@ec-lyon.fr
c© EDP Sciences, SMAI 1999
http://arxiv.org/abs/0704.0783v2
2 TITLE WILL BE SET BY THE PUBLISHER
For finite volume schemes, numerous computations have been conducted ( [12] and [1] for example). However,
few mathematical results are available in this case. Let us cite Eymard and Herbin [6] and Eymard,
Latché and Herbin [7]. In order to deal with the incompressibility constraint (1.2), these works use a
penalization method. Another way is to use the projection methods which have been introduced by Chorin [4]
and Temam [13]. This is the case in Faure [8] where the mesh is made of squares. In Zimmermann [14] the
mesh is made of triangles, so that more complex geometries can be considered. In the present paper the mesh
is also made of triangles, but we consider a different discretization for the pressure. It leads to a linear system
with a better-conditioned matrix. The layout of the article is the following. We first introduce in section 2
the discrete setting. We state (section 2.1) some notations and hypotheses on the mesh. We define (section
2.2) the spaces we use to approximate the velocity and pressure. We define also (section 2.3) the operators
we use to approximate the differential operators in (1.1)–(1.2). Combining this with a projection method, we
build the scheme in section 3. In order to provide a mathematical analysis, we show in section 4 that the
differential operators in (1.1)–(1.2) and their discrete counterparts share similar properties. In particular, the
discrete operators for the gradient and the divergence are adjoint. The discrete operator for the convection term
is positive, stable and consistent. The discrete operator for the divergence satisfy an inf-sup (Babuška-Brezzi)
condition. From these properties we deduce in section 5 the stability of the scheme.
We conclude with some notations. The spaces (L2, |.|) and (L∞, ‖.‖∞) are the usual Lebesgue spaces and we
set L20 = {q ∈ L
q(x) dx = 0}. Their vectorial counterparts are (L2, |.|) and (L∞, ‖.‖∞) with L
2 = (L2)2
and L∞ = (L∞). For k ∈ N∗, (Hk, ‖ ·‖k) is the usual Sobolev space. Its vectorial counterpart is (H
k, ‖.‖k) with
Hk = (Hk)2. For k = 1, the functions of H1 with a null trace on the boundary form the space H10. Also, we set
∇u = (∇u1,∇u2)
T if u = (u1, u2) ∈ H
1. If X ⊂ L2 is a Banach space, we define C(0, T ;X) (resp. L2(0, T ;X))
as the set of the applications g : [0, T ] → X such that t → |g(t)| is continuous (resp. square integrable). The
norm ‖.‖C(0,T ;X) is defined by ‖g‖C(0,T ;X) = sups∈[0,T ] |g(s)|. In all calculations, C is a generic positive constant,
depending only on Ω, u0 and f .
2. Discrete setting
First, we introduce the spaces and the operators needed to build the scheme.
2.1. The mesh
Let Th be a triangular mesh of Ω. The circumscribed circle of a triangle K ∈ Th is centered at xK and has
the diameter hK . We set h = maxK∈Th hK . We assume that all the interior angles of the triangles of the mesh
are less than π
, so that xK ∈ K. The set of the edges of the triangle K ∈ Th is EK . The symbol nK,σ denotes
the unit vector normal to an edge σ ∈ EK and pointing outward K. We denote by Eh the set of the edges of the
mesh. We distinguish the subset E inth ⊂ Eh (resp. E
h ) of the edges located inside Ω (resp. on ∂Ω). The middle
of an edge σ ∈ Eh is xσ and its length |σ|. For each edge σ ∈ E
h , let Kσ and Lσ be the two triangles having
σ in common. We set dσ = d(xKσ ,xLσ ). For all σ ∈ E
h , only the triangle Kσ located inside Ω is defined and
we set dσ = d(xKσ ,xσ). Then for all σ ∈ Eh we set τσ =
. As in [5] we assume the following on the mesh:
there exists C > 0 such that
∀σ ∈ Eh , dσ ≥ C |σ| and |σ| ≥ C h.
It implies that there exists C > 0 such that
∀σ ∈ E inth , τσ = |σ|/dσ ≥ C. (2.1)
2.2. The discrete spaces
We first define
P0 = {q ∈ L
2 ; ∀K ∈ Th, q|K is a constant} , P0 = (P0)
TITLE WILL BE SET BY THE PUBLISHER 3
For the sake of concision, we set for all qh ∈ P0 (resp. vh ∈ P0) and all triangle K ∈ Th: qK = qh|K (resp.
vK = vh|K). Although P0 6⊂ H
1, we define the discrete equivalent of a H1 norm as follows. For all vh ∈ P0
we set
‖vh‖h =
σ∈Eint
τσ |vLσ − vKσ |
σ∈Eext
τσ |vKσ |
. (2.2)
We have [5] a Poincaré-like inequality: there exists C > 0 such that for all vh ∈ P0
|vh| ≤ C ‖vh‖h. (2.3)
We also have [14] an inverse inequality: there exists C > 0 such that for all vh ∈ P0
h ‖vh‖h ≤ C |vh|. (2.4)
From the norm ‖.‖h we deduce a dual norm. For all vh ∈ P0 we set
‖vh‖−1,h = sup
(vh,ψh)
‖ψh‖h
. (2.5)
For all uh ∈ P0 and vh ∈ P0 we have (uh,vh) ≤ ‖uh‖−1,h ‖vh‖h. We define the projection operator ΠP0 :
L2 → P0 as follows. For all w ∈ L
2, ΠP0w ∈ P0 is given by
∀K ∈ Th , (ΠP0w)|K =
w(x) dx. (2.6)
We easily check that for all w ∈ L2 and vh ∈ P0 we have (ΠP0w,vh) = (w,vh). We deduce from this that ΠP0
is stable for the L2 norm. We define also the operator Π̃P0 : H
2 → P0. For all w ∈ H
2, Π̃P0w ∈ P0 is given by
∀K ∈ Th , Π̃P0w|K = w(xK).
According to the Sobolev embedding theorem, w ∈ H2 is a.e. equal to a continuous function. Therefore the
definition above makes sense. We introduce also the finite element spaces
P d1 = {v ∈ L
2 ; ∀K ∈ Th, v|K is affine} ,
Pnc1 = {vh ∈ P
1 ; ∀σ ∈ E
h , vh|Kσ (xσ) = vh|Lσ(xσ) ,
Pc1 = {vh ∈ (P
2 ; vh is continuous and vh|∂Ω = 0}.
We have Pc1 ⊂ H
0. We define ΠPc1 : H
0 → P
1. For all v = (v1, v2) ∈ H
0, ΠPc1v = (v
h) ∈ P
1 is given by
∀φh = (φ
h) ∈ P
∇vih,∇φ
∇vi,∇φ
The operator ΠPc
is stable for the H1 norm. One checks ( [2] p. 110) that there exists C > 0 such that for all
v ∈ H1
|v −ΠPc
v| ≤ C h ‖v‖1. (2.7)
Let us address now the space Pnc1 . If qh ∈ P
1 , we have usually ∇qh 6∈ L
2. Thus we define the operator
∇h : P
1 → P0 by setting for all qh ∈ P0 and all triangle K ∈ Th
∇hqh|K =
∇qh dx. (2.8)
4 TITLE WILL BE SET BY THE PUBLISHER
The associated norm is defined by
‖qh‖1,h =
2 + |∇hqh|
We have a Poincaré-like inequality : there exists C > 0 such that for all qh ∈ P
1 ∩ L
|qh| ≤ C |∇hqh|. (2.9)
We define the projection operator ΠPnc
. For all q ∈ H1, ΠPnc
q is given by
∀σ ∈ Eh ,
(ΠPnc
q) dσ =
q dσ.
One checks ( [2] p.110) that there exists C > 0 such that
|p−ΠPnc
p| ≤ C h ‖p‖1 ,
∣∣∣∇̃h(p−ΠPnc
∣∣∣ ≤ C ‖p‖1. (2.10)
Finally, we use the Raviart-Thomas spaces (see [3])
= {vh ∈ P
1 ; ∀σ ∈ EK , vh|K · nK,σ is a constant, and vh · n|∂Ω = 0} ,
RT0 = {vh ∈ RT
; ∀K ∈ Th, ∀σ ∈ EK , vh|Kσ · nKσ ,σ = vh|Lσ · nKσ ,σ}.
For all vh ∈ RT0, K ∈ Th and σ ∈ EK we set (vh ·nK,σ)σ = vh|K ·nK,σ. We define the operator ΠRT0 : H
RT0. For all v ∈ H
1, ΠRT0v ∈ RT0 is given by
∀K ∈ Th , ∀σ ∈ EK , (ΠRT0v · nK,σ)σ =
v dσ. (2.11)
2.3. The discrete operators
The equations (1.1)–(1.2) use the differential operators gradient, divergence and laplacian. Using the spaces
of section 2.2, we define their discrete counterparts. The discrete gradient ∇h : P
1 → P0 is defined by (2.8).
The discrete divergence operator divh : P0 → P
1 is built so that it is adjoint to the operator ∇h. We set for
all vh ∈ P0 and all triangle K ∈ Th
∀σ ∈ E inth , (divh vh)(xσ) =
3 |σ|
|Kσ|+ |Lσ|
(vLσ − vKσ ) · nK,σ ;
∀σ ∈ Eexth , (divh vh)(xσ) = −
3 |σ|
|Kσ|+ |Lσ|
vKσ · nK,σ. (2.12)
The first discrete laplacian ∆h : P
1 → P
1 ensures that the incompressibility constraint (1.2) is satisfied in a
discrete sense (see the proof of proposition 3.1 below). We set for all qh ∈ P
∆hqh = divh(∇hqh).
The second discrete laplacian ∆̃h : P0 → P0 is the usual operator in finite volume schemes [5]. We set for all
vh ∈ P0 and all triangle K ∈ Th
∆̃hvh|K =
σ∈EK∩E
τσ (vLσ − vKσ )−
σ∈EK∩E
τσ vKσ .
TITLE WILL BE SET BY THE PUBLISHER 5
In order to approximate the term (u · ∇)u in (1.1) we define a bilinear form b̃h : RT0 × P0 → P0 using the
well-known upwind scheme [5]. For all uh ∈ P0, vh ∈ P0, and all triangle K ∈ Th we set
b̃h(uh,vh)
σ∈EK∩E
(u · nK,σ)
σ vK + (u · nK,σ)
σ vLσ
. (2.13)
We have set a+ = max(a, 0), a− = min(a, 0) for all a ∈ R. Lastly, we define the trilinear form bh : RT0 ×P0 ×
P0 → R
2 as follows. For all uh ∈ RT0, vh ∈ P0, wh ∈ P0, we set
bh(uh,vh,wh) =
|K|wK · b̃h(uh,vh)
. (2.14)
3. The scheme
In order to deal with the incompressibility constraint (1.2) we use a projection method. This kind of method
has been introduced by Chorin [4] and Temam [13]. The basic idea is the following. The time interval [0, T ]
is split with a time step k: [0, T ] =
n=0[tn, tn+1] with N ∈ N
∗ and tn = n k for all n ∈ {0, . . . , N}. For all
m ∈ {2, . . . , N}, we compute (see equation (3.2) below) a first velocity field ũmh ≃ u(tm) using only equation
(1.1). We use a second-order BDF scheme for the discretization in time. We then project ũmh (see equation (3.4)
below) over a subspace of P0. We get a a pressure field p
h ≃ p(tm) and a second velocity field u
h ≃ u(tm),
which fulfills the incompressibility constraint (1.2) in a discrete sense. The algorithm goes as follows. For all
m ∈ {0, . . . , N}, we set fmh = ΠP0 f(tm). Since the operator ΠP0 is stable for the L
2-norm we get
|fmh | = |ΠP0 f(tm)| ≤ |f(tm)| ≤ ‖f‖C(0,T ;L2). (3.1)
We start with the initial values
u0h ∈ P0 ∩RT0 , u
h ∈ P0 ∩RT0 p
h ∈ P0 ∩ L
For all n ∈ {1, . . . , N}, (ũn+1h , p
h ) is deduced from (ũ
h , p
h) as follows.
• ũn+1h ∈ P0 is given by
3 ũn+1h − 4u
h + u
∆̃hũ
h + b̃h(2u
h − u
h , ũ
h ) +∇hp
h = f
h , (3.2)
• pn+1h ∈ P
1 ∩ L
0 is the solution of
h − p
divh ũ
h , (3.3)
• un+1h ∈ P0 is deduced by
un+1h = ũ
h − p
h). (3.4)
Existence and unicity of a solution to equation (3.2) is classical ( [5] for example). The convection term in (3.2)
is well defined thanks to the following result.
Proposition 3.1. For all m ∈ {0, . . . , N} we have umh ∈ RT0 .
Proof. If m ∈ {0, 1} the result holds by definition. If m ∈ {2, . . . , N} we apply the operator divh to (3.3) and
compare with (3.4). We get divh u
h = 0. Using definition (2.12) we get u
h ∈ RT0.
6 TITLE WILL BE SET BY THE PUBLISHER
Let us show that equation (3.3) also has a unique solution. Let qh ∈ P
1 ∩ L
0 such that ∆hqh = 0. According
to proposition 4.4 we have for all qh ∈ P0
−(∆hqh, qh) = −
divh(∇hqh), qh
= (∇hqh,∇hqh) = |∇hqh|
Therefore we have ∇hqh = 0, so that qh = 0 since qh ∈ L
0. We have thus proved the unicity of a solution for
(3.3). It is also the case for the associated linear system. It implies that this linear system has indeed a solution.
Hence it is also the case for equation (3.3). Note finally that since umh ∈ P0 ∩RT0, we have divu
h = 0 for all
m ∈ {0, . . . , N}. Hence the incompressibility condition (1.2) is fulfilled.
4. Properties of the discrete operators
We show that the differential operators in (1.1)–(1.2) and the operators defined in section 2.3 share similar
properties.
4.1. Properties of the discrete convective term
We define b̃ : H1 ×H1 → L2. For all u ∈ H1 and v = (v1, v2) ∈ H
1 we set b̃(u,v) =
div(v1 u), div(v2 u)
We show that the operator b̃h is a consistent approximation of b̃.
Proposition 4.1. There exists a constant C > 0 such that for all v ∈ H2 and all u ∈ H2 ∩ H10 satisfying
divu = 0
‖ΠP0 b̃(u,v)− b̃h(ΠRT0u, Π̃P0v)‖−1,h ≤ C h ‖u‖2 ‖v‖1.
Proof. We set uh = ΠRT0u and vh = Π̃P0v. Let K ∈ Th. According to the divergence formula and (2.6) we
ΠP0 b̃(u,v)|K =
σ∈EK∩E
v (u · n) dσ.
On the other hand, let us rewrite b̃h(uh,vh). Let σ ∈ EK ∩ E
h . Setting
vK,Lσ =
vK si (uh · nK,σ)σ ≥ 0
vLσ si (uh · nK,σ)σ < 0
one checks that vK (uh · nK,σ)
σ + vLσ (uh · nK,σ)
σ = vK,Lσ (uh · nK,σ)σ. Using (2.11), we deduce from (2.13)
b̃h(uh,vh)|K =
σ∈EK∩E
vK,Lσ (uh · nK,σ) dσ.
ΠP0 b̃(u,v)− b̃h(uh,vh)
σ∈EK∩E
(v − vK,Lσ) (uh · n) dσ.
Let ψh ∈ P0. We have
ΠP0 b̃(u,v) − b̃h(uh,vh),ψh
σ∈EK∩E
(v − vK,Lσ) (uh · n) dσ
σ∈Eint
(ψKσ −ψLσ)
(v − vKσ ,Lσ) (uh · n) dσ.
TITLE WILL BE SET BY THE PUBLISHER 7
Let σ ∈ E inth . We consider the quadrilateral Dσ defined by xKσ , xLσ and the vertex of σ. We set
DK,Lσ =
Dσ ∩K si (uh · nK,σ)σ ≥ 0
Dσ ∩ Lσ si (uh · nK,σ)σ < 0
Using a Taylor expansion and a density argument (see [14]) one checks that
|v − vKσ ,Lσ | dσ ≤ C h
DKσ,Lσ
|∇v (y)|2 dy
ΠP0 b̃(u,v) − b̃h(ΠRT0u, Π̃P0v),ψh
≤ C h ‖u‖H2
σ∈Eint
|ψLσ −ψKσ |
σ∈Eint
DKσ,Lσ
|∇v (y)|2 dy
so that
ΠP0b̃(u,v)− b̃h(ΠRT0u, Π̃P0v),ψh
)∣∣∣ ≤ C h ‖u‖H2 ‖ψh‖1,h ‖v‖1. Using then definition (2.5), we get
the result.
Let v ∈ L∞ ∩H1 and u ∈ H1 with divu ≥ 0 a.e. in Ω. Integrating by parts one checks that
v · b̃(u,v) dx =
divu dx ≥ 0. The operator bh shares a similar property.
Proposition 4.2. Let uh ∈ RT0 such that divuh ≥ 0. For all vh ∈ P0 we have
bh(uh,vh,vh) ≥ 0.
Proof. Remember that for all edges σ ∈ E inth , two triangles Kσ et Lσ share σ as an edge. We denote by Kσ
the one such that uσ · nKσ ,σ ≥ 0. Using the algebraic identity 2 a (a− b) = a
2 − b2 + (a − b)2 we deduce from
(2.14)
2 bh(uh,vh,vh) = 2
σ∈Eint
|σ|vKσ · (vKσ − vLσ) (uh · nKσ,σ)
σ∈Eint
|vKσ|
2 − |vLσ |
2 + |vKσ − vLσ |
(uh · nKσ,σ)
so that 2 bh(uh,vh,vh) ≥
σ∈Eint
|vKσ|
2 − |vLσ |
(uh · nKσ,σ). This sum can be written as a sum over
the triangles of the mesh. We get
2 bh(uh,vh,vh) ≥
|vKσ |
σ∈EK∩E
|σ| (uh · nKσ,σ).
Using finally the divergence formula we get
2 bh(uh,vh,vh) ≥
|K| |vK |
divuh dx ≥ 0.
The following result states that the operator bh is stable for suitable norms.
8 TITLE WILL BE SET BY THE PUBLISHER
Proposition 4.3. There exists a constant C > 0 such that for all vh ∈ P0, wh ∈ P0, uh ∈ P0 satisfying
divuh = 0
|bh(uh,vh,vh)| ≤ C |uh| ‖vh‖h ‖vh‖h.
Proof. For all triangle K ∈ Th and all edge σ ∈ EK ∩ E
h , we have
(uh · nK,σ)
σ vK + (uh · nK,σ)
σ vLσ = (uh · nK,σ)σ vK − |(uh · nK,σ)σ| (vLσ − vK).
Using this splitting, we deduce from (2.14) bh(uh,vh,wh) = S1 + S2 with
vK ·wK
σ∈EK∩E
|σ| (uh · nK,σ)σ ,
S2 = −
σ∈EK∩E
|σ| |(uh · nK,σ)σ| (vLσ − vK).
By writing the sum over the edges as a sum over the triangles we have
S2 = −
σ∈Eint
|σ| |(uh · nK,σ)σ| (vLσ − vK) · (wLσ −wK).
Using the Cauchy-Schwarz inequality we get
|S2| ≤ h ‖uh‖∞
σ∈Eint
|vLσ − vKσ |
1/2 
σ∈Eint
|wLσ −wKσ |
Since uh ∈ RT0 we have [5] the inverse inequality h ‖uh‖∞ ≤ C |uh|. Using (2.1) and (2.2) we get
σ∈Eint
|vLσ − vKσ |
2 ≤ C
σ∈Eint
τσ |vLσ − vKσ |
2 ≤ C ‖vh‖
and in a similar way
σ∈Eint
|wLσ −wKσ |
2 ≤ C ‖wh‖
h. Thus |S2| ≤ C |uh| ‖vh‖h ‖wh‖h. On the other hand,
according to the divergence formula
|K| (vK ·wK)
divuh dx = 0.
By gathering the estimates for S1 and S2 we get the result.
4.2. Properties of the discrete divergence
The operators gradient and divergence are adjoint: if q ∈ H1 , v ∈ H1 with v · n|∂Ω = 0, we get (v,∇q) =
−(q, divv) by integrating by parts. For ∇h and divh we state the following.
Proposition 4.4. For all vh ∈ P0 and qh ∈ P
1 we have: (vh,∇hqh) = −(qh, divh vh).
Proof. According to (2.8)
(vh,∇hqh) =
|K|vK · ∇hqh|K =
|σ|qh(xσ)nK,σ
TITLE WILL BE SET BY THE PUBLISHER 9
By writing this sum as a sum over the edges we get
(vh,∇hqh) = −
σ∈Eint
|σ| qh(xσ) (vLσ − vKσ) · nKσ,σ +
σ∈Eext
|σ| qh(xσ)vKσ · nKσ,σ. (4.1)
On the other hand, using a quadrature formula
−(qh, divh vh) = −
qh(xσ) (divh vh)(xσ).
By writing this sum as a sum over the edges of the mesh we get
−(qh, divh vh) = −
σ∈Eint
( |Kσ|
qh(xσ) (divhvh)(xσ)−
σ∈Eext
qh(xσ) (divh vh)(xσ).
Using definition (2.12) and comparing with (4.1) we get the result.
The divergence operator and the spaces L20, H
0 satisfy the following property, called inf-sup (or Babuška-Brezzi)
condition (see [9] for example). There exists a constant C > 0 such that
(q, divv)
‖v‖1|q|
≥ C. (4.2)
We will now show that the operator divh and the spaces P0 ∩ L
0, P0 satisfy an analogous property. The proof
uses the following lemma.
Lemma 4.1. There exists a constant C > 0 such that
∀ qh ∈ P
1 ∩ L
0 , sup
vh∈P0\{0}
(qh, divh vh)
‖vh‖h
≥ C h ‖qh‖1,h.
Proof. If qh = 0 the result is trivial. Let qh ∈ P
0\{0}. Let vh = ∇hqh ∈ P0\{0}. Using proposition 4.4
we have
−(qh, divhvh) = (vh,∇hqh) = |∇hqh|
2 = |∇hqh| |vh|.
Using (2.3) and (2.4) we get −(qh, divhvh) ≥ C h ‖qh‖1,h ‖vh‖h.
We now state the result.
Proposition 4.5. There exists a constant C > 0 such that for all qh ∈ P
1 ∩ L
vh∈P0\{0}
(qh, divh vh)
‖vh‖h
≥ C |qh|.
Proof. If qh = 0 the result is trivial. Let qh ∈ P
0\{0}. According to (4.2) there exists v ∈ H
0 such that
divv = −qh and ‖v‖1 ≤ C |qh|. (4.3)
We set vh = ΠPc
v. We want to estimate −
qh, divh(ΠP0vh)
. Since ∇hqh ∈ P0 we deduce from proposition
qh, divh(ΠP0vh)
= (ΠP0vh,∇hqh) = (vh,∇hqh).
By splitting the last term we get
qh, divh(ΠP0vh)
= (v,∇hqh)− (v − vh,∇hqh). (4.4)
10 TITLE WILL BE SET BY THE PUBLISHER
We bound the right-hand side of (4.4). Using (2.7) and (4.3) we have
|v − vh| = |v −ΠPc
v| ≤ C h ‖v‖1 ≤ C h |qh|.
Thus, using the Cauchy-Schwarz inequality, we get
|(v − vh,∇hqh)| ≤ C h |qh| |∇hqh| ≤ C h |qh| ‖qh‖1,h.
We estimate the other term as follows. Integrating by parts we get
(v,∇hqh) = −(qh, divv) +
qh (v · nK,σ) dσ.
We have −(qh, divv) = |qh|
2 thanks to (4.3). On the other hand
qh (v · nK,σ) dσ =
σ∈Eint
qh (v · nKσ,σ) dσ
since v|∂Ω = 0. Using [2] p.269 and (4.3) we have
∣∣∣∣∣
qh (v · nK,σ) dσ
∣∣∣∣∣
≤ C h ‖v‖1 ‖qh‖1,h ≤ C h |qh| ‖qh‖1,h.
Hence we get (v,∇hqh) ≥ (|qh| − C h ‖qh‖1,h) |qh|. Thus we deduce from (4.4)
qh, divh(ΠP0vh)
≥ (|qh| − C h ‖qh‖1,h) |qh|. (4.5)
We now introduce the norm ‖.‖h. We have vh = ΠPc
v ∈ Pc1 ⊂ H
1. From [5] p. 776 we deduce ‖ΠP0vh‖h ≤
C ‖vh‖1. Since ΠPc
is stable for the H1 norm, using (4.3), we get
‖vh‖1 = ‖ΠPc
v‖1 ≤ ‖v‖1 ≤ C |qh|.
Therefore ‖ΠP0vh‖h ≤ C |qh|. Using this inequality in (4.5) we obtain that there exists C1 > 0 and C2 > 0
such that
qh, divh(ΠP0vh)
≥ (C1 |qh| − C2 h ‖qh‖1,h) ‖ΠP0vh‖h.
We deduce from this
vh∈P0\{0}
(qh, divh vh)
‖vh‖h
≥ C1 |qh| − C2 h ‖qh‖1,h.
Let us combine this result with lemma 4.1. Since
∀ t ≥ 0 , max
C t , C1 |qh| − C2 t
C + C2
|qh| ,
we finally get the result.
4.3. Properties of the discrete laplacian
We recall from [14] the coercivity of the laplacian operator.
Proposition 4.6. For all uh ∈ P0 and vh ∈ P0 we have
−(∆̃huh,uh) = ‖uh‖
h , −(∆̃huh,vh) ≤ ‖uh‖h ‖vh‖h.
TITLE WILL BE SET BY THE PUBLISHER 11
5. Stability of the scheme
We first prove an estimate for the computed velocity (theorem 5.1). We show a similar result for the
increments in time (lemma 5.2). Using the inf-sup condition (proposition 4.5), we infer from it some estimates
on the pressure (theorem 5.2).
Lemma 5.1. For all m ∈ {0, . . . , N} et n ∈ {0, . . . , N} we have
(umh ,∇hp
h) = 0 , |u
2 − |ũmh |
2 + |umh − ũ
2 = 0.
Proof. First, using propositions 3.1 and 4.4, we get (umh ,∇hp
h) = −(p
h, divhu
h ) = 0. Also, we deduce from
(3.4)
2 (umh ,u
h − ũ
h ) = −
umh ,∇h(p
h − p
Using the algebraic identity 2 a (a− b) = a2 − b2 + (a− b)2 we get
2 (umh ,u
h − ũ
h ) = |u
2 − |ũmh |
2 + |umh − ũ
2 = 0.
We introduce the following hypothesis on the initial data.
(H1) There exists C > 0 such that |u0h|+ |u
h|+ k|∇hp
h| ≤ C.
Hypothesis (H1) is fulfilled if we set u0h = ΠRT0u0 and we use a semi-implicit Euler scheme to compute u
We have the following stability result.
Theorem 5.1. We assume that the initial values of the scheme fulfill (H1). For all m ∈ {2, . . . , N} we have
|umh |
2 + k
‖ũnh‖
h ≤ C. (5.1)
Proof. Let m ∈ {2, . . . , N} and n ∈ {1, . . . ,m− 1}. Taking the scalar product of (3.2) with 4 k ũn+1h we get
3 ũn+1h − 4u
h + u
, 4 k ũn+1h
(∆̃hũ
h , ũ
+4 k bh(2u
h − u
h , ũ
h , ũ
h ) + 4 k (∇hp
h, ũ
h ) = 4 k (f
h , ũ
h ). (5.2)
First of all, using lemma 5.1 and proceeding as in [10], we get
ũn+1h ,
3 ũn+1h − 4u
h + u
= |un+1h |
2 − |unh|
2 + |2un+1h − u
2 − |2unh − u
+ |un+1h − 2u
h + u
2 + 6 |ũn+1h − u
According to proposition 4.6 we have − 4 k
(∆̃hũ
h , ũ
h ) =
‖ũn+1h ‖
h. Also, according to lemma 5.1 and
(3.4)
4 k (∇hp
h, ũ
h ) = 4 k (∇hp
h, ũ
h − u
(|∇pn+1h |
2 − |∇pnh|
2 − |∇pn+1h −∇p
12 TITLE WILL BE SET BY THE PUBLISHER
Multiplying equation (3.4) by 4 k∇h(p
h − p
h) and using the Young inequality we get
|∇(pn+1h − p
2 ≤ 3 |un+1h − ũ
According to proposition 4.2, we have 4 k bh(2u
h − u
h , ũ
h , ũ
h ) ≥ 0. At last using the Cauchy-Schwarz
inequality, (2.3) and (3.1) we have
4 k (fn+1h , ũ
h ) ≤ 4 k |f
h | |ũ
h | ≤ C k ‖f‖C(0,T ;L2) ‖ũ
h ‖h.
Using the Young inequality we get
4 k (fn+1h , ũ
h ) ≤ 3 k ‖ũ
h + C k ‖f‖
C(0,T ;L2).
Thus we deduce from (5.2)
|un+1h |
2 − |unh |
2 + |2un+1h − u
2 − |2unh − u
2 + |un+1h − 2u
h + u
+3 |ũn+1h − u
2 + k ‖ũn+1h ‖
(|∇hp
2 − |∇hp
2) ≤ C k.
Summing from n = 1 to m− 1 we have
|umh |
2 + |2umh − u
2 + 3
|ũn+1h − u
2 + k
‖ũn+1h ‖
≤ C + 4 |u1h|
2 + |2u1h − u
2 + k2 |∇hp
Using hypothesis (H1) we get (5.1).
We now want to estimate the computed pressure. From now on, we make the following hypothesis on the data
f ∈ C(0, T ;L2) , ft ∈ L
2(0, T ;L2) , u0 ∈ H
2 ∩H10 , divu0 = 0.
One shows that if the data u0 and f fulfill a compatibility condition [11] there exists a solution (u, p) to the
equations (1.1)–(1.2) such that
u ∈ C(0, T ;H2) , ut ∈ C(0, T ;L
2) , ∇p ∈ C(0, T ;L2).
We introduce the following hypothesis on the initial values of the scheme: there exists a constant C > 0 such
(H2) |u0h − u0|+
‖u1h − u(t1)‖∞ + |p
h − p(t1)| ≤ C h , |u
h − u
h| ≤ C k.
One checks easily that this hypothesis implies (H1). We have the following result.
Lemma 5.2. We assume that the initial values of the scheme fulfill (H2). Then there exists a constant C > 0
such that for all m ∈ {1, . . . , N}
|umh − u
h | ≤ C.
Proof. Using proposition 4.1 one proceeds as in [14]. The difference lies in the way we bound the term ∇hp
We use the splitting
p1h = (p
h −ΠPnc1 p(t1)) + (ΠP
p(t1)− p(t1)) + p(t1).
TITLE WILL BE SET BY THE PUBLISHER 13
Using an inverse inequality [2] we have
p1h −ΠPnc1 p(t1)
)∣∣ ≤
∣∣p1h −ΠPnc1 p(t1)
(∣∣p1h − p(t1)
∣∣p(t1)−ΠPnc
p(t1)
∣∣) .
Using (2.10) and hypothesis (H2) we get
p1h −ΠPnc1 p(t1)
)∣∣ ≤ C ‖p(t1)‖1 ≤ C ‖p‖C(0,T ;H1).
According to (2.10) we also have
∣∣∇h(p(t1)−ΠPnc
p(t1))
∣∣ ≤ C ‖p(t1)‖1 ≤ C ‖p‖C(0,T ;H1). Lastly |∇p(t1)| ≤
‖p‖C(0,T ;H1). Thus we get |∇hp
h| ≤ C.
Theorem 5.2. We assume that the initial values of the scheme fulfull (H2). There exists a constant C > 0
such that for all m ∈ {2, . . . , N}
|pnh|
2 ≤ C.
Proof. Let m ∈ {2, . . . , N}. We set n = m− 1. Using the inf-sup condition (4.5) and proposition 4.4, we get
that there exists vh ∈ P0\{0} such that
C ‖vh‖h |p
h | ≤ −(p
h , divh vh) = (∇hp
h ,vh). (5.3)
Plugging (3.4) into (3.2) we have
h = −
3un+1h − 4u
h + u
∆̃hũ
h − b̃h(2u
h − u
h , ũ
h ) + f
so that
h ,vh) = −
3un+1h − 4u
h + u
∆̃hũ
h ,vh
− bh(2u
h − u
h , ũ
h ,vh) + (f
h ,vh).
Thanks to proposition 4.3 and theorem 5.1 we have
∣∣bh(2unh − u
h , ũ
h ,vh)
2 |unh|+ |u
‖ũn+1h ‖h ‖vh‖h ≤ C ‖ũ
h ‖h ‖vh‖h.
According to proposition 4.6 we have
∆̃hũ
h ,vh
≤ ‖ũn+1h ‖h ‖vh‖h. Using the Cauchy-Schwarz inequality,
(2.3) and (3.1) we have
(fn+1h ,vh) ≤ |f
h | |vh| ≤ C |vh| ≤ C ‖vh‖h
and in a similar way
3un+1h − 4u
h + u
)∣∣∣∣ ≤ C
3un+1h − 4u
h + u
∣∣∣∣ ‖vh‖h.
Thus we get
h ,vh) ≤ C + C
|3un+1h − 4u
h + u
+ ‖ũn+1h ‖h
‖vh‖h.
By comparing with (5.3) we get
|pn+1h | ≤ C + C
|3un+1h − 4u
h + u
+ ‖ũn+1h ‖h
14 TITLE WILL BE SET BY THE PUBLISHER
Squaring and summing from n = 1 to m− 1 we obtain
|pnh|
2 ≤ C + C k
|3un+1h − 4u
h + u
+ C k
‖ũn+1h ‖
The last term on the right-hand side is bounded, thanks to theorem 5.1. And since
3un+1h − 4u
h + u
h = 3(u
h − u
h)− (u
h − u
h ) = 3 δu
h − δu
we deduce from lemma 5.2
|3un+1h − 4u
h + u
≤ C k
|δunh |
References
[1] S. Boivin, F. Cayre, J. M Herard, A finite volume method to solve the Navier-Stokes equations for incompressible flows on
unstructured meshes, Int. J. Therm. Sci. 39 (2000) 806–825.
[2] S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods, Springer, 2002.
[3] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, 1991.
[4] J. Chorin, On the convergence of discrete approximations to the Navier-Stokes equations, Math. Comp. 23 (1969) 341–353.
[5] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, P.G. Ciarlet and J.L. Lions eds, North-Holland, 2000.
[6] R. Eymard and R. Herbin, A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space
dimensions, Int.J. Finite Volumes (2005).
[7] R. Eymard, J. C. Latché and R. Herbin, Convergence analysis of a colocated finite volume scheme for the incompressible
Navier-Stokes equations on general 2 or 3D meshes, SIAM J. Numer. Anal. 45(1) (2007) 1–36.
[8] S. Faure, Stability of a colocated finite volume scheme for the Navier-Stokes equations, Num. Meth. PDE 21(2) (2005) 242–271.
[9] V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes equations: Theory and Algorithms, Springer, 1986.
[10] J. L. Guermond, Some implementations of projection methods for Navier-Stokes equations, M2AN 30(5) (1996) 637–667.
[11] J. G. Heywood and R. Rannacher, Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of
solutions and second-order error estimates for spatial discretization, SIAM J. Numer. Anal. 19(26) (1982) 275–311.
[12] D. Kim and H. Choi, A second-order time-accurate finite volume method for unsteady incompressible flow on hybrid unstruc-
tured grids, J. Comp. Phys. 162 (2000) 411–428.
[13] R. Temam, Sur l’approximation de la solution des équations de Navier-Stokes par la méthode de pas fractionnaires II, Arch.
Rat. Mech. Anal. 33 (1969) 377–385.
[14] S. Zimmermann, Stability of a colocated finite volume for the incompressible Navier-Stokes equations, arXiv :0704.0772 (2006).
[15] S. Zimmermann, Étude et implémentation de méthodes de volumes finis pour les fluides incompressibles, PhD, Blaise Pascal
University, France (2006).

  We introduce a finite volume scheme for the two-dimensional incompressible
Navier-Stokes equations. We use a triangular mesh. The unknowns for the
velocity and pressure are respectively piecewise constant and affine. We use a
projection method to deal with the incompressibility constraint. We show that
the differential operators in the Navier-Stokes equations and their discrete
counterparts share similar properties. In particular we state an inf-sup
(Babuska-Brezzi) condition. Using these properties we infer the stability of
the scheme.

Introduction
We consider the flow of an incompressible fluid in a polyhedral set Ω ⊂ R2 during the time interval [0, T ].
The velocity field u : Ω× [0, T ] → R2 and the pressure field p : Ω× [0, T ] → R satisfy the Navier-Stokes equations
∆u+ (u · ∇)u+∇p = f , (1.1)
div u = 0 , (1.2)
with the boundary and initial conditions
u|∂Ω = 0 , u|t=0 = u0.
The terms ∆u and (u·∇)u are associated with the physical phenomena of diffusion and convection, respectively.
The Reynolds number Re measures the influence of convection in the flow. For equations (1.1)–(1.2), finite
element and finite difference methods are well known and mathematical studies are available (see [9] for example).
Keywords and phrases: Incompressible fluids, Navier-Stokes equations, projection methods, finite volume.
1 17 rue Barrème - 69006 LYON. e-mail: Sebastien.Zimmermann@ec-lyon.fr
c© EDP Sciences, SMAI 1999
http://arxiv.org/abs/0704.0783v2
2 TITLE WILL BE SET BY THE PUBLISHER
For finite volume schemes, numerous computations have been conducted ( [12] and [1] for example). However,
few mathematical results are available in this case. Let us cite Eymard and Herbin [6] and Eymard,
Latché and Herbin [7]. In order to deal with the incompressibility constraint (1.2), these works use a
penalization method. Another way is to use the projection methods which have been introduced by Chorin [4]
and Temam [13]. This is the case in Faure [8] where the mesh is made of squares. In Zimmermann [14] the
mesh is made of triangles, so that more complex geometries can be considered. In the present paper the mesh
is also made of triangles, but we consider a different discretization for the pressure. It leads to a linear system
with a better-conditioned matrix. The layout of the article is the following. We first introduce in section 2
the discrete setting. We state (section 2.1) some notations and hypotheses on the mesh. We define (section
2.2) the spaces we use to approximate the velocity and pressure. We define also (section 2.3) the operators
we use to approximate the differential operators in (1.1)–(1.2). Combining this with a projection method, we
build the scheme in section 3. In order to provide a mathematical analysis, we show in section 4 that the
differential operators in (1.1)–(1.2) and their discrete counterparts share similar properties. In particular, the
discrete operators for the gradient and the divergence are adjoint. The discrete operator for the convection term
is positive, stable and consistent. The discrete operator for the divergence satisfy an inf-sup (Babuška-Brezzi)
condition. From these properties we deduce in section 5 the stability of the scheme.
We conclude with some notations. The spaces (L2, |.|) and (L∞, ‖.‖∞) are the usual Lebesgue spaces and we
set L20 = {q ∈ L
q(x) dx = 0}. Their vectorial counterparts are (L2, |.|) and (L∞, ‖.‖∞) with L
2 = (L2)2
and L∞ = (L∞). For k ∈ N∗, (Hk, ‖ ·‖k) is the usual Sobolev space. Its vectorial counterpart is (H
k, ‖.‖k) with
Hk = (Hk)2. For k = 1, the functions of H1 with a null trace on the boundary form the space H10. Also, we set
∇u = (∇u1,∇u2)
T if u = (u1, u2) ∈ H
1. If X ⊂ L2 is a Banach space, we define C(0, T ;X) (resp. L2(0, T ;X))
as the set of the applications g : [0, T ] → X such that t → |g(t)| is continuous (resp. square integrable). The
norm ‖.‖C(0,T ;X) is defined by ‖g‖C(0,T ;X) = sups∈[0,T ] |g(s)|. In all calculations, C is a generic positive constant,
depending only on Ω, u0 and f .
2. Discrete setting
First, we introduce the spaces and the operators needed to build the scheme.
2.1. The mesh
Let Th be a triangular mesh of Ω. The circumscribed circle of a triangle K ∈ Th is centered at xK and has
the diameter hK . We set h = maxK∈Th hK . We assume that all the interior angles of the triangles of the mesh
are less than π
, so that xK ∈ K. The set of the edges of the triangle K ∈ Th is EK . The symbol nK,σ denotes
the unit vector normal to an edge σ ∈ EK and pointing outward K. We denote by Eh the set of the edges of the
mesh. We distinguish the subset E inth ⊂ Eh (resp. E
h ) of the edges located inside Ω (resp. on ∂Ω). The middle
of an edge σ ∈ Eh is xσ and its length |σ|. For each edge σ ∈ E
h , let Kσ and Lσ be the two triangles having
σ in common. We set dσ = d(xKσ ,xLσ ). For all σ ∈ E
h , only the triangle Kσ located inside Ω is defined and
we set dσ = d(xKσ ,xσ). Then for all σ ∈ Eh we set τσ =
. As in [5] we assume the following on the mesh:
there exists C > 0 such that
∀σ ∈ Eh , dσ ≥ C |σ| and |σ| ≥ C h.
It implies that there exists C > 0 such that
∀σ ∈ E inth , τσ = |σ|/dσ ≥ C. (2.1)
2.2. The discrete spaces
We first define
P0 = {q ∈ L
2 ; ∀K ∈ Th, q|K is a constant} , P0 = (P0)
TITLE WILL BE SET BY THE PUBLISHER 3
For the sake of concision, we set for all qh ∈ P0 (resp. vh ∈ P0) and all triangle K ∈ Th: qK = qh|K (resp.
vK = vh|K). Although P0 6⊂ H
1, we define the discrete equivalent of a H1 norm as follows. For all vh ∈ P0
we set
‖vh‖h =
σ∈Eint
τσ |vLσ − vKσ |
σ∈Eext
τσ |vKσ |
. (2.2)
We have [5] a Poincaré-like inequality: there exists C > 0 such that for all vh ∈ P0
|vh| ≤ C ‖vh‖h. (2.3)
We also have [14] an inverse inequality: there exists C > 0 such that for all vh ∈ P0
h ‖vh‖h ≤ C |vh|. (2.4)
From the norm ‖.‖h we deduce a dual norm. For all vh ∈ P0 we set
‖vh‖−1,h = sup
(vh,ψh)
‖ψh‖h
. (2.5)
For all uh ∈ P0 and vh ∈ P0 we have (uh,vh) ≤ ‖uh‖−1,h ‖vh‖h. We define the projection operator ΠP0 :
L2 → P0 as follows. For all w ∈ L
2, ΠP0w ∈ P0 is given by
∀K ∈ Th , (ΠP0w)|K =
w(x) dx. (2.6)
We easily check that for all w ∈ L2 and vh ∈ P0 we have (ΠP0w,vh) = (w,vh). We deduce from this that ΠP0
is stable for the L2 norm. We define also the operator Π̃P0 : H
2 → P0. For all w ∈ H
2, Π̃P0w ∈ P0 is given by
∀K ∈ Th , Π̃P0w|K = w(xK).
According to the Sobolev embedding theorem, w ∈ H2 is a.e. equal to a continuous function. Therefore the
definition above makes sense. We introduce also the finite element spaces
P d1 = {v ∈ L
2 ; ∀K ∈ Th, v|K is affine} ,
Pnc1 = {vh ∈ P
1 ; ∀σ ∈ E
h , vh|Kσ (xσ) = vh|Lσ(xσ) ,
Pc1 = {vh ∈ (P
2 ; vh is continuous and vh|∂Ω = 0}.
We have Pc1 ⊂ H
0. We define ΠPc1 : H
0 → P
1. For all v = (v1, v2) ∈ H
0, ΠPc1v = (v
h) ∈ P
1 is given by
∀φh = (φ
h) ∈ P
∇vih,∇φ
∇vi,∇φ
The operator ΠPc
is stable for the H1 norm. One checks ( [2] p. 110) that there exists C > 0 such that for all
v ∈ H1
|v −ΠPc
v| ≤ C h ‖v‖1. (2.7)
Let us address now the space Pnc1 . If qh ∈ P
1 , we have usually ∇qh 6∈ L
2. Thus we define the operator
∇h : P
1 → P0 by setting for all qh ∈ P0 and all triangle K ∈ Th
∇hqh|K =
∇qh dx. (2.8)
4 TITLE WILL BE SET BY THE PUBLISHER
The associated norm is defined by
‖qh‖1,h =
2 + |∇hqh|
We have a Poincaré-like inequality : there exists C > 0 such that for all qh ∈ P
1 ∩ L
|qh| ≤ C |∇hqh|. (2.9)
We define the projection operator ΠPnc
. For all q ∈ H1, ΠPnc
q is given by
∀σ ∈ Eh ,
(ΠPnc
q) dσ =
q dσ.
One checks ( [2] p.110) that there exists C > 0 such that
|p−ΠPnc
p| ≤ C h ‖p‖1 ,
∣∣∣∇̃h(p−ΠPnc
∣∣∣ ≤ C ‖p‖1. (2.10)
Finally, we use the Raviart-Thomas spaces (see [3])
= {vh ∈ P
1 ; ∀σ ∈ EK , vh|K · nK,σ is a constant, and vh · n|∂Ω = 0} ,
RT0 = {vh ∈ RT
; ∀K ∈ Th, ∀σ ∈ EK , vh|Kσ · nKσ ,σ = vh|Lσ · nKσ ,σ}.
For all vh ∈ RT0, K ∈ Th and σ ∈ EK we set (vh ·nK,σ)σ = vh|K ·nK,σ. We define the operator ΠRT0 : H
RT0. For all v ∈ H
1, ΠRT0v ∈ RT0 is given by
∀K ∈ Th , ∀σ ∈ EK , (ΠRT0v · nK,σ)σ =
v dσ. (2.11)
2.3. The discrete operators
The equations (1.1)–(1.2) use the differential operators gradient, divergence and laplacian. Using the spaces
of section 2.2, we define their discrete counterparts. The discrete gradient ∇h : P
1 → P0 is defined by (2.8).
The discrete divergence operator divh : P0 → P
1 is built so that it is adjoint to the operator ∇h. We set for
all vh ∈ P0 and all triangle K ∈ Th
∀σ ∈ E inth , (divh vh)(xσ) =
3 |σ|
|Kσ|+ |Lσ|
(vLσ − vKσ ) · nK,σ ;
∀σ ∈ Eexth , (divh vh)(xσ) = −
3 |σ|
|Kσ|+ |Lσ|
vKσ · nK,σ. (2.12)
The first discrete laplacian ∆h : P
1 → P
1 ensures that the incompressibility constraint (1.2) is satisfied in a
discrete sense (see the proof of proposition 3.1 below). We set for all qh ∈ P
∆hqh = divh(∇hqh).
The second discrete laplacian ∆̃h : P0 → P0 is the usual operator in finite volume schemes [5]. We set for all
vh ∈ P0 and all triangle K ∈ Th
∆̃hvh|K =
σ∈EK∩E
τσ (vLσ − vKσ )−
σ∈EK∩E
τσ vKσ .
TITLE WILL BE SET BY THE PUBLISHER 5
In order to approximate the term (u · ∇)u in (1.1) we define a bilinear form b̃h : RT0 × P0 → P0 using the
well-known upwind scheme [5]. For all uh ∈ P0, vh ∈ P0, and all triangle K ∈ Th we set
b̃h(uh,vh)
σ∈EK∩E
(u · nK,σ)
σ vK + (u · nK,σ)
σ vLσ
. (2.13)
We have set a+ = max(a, 0), a− = min(a, 0) for all a ∈ R. Lastly, we define the trilinear form bh : RT0 ×P0 ×
P0 → R
2 as follows. For all uh ∈ RT0, vh ∈ P0, wh ∈ P0, we set
bh(uh,vh,wh) =
|K|wK · b̃h(uh,vh)
. (2.14)
3. The scheme
In order to deal with the incompressibility constraint (1.2) we use a projection method. This kind of method
has been introduced by Chorin [4] and Temam [13]. The basic idea is the following. The time interval [0, T ]
is split with a time step k: [0, T ] =
n=0[tn, tn+1] with N ∈ N
∗ and tn = n k for all n ∈ {0, . . . , N}. For all
m ∈ {2, . . . , N}, we compute (see equation (3.2) below) a first velocity field ũmh ≃ u(tm) using only equation
(1.1). We use a second-order BDF scheme for the discretization in time. We then project ũmh (see equation (3.4)
below) over a subspace of P0. We get a a pressure field p
h ≃ p(tm) and a second velocity field u
h ≃ u(tm),
which fulfills the incompressibility constraint (1.2) in a discrete sense. The algorithm goes as follows. For all
m ∈ {0, . . . , N}, we set fmh = ΠP0 f(tm). Since the operator ΠP0 is stable for the L
2-norm we get
|fmh | = |ΠP0 f(tm)| ≤ |f(tm)| ≤ ‖f‖C(0,T ;L2). (3.1)
We start with the initial values
u0h ∈ P0 ∩RT0 , u
h ∈ P0 ∩RT0 p
h ∈ P0 ∩ L
For all n ∈ {1, . . . , N}, (ũn+1h , p
h ) is deduced from (ũ
h , p
h) as follows.
• ũn+1h ∈ P0 is given by
3 ũn+1h − 4u
h + u
∆̃hũ
h + b̃h(2u
h − u
h , ũ
h ) +∇hp
h = f
h , (3.2)
• pn+1h ∈ P
1 ∩ L
0 is the solution of
h − p
divh ũ
h , (3.3)
• un+1h ∈ P0 is deduced by
un+1h = ũ
h − p
h). (3.4)
Existence and unicity of a solution to equation (3.2) is classical ( [5] for example). The convection term in (3.2)
is well defined thanks to the following result.
Proposition 3.1. For all m ∈ {0, . . . , N} we have umh ∈ RT0 .
Proof. If m ∈ {0, 1} the result holds by definition. If m ∈ {2, . . . , N} we apply the operator divh to (3.3) and
compare with (3.4). We get divh u
h = 0. Using definition (2.12) we get u
h ∈ RT0.
6 TITLE WILL BE SET BY THE PUBLISHER
Let us show that equation (3.3) also has a unique solution. Let qh ∈ P
1 ∩ L
0 such that ∆hqh = 0. According
to proposition 4.4 we have for all qh ∈ P0
−(∆hqh, qh) = −
divh(∇hqh), qh
= (∇hqh,∇hqh) = |∇hqh|
Therefore we have ∇hqh = 0, so that qh = 0 since qh ∈ L
0. We have thus proved the unicity of a solution for
(3.3). It is also the case for the associated linear system. It implies that this linear system has indeed a solution.
Hence it is also the case for equation (3.3). Note finally that since umh ∈ P0 ∩RT0, we have divu
h = 0 for all
m ∈ {0, . . . , N}. Hence the incompressibility condition (1.2) is fulfilled.
4. Properties of the discrete operators
We show that the differential operators in (1.1)–(1.2) and the operators defined in section 2.3 share similar
properties.
4.1. Properties of the discrete convective term
We define b̃ : H1 ×H1 → L2. For all u ∈ H1 and v = (v1, v2) ∈ H
1 we set b̃(u,v) =
div(v1 u), div(v2 u)
We show that the operator b̃h is a consistent approximation of b̃.
Proposition 4.1. There exists a constant C > 0 such that for all v ∈ H2 and all u ∈ H2 ∩ H10 satisfying
divu = 0
‖ΠP0 b̃(u,v)− b̃h(ΠRT0u, Π̃P0v)‖−1,h ≤ C h ‖u‖2 ‖v‖1.
Proof. We set uh = ΠRT0u and vh = Π̃P0v. Let K ∈ Th. According to the divergence formula and (2.6) we
ΠP0 b̃(u,v)|K =
σ∈EK∩E
v (u · n) dσ.
On the other hand, let us rewrite b̃h(uh,vh). Let σ ∈ EK ∩ E
h . Setting
vK,Lσ =
vK si (uh · nK,σ)σ ≥ 0
vLσ si (uh · nK,σ)σ < 0
one checks that vK (uh · nK,σ)
σ + vLσ (uh · nK,σ)
σ = vK,Lσ (uh · nK,σ)σ. Using (2.11), we deduce from (2.13)
b̃h(uh,vh)|K =
σ∈EK∩E
vK,Lσ (uh · nK,σ) dσ.
ΠP0 b̃(u,v)− b̃h(uh,vh)
σ∈EK∩E
(v − vK,Lσ) (uh · n) dσ.
Let ψh ∈ P0. We have
ΠP0 b̃(u,v) − b̃h(uh,vh),ψh
σ∈EK∩E
(v − vK,Lσ) (uh · n) dσ
σ∈Eint
(ψKσ −ψLσ)
(v − vKσ ,Lσ) (uh · n) dσ.
TITLE WILL BE SET BY THE PUBLISHER 7
Let σ ∈ E inth . We consider the quadrilateral Dσ defined by xKσ , xLσ and the vertex of σ. We set
DK,Lσ =
Dσ ∩K si (uh · nK,σ)σ ≥ 0
Dσ ∩ Lσ si (uh · nK,σ)σ < 0
Using a Taylor expansion and a density argument (see [14]) one checks that
|v − vKσ ,Lσ | dσ ≤ C h
DKσ,Lσ
|∇v (y)|2 dy
ΠP0 b̃(u,v) − b̃h(ΠRT0u, Π̃P0v),ψh
≤ C h ‖u‖H2
σ∈Eint
|ψLσ −ψKσ |
σ∈Eint
DKσ,Lσ
|∇v (y)|2 dy
so that
ΠP0b̃(u,v)− b̃h(ΠRT0u, Π̃P0v),ψh
)∣∣∣ ≤ C h ‖u‖H2 ‖ψh‖1,h ‖v‖1. Using then definition (2.5), we get
the result.
Let v ∈ L∞ ∩H1 and u ∈ H1 with divu ≥ 0 a.e. in Ω. Integrating by parts one checks that
v · b̃(u,v) dx =
divu dx ≥ 0. The operator bh shares a similar property.
Proposition 4.2. Let uh ∈ RT0 such that divuh ≥ 0. For all vh ∈ P0 we have
bh(uh,vh,vh) ≥ 0.
Proof. Remember that for all edges σ ∈ E inth , two triangles Kσ et Lσ share σ as an edge. We denote by Kσ
the one such that uσ · nKσ ,σ ≥ 0. Using the algebraic identity 2 a (a− b) = a
2 − b2 + (a − b)2 we deduce from
(2.14)
2 bh(uh,vh,vh) = 2
σ∈Eint
|σ|vKσ · (vKσ − vLσ) (uh · nKσ,σ)
σ∈Eint
|vKσ|
2 − |vLσ |
2 + |vKσ − vLσ |
(uh · nKσ,σ)
so that 2 bh(uh,vh,vh) ≥
σ∈Eint
|vKσ|
2 − |vLσ |
(uh · nKσ,σ). This sum can be written as a sum over
the triangles of the mesh. We get
2 bh(uh,vh,vh) ≥
|vKσ |
σ∈EK∩E
|σ| (uh · nKσ,σ).
Using finally the divergence formula we get
2 bh(uh,vh,vh) ≥
|K| |vK |
divuh dx ≥ 0.
The following result states that the operator bh is stable for suitable norms.
8 TITLE WILL BE SET BY THE PUBLISHER
Proposition 4.3. There exists a constant C > 0 such that for all vh ∈ P0, wh ∈ P0, uh ∈ P0 satisfying
divuh = 0
|bh(uh,vh,vh)| ≤ C |uh| ‖vh‖h ‖vh‖h.
Proof. For all triangle K ∈ Th and all edge σ ∈ EK ∩ E
h , we have
(uh · nK,σ)
σ vK + (uh · nK,σ)
σ vLσ = (uh · nK,σ)σ vK − |(uh · nK,σ)σ| (vLσ − vK).
Using this splitting, we deduce from (2.14) bh(uh,vh,wh) = S1 + S2 with
vK ·wK
σ∈EK∩E
|σ| (uh · nK,σ)σ ,
S2 = −
σ∈EK∩E
|σ| |(uh · nK,σ)σ| (vLσ − vK).
By writing the sum over the edges as a sum over the triangles we have
S2 = −
σ∈Eint
|σ| |(uh · nK,σ)σ| (vLσ − vK) · (wLσ −wK).
Using the Cauchy-Schwarz inequality we get
|S2| ≤ h ‖uh‖∞
σ∈Eint
|vLσ − vKσ |
1/2 
σ∈Eint
|wLσ −wKσ |
Since uh ∈ RT0 we have [5] the inverse inequality h ‖uh‖∞ ≤ C |uh|. Using (2.1) and (2.2) we get
σ∈Eint
|vLσ − vKσ |
2 ≤ C
σ∈Eint
τσ |vLσ − vKσ |
2 ≤ C ‖vh‖
and in a similar way
σ∈Eint
|wLσ −wKσ |
2 ≤ C ‖wh‖
h. Thus |S2| ≤ C |uh| ‖vh‖h ‖wh‖h. On the other hand,
according to the divergence formula
|K| (vK ·wK)
divuh dx = 0.
By gathering the estimates for S1 and S2 we get the result.
4.2. Properties of the discrete divergence
The operators gradient and divergence are adjoint: if q ∈ H1 , v ∈ H1 with v · n|∂Ω = 0, we get (v,∇q) =
−(q, divv) by integrating by parts. For ∇h and divh we state the following.
Proposition 4.4. For all vh ∈ P0 and qh ∈ P
1 we have: (vh,∇hqh) = −(qh, divh vh).
Proof. According to (2.8)
(vh,∇hqh) =
|K|vK · ∇hqh|K =
|σ|qh(xσ)nK,σ
TITLE WILL BE SET BY THE PUBLISHER 9
By writing this sum as a sum over the edges we get
(vh,∇hqh) = −
σ∈Eint
|σ| qh(xσ) (vLσ − vKσ) · nKσ,σ +
σ∈Eext
|σ| qh(xσ)vKσ · nKσ,σ. (4.1)
On the other hand, using a quadrature formula
−(qh, divh vh) = −
qh(xσ) (divh vh)(xσ).
By writing this sum as a sum over the edges of the mesh we get
−(qh, divh vh) = −
σ∈Eint
( |Kσ|
qh(xσ) (divhvh)(xσ)−
σ∈Eext
qh(xσ) (divh vh)(xσ).
Using definition (2.12) and comparing with (4.1) we get the result.
The divergence operator and the spaces L20, H
0 satisfy the following property, called inf-sup (or Babuška-Brezzi)
condition (see [9] for example). There exists a constant C > 0 such that
(q, divv)
‖v‖1|q|
≥ C. (4.2)
We will now show that the operator divh and the spaces P0 ∩ L
0, P0 satisfy an analogous property. The proof
uses the following lemma.
Lemma 4.1. There exists a constant C > 0 such that
∀ qh ∈ P
1 ∩ L
0 , sup
vh∈P0\{0}
(qh, divh vh)
‖vh‖h
≥ C h ‖qh‖1,h.
Proof. If qh = 0 the result is trivial. Let qh ∈ P
0\{0}. Let vh = ∇hqh ∈ P0\{0}. Using proposition 4.4
we have
−(qh, divhvh) = (vh,∇hqh) = |∇hqh|
2 = |∇hqh| |vh|.
Using (2.3) and (2.4) we get −(qh, divhvh) ≥ C h ‖qh‖1,h ‖vh‖h.
We now state the result.
Proposition 4.5. There exists a constant C > 0 such that for all qh ∈ P
1 ∩ L
vh∈P0\{0}
(qh, divh vh)
‖vh‖h
≥ C |qh|.
Proof. If qh = 0 the result is trivial. Let qh ∈ P
0\{0}. According to (4.2) there exists v ∈ H
0 such that
divv = −qh and ‖v‖1 ≤ C |qh|. (4.3)
We set vh = ΠPc
v. We want to estimate −
qh, divh(ΠP0vh)
. Since ∇hqh ∈ P0 we deduce from proposition
qh, divh(ΠP0vh)
= (ΠP0vh,∇hqh) = (vh,∇hqh).
By splitting the last term we get
qh, divh(ΠP0vh)
= (v,∇hqh)− (v − vh,∇hqh). (4.4)
10 TITLE WILL BE SET BY THE PUBLISHER
We bound the right-hand side of (4.4). Using (2.7) and (4.3) we have
|v − vh| = |v −ΠPc
v| ≤ C h ‖v‖1 ≤ C h |qh|.
Thus, using the Cauchy-Schwarz inequality, we get
|(v − vh,∇hqh)| ≤ C h |qh| |∇hqh| ≤ C h |qh| ‖qh‖1,h.
We estimate the other term as follows. Integrating by parts we get
(v,∇hqh) = −(qh, divv) +
qh (v · nK,σ) dσ.
We have −(qh, divv) = |qh|
2 thanks to (4.3). On the other hand
qh (v · nK,σ) dσ =
σ∈Eint
qh (v · nKσ,σ) dσ
since v|∂Ω = 0. Using [2] p.269 and (4.3) we have
∣∣∣∣∣
qh (v · nK,σ) dσ
∣∣∣∣∣
≤ C h ‖v‖1 ‖qh‖1,h ≤ C h |qh| ‖qh‖1,h.
Hence we get (v,∇hqh) ≥ (|qh| − C h ‖qh‖1,h) |qh|. Thus we deduce from (4.4)
qh, divh(ΠP0vh)
≥ (|qh| − C h ‖qh‖1,h) |qh|. (4.5)
We now introduce the norm ‖.‖h. We have vh = ΠPc
v ∈ Pc1 ⊂ H
1. From [5] p. 776 we deduce ‖ΠP0vh‖h ≤
C ‖vh‖1. Since ΠPc
is stable for the H1 norm, using (4.3), we get
‖vh‖1 = ‖ΠPc
v‖1 ≤ ‖v‖1 ≤ C |qh|.
Therefore ‖ΠP0vh‖h ≤ C |qh|. Using this inequality in (4.5) we obtain that there exists C1 > 0 and C2 > 0
such that
qh, divh(ΠP0vh)
≥ (C1 |qh| − C2 h ‖qh‖1,h) ‖ΠP0vh‖h.
We deduce from this
vh∈P0\{0}
(qh, divh vh)
‖vh‖h
≥ C1 |qh| − C2 h ‖qh‖1,h.
Let us combine this result with lemma 4.1. Since
∀ t ≥ 0 , max
C t , C1 |qh| − C2 t
C + C2
|qh| ,
we finally get the result.
4.3. Properties of the discrete laplacian
We recall from [14] the coercivity of the laplacian operator.
Proposition 4.6. For all uh ∈ P0 and vh ∈ P0 we have
−(∆̃huh,uh) = ‖uh‖
h , −(∆̃huh,vh) ≤ ‖uh‖h ‖vh‖h.
TITLE WILL BE SET BY THE PUBLISHER 11
5. Stability of the scheme
We first prove an estimate for the computed velocity (theorem 5.1). We show a similar result for the
increments in time (lemma 5.2). Using the inf-sup condition (proposition 4.5), we infer from it some estimates
on the pressure (theorem 5.2).
Lemma 5.1. For all m ∈ {0, . . . , N} et n ∈ {0, . . . , N} we have
(umh ,∇hp
h) = 0 , |u
2 − |ũmh |
2 + |umh − ũ
2 = 0.
Proof. First, using propositions 3.1 and 4.4, we get (umh ,∇hp
h) = −(p
h, divhu
h ) = 0. Also, we deduce from
(3.4)
2 (umh ,u
h − ũ
h ) = −
umh ,∇h(p
h − p
Using the algebraic identity 2 a (a− b) = a2 − b2 + (a− b)2 we get
2 (umh ,u
h − ũ
h ) = |u
2 − |ũmh |
2 + |umh − ũ
2 = 0.
We introduce the following hypothesis on the initial data.
(H1) There exists C > 0 such that |u0h|+ |u
h|+ k|∇hp
h| ≤ C.
Hypothesis (H1) is fulfilled if we set u0h = ΠRT0u0 and we use a semi-implicit Euler scheme to compute u
We have the following stability result.
Theorem 5.1. We assume that the initial values of the scheme fulfill (H1). For all m ∈ {2, . . . , N} we have
|umh |
2 + k
‖ũnh‖
h ≤ C. (5.1)
Proof. Let m ∈ {2, . . . , N} and n ∈ {1, . . . ,m− 1}. Taking the scalar product of (3.2) with 4 k ũn+1h we get
3 ũn+1h − 4u
h + u
, 4 k ũn+1h
(∆̃hũ
h , ũ
+4 k bh(2u
h − u
h , ũ
h , ũ
h ) + 4 k (∇hp
h, ũ
h ) = 4 k (f
h , ũ
h ). (5.2)
First of all, using lemma 5.1 and proceeding as in [10], we get
ũn+1h ,
3 ũn+1h − 4u
h + u
= |un+1h |
2 − |unh|
2 + |2un+1h − u
2 − |2unh − u
+ |un+1h − 2u
h + u
2 + 6 |ũn+1h − u
According to proposition 4.6 we have − 4 k
(∆̃hũ
h , ũ
h ) =
‖ũn+1h ‖
h. Also, according to lemma 5.1 and
(3.4)
4 k (∇hp
h, ũ
h ) = 4 k (∇hp
h, ũ
h − u
(|∇pn+1h |
2 − |∇pnh|
2 − |∇pn+1h −∇p
12 TITLE WILL BE SET BY THE PUBLISHER
Multiplying equation (3.4) by 4 k∇h(p
h − p
h) and using the Young inequality we get
|∇(pn+1h − p
2 ≤ 3 |un+1h − ũ
According to proposition 4.2, we have 4 k bh(2u
h − u
h , ũ
h , ũ
h ) ≥ 0. At last using the Cauchy-Schwarz
inequality, (2.3) and (3.1) we have
4 k (fn+1h , ũ
h ) ≤ 4 k |f
h | |ũ
h | ≤ C k ‖f‖C(0,T ;L2) ‖ũ
h ‖h.
Using the Young inequality we get
4 k (fn+1h , ũ
h ) ≤ 3 k ‖ũ
h + C k ‖f‖
C(0,T ;L2).
Thus we deduce from (5.2)
|un+1h |
2 − |unh |
2 + |2un+1h − u
2 − |2unh − u
2 + |un+1h − 2u
h + u
+3 |ũn+1h − u
2 + k ‖ũn+1h ‖
(|∇hp
2 − |∇hp
2) ≤ C k.
Summing from n = 1 to m− 1 we have
|umh |
2 + |2umh − u
2 + 3
|ũn+1h − u
2 + k
‖ũn+1h ‖
≤ C + 4 |u1h|
2 + |2u1h − u
2 + k2 |∇hp
Using hypothesis (H1) we get (5.1).
We now want to estimate the computed pressure. From now on, we make the following hypothesis on the data
f ∈ C(0, T ;L2) , ft ∈ L
2(0, T ;L2) , u0 ∈ H
2 ∩H10 , divu0 = 0.
One shows that if the data u0 and f fulfill a compatibility condition [11] there exists a solution (u, p) to the
equations (1.1)–(1.2) such that
u ∈ C(0, T ;H2) , ut ∈ C(0, T ;L
2) , ∇p ∈ C(0, T ;L2).
We introduce the following hypothesis on the initial values of the scheme: there exists a constant C > 0 such
(H2) |u0h − u0|+
‖u1h − u(t1)‖∞ + |p
h − p(t1)| ≤ C h , |u
h − u
h| ≤ C k.
One checks easily that this hypothesis implies (H1). We have the following result.
Lemma 5.2. We assume that the initial values of the scheme fulfill (H2). Then there exists a constant C > 0
such that for all m ∈ {1, . . . , N}
|umh − u
h | ≤ C.
Proof. Using proposition 4.1 one proceeds as in [14]. The difference lies in the way we bound the term ∇hp
We use the splitting
p1h = (p
h −ΠPnc1 p(t1)) + (ΠP
p(t1)− p(t1)) + p(t1).
TITLE WILL BE SET BY THE PUBLISHER 13
Using an inverse inequality [2] we have
p1h −ΠPnc1 p(t1)
)∣∣ ≤
∣∣p1h −ΠPnc1 p(t1)
(∣∣p1h − p(t1)
∣∣p(t1)−ΠPnc
p(t1)
∣∣) .
Using (2.10) and hypothesis (H2) we get
p1h −ΠPnc1 p(t1)
)∣∣ ≤ C ‖p(t1)‖1 ≤ C ‖p‖C(0,T ;H1).
According to (2.10) we also have
∣∣∇h(p(t1)−ΠPnc
p(t1))
∣∣ ≤ C ‖p(t1)‖1 ≤ C ‖p‖C(0,T ;H1). Lastly |∇p(t1)| ≤
‖p‖C(0,T ;H1). Thus we get |∇hp
h| ≤ C.
Theorem 5.2. We assume that the initial values of the scheme fulfull (H2). There exists a constant C > 0
such that for all m ∈ {2, . . . , N}
|pnh|
2 ≤ C.
Proof. Let m ∈ {2, . . . , N}. We set n = m− 1. Using the inf-sup condition (4.5) and proposition 4.4, we get
that there exists vh ∈ P0\{0} such that
C ‖vh‖h |p
h | ≤ −(p
h , divh vh) = (∇hp
h ,vh). (5.3)
Plugging (3.4) into (3.2) we have
h = −
3un+1h − 4u
h + u
∆̃hũ
h − b̃h(2u
h − u
h , ũ
h ) + f
so that
h ,vh) = −
3un+1h − 4u
h + u
∆̃hũ
h ,vh
− bh(2u
h − u
h , ũ
h ,vh) + (f
h ,vh).
Thanks to proposition 4.3 and theorem 5.1 we have
∣∣bh(2unh − u
h , ũ
h ,vh)
2 |unh|+ |u
‖ũn+1h ‖h ‖vh‖h ≤ C ‖ũ
h ‖h ‖vh‖h.
According to proposition 4.6 we have
∆̃hũ
h ,vh
≤ ‖ũn+1h ‖h ‖vh‖h. Using the Cauchy-Schwarz inequality,
(2.3) and (3.1) we have
(fn+1h ,vh) ≤ |f
h | |vh| ≤ C |vh| ≤ C ‖vh‖h
and in a similar way
3un+1h − 4u
h + u
)∣∣∣∣ ≤ C
3un+1h − 4u
h + u
∣∣∣∣ ‖vh‖h.
Thus we get
h ,vh) ≤ C + C
|3un+1h − 4u
h + u
+ ‖ũn+1h ‖h
‖vh‖h.
By comparing with (5.3) we get
|pn+1h | ≤ C + C
|3un+1h − 4u
h + u
+ ‖ũn+1h ‖h
14 TITLE WILL BE SET BY THE PUBLISHER
Squaring and summing from n = 1 to m− 1 we obtain
|pnh|
2 ≤ C + C k
|3un+1h − 4u
h + u
+ C k
‖ũn+1h ‖
The last term on the right-hand side is bounded, thanks to theorem 5.1. And since
3un+1h − 4u
h + u
h = 3(u
h − u
h)− (u
h − u
h ) = 3 δu
h − δu
we deduce from lemma 5.2
|3un+1h − 4u
h + u
≤ C k
|δunh |
References
[1] S. Boivin, F. Cayre, J. M Herard, A finite volume method to solve the Navier-Stokes equations for incompressible flows on
unstructured meshes, Int. J. Therm. Sci. 39 (2000) 806–825.
[2] S. C. Brenner and L. R. Scott, The mathematical theory of finite element methods, Springer, 2002.
[3] F. Brezzi and M. Fortin, Mixed and hybrid finite element methods, Springer-Verlag, 1991.
[4] J. Chorin, On the convergence of discrete approximations to the Navier-Stokes equations, Math. Comp. 23 (1969) 341–353.
[5] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods, P.G. Ciarlet and J.L. Lions eds, North-Holland, 2000.
[6] R. Eymard and R. Herbin, A staggered finite volume scheme on general meshes for the Navier-Stokes equations in two space
dimensions, Int.J. Finite Volumes (2005).
[7] R. Eymard, J. C. Latché and R. Herbin, Convergence analysis of a colocated finite volume scheme for the incompressible
Navier-Stokes equations on general 2 or 3D meshes, SIAM J. Numer. Anal. 45(1) (2007) 1–36.
[8] S. Faure, Stability of a colocated finite volume scheme for the Navier-Stokes equations, Num. Meth. PDE 21(2) (2005) 242–271.
[9] V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes equations: Theory and Algorithms, Springer, 1986.
[10] J. L. Guermond, Some implementations of projection methods for Navier-Stokes equations, M2AN 30(5) (1996) 637–667.
[11] J. G. Heywood and R. Rannacher, Finite element approximation of the nonstationary Navier-Stokes problem. I. Regularity of
solutions and second-order error estimates for spatial discretization, SIAM J. Numer. Anal. 19(26) (1982) 275–311.
[12] D. Kim and H. Choi, A second-order time-accurate finite volume method for unsteady incompressible flow on hybrid unstruc-
tured grids, J. Comp. Phys. 162 (2000) 411–428.
[13] R. Temam, Sur l’approximation de la solution des équations de Navier-Stokes par la méthode de pas fractionnaires II, Arch.
Rat. Mech. Anal. 33 (1969) 377–385.
[14] S. Zimmermann, Stability of a colocated finite volume for the incompressible Navier-Stokes equations, arXiv :0704.0772 (2006).
[15] S. Zimmermann, Étude et implémentation de méthodes de volumes finis pour les fluides incompressibles, PhD, Blaise Pascal
University, France (2006).
