Stability of a colocated finite volume scheme for the incompressible
  Navier-Stokes equations
arXiv:0704.0772v1  [math.NA]  5 Apr 2007
INTERNATIONAL JOURNAL OF c© 2007 Institute for Scientific
NUMERICAL ANALYSIS AND MODELING Computing and Information
Volume 1, Number 1, Pages 1–18
STABILITY OF A COLOCATED FINITE VOLUME SCHEME
FOR THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS
SÉBASTIEN ZIMMERMANN
(Communicated by Jean-Luc Guermond)
Abstract. We introduce a finite volume scheme for the two-dimensional incom-
pressible Navier-Stokes equations. We use a triangular mesh. The unknowns
for the velocity and pressure are both piecewise constant (colocated scheme).
We use a projection (fractional-step) method to deal with the incompressibility
constraint. We prove that the differential operators in the Navier-Stokes equa-
tions and their discrete counterparts share similar properties. In particular, we
state an inf-sup (Babuška-Brezzi) condition. We infer from it the stability of
the scheme.
Key Words. Incompressible fluids, Navier-Stokes equations, projection meth-
ods and finite volume.
1. Introduction
We consider the flow of an incompressible fluid in a open bounded 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 condition
u|∂Ω = 0 , u|t=0 = u0.
The terms ∆u and (u ·∇)u are respectively associated with the physical phenom-
ena of diffusion and convection. The Reynolds number Re measures the influence
of convection in the flow. For equations (1.1)–(1.2), finite element and finite dif-
ference methods are well known and mathematical studies are available (see [10]
for example). Numerous computations have also been conducted with finite vol-
ume schemes (e.g. [14] and [1]). However, in this case, few mathematical results
are available. Let us cite Eymard and Herbin [7] and Eymard, Latché and
Herbin [8]. 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 [15]. This is the case in Faure
[9]. In this work, however, the mesh is made of squares, so that the geometry of
the problem is limited. Therefore, we introduce in what follows a finite volume
scheme on triangular meshes for equations (1.1)–(1.2), using a projection method.
An interesting feature of this scheme is that the unknowns for the velocity and
Received by the editors April 1, 2007 and, in revised form, April 1, 2007.
2000 Mathematics Subject Classification. 76M12, 76B99.
http://arxiv.org/abs/0704.0772v1
2 S. ZIMMERMANN
pressure are both piecewise constant (colocated scheme). It leads to an economic
computer storage, and allows an easy generalization of the scheme to the 3D case.
The layout of the article is the following. We first introduce (section 2) some no-
tations 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 for the scheme, we prove 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.
Also, the discrete gradient operator is a consistent approximation of its continuous
counterpart. The discrete operator for the convection term is positive, stable and
consistent. The discrete operator for the divergence satisfies 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 ∈ L2 ;
q(x) dx = 0}. Their vectorial
counterparts are (L2, |.|) and (L∞, ‖.‖∞) with L2 = (L2)2 and L∞ = (L∞)2. For
k ∈ N∗, (Hk, ‖·‖k) is the usual Sobolev space. Its vectorial counterpart is (Hk, ‖.‖k)
with Hk = (Hk)2. For k = 1, the functions of H1 with a null trace on the
boundary form the spaceH10. Also, we set ∇u = (∇u1,∇u2)T if u = (u1, u2) ∈ H1.
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 continous (resp.
square integrable). The norms ‖.‖C(0,T ;X) and ‖.‖L2(0,T ;X) are defined respectively
by ‖g‖C(0,T ;X) = supt∈[0,T ] |g(t)| and ‖g‖L2(0,T ;X) =
|g(t)|2 ds
. 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 Ω: Ω = ∪K∈ThK. For each
triangle K ∈ Th, we denote by |K| its area and EK the set of his edges. If σ ∈ EK ,
nK,σ is the unit vector normal to σ pointing outward of K.
The set of edges of the mesh is Eh = ∪K∈ThEK . The length of an edge σ ∈ Eh is |σ|
and its middle point xσ. The set of edges located inside Ω (resp. on its boundary)
is E inth (resp. Eexth ): Eh = E inth ∪ Eexth . If σ ∈ E inth , Kσ and Lσ are the triangles
sharing σ as an edge. If σ ∈ Eexth , only the triangle Kσ inside Ω is defined.
We denote by xK the circumcenter of a triangle K. We assume that the measure
of all interior angles of the triangles of the mesh are below π
, so that xK ∈ K. If
σ ∈ E inth (resp. σ ∈ Eexth ) we set dσ = d(xKσ ,xLσ ) (resp. dσ = d(xσ ,xKσ)). We
define for all edge σ ∈ Eh
(2.1) τσ =
The maximum circumradius of the triangles of the mesh is h. We assume ([6] p.
776) that there exists C > 0 such that
∀σ ∈ Eh, d(xKσ , σ) ≥ C|σ| and |σ| ≥ Ch.
It implies that there exists C > 0 such that
(2.2) ∀σ ∈ Eh , τσ ≥ C ,
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 3
and for all triangles K ∈ Th we have (with σ ∈ EK and hK,σ the matching altitude)
(2.3) |K| = 1
|σ|hK,σ ≥
|σ| d(xK ,xσ) ≥ C h2.
Lastly, if K ∈ Th and L ∈ Th are two triangles sharing the edge σ ∈ E inth , we define
αK,L =
d(xL,xσ)
d(xK ,xL)
Let us notice that αK,L ∈ [0, 1] and αK,L + αL,K = 1.
2.2. The discrete spaces. We first define
P0 = {q ∈ L2 ; ∀K ∈ Th, q|K is a constant} , P0 = (P0)2.
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⊂ H1, we define the discrete
equivalent of a H1 norm as follows. For all vh ∈ P0 we set
(2.4) ‖vh‖h =
σ∈Eint
τσ |vLσ − vKσ |2 +
σ∈Eext
τσ |vKσ |2
where τσ is given by (2.1). We have [6] a Poincaré-like inequality for P0: there
exists C > 0 such that for all vh ∈ P0
(2.5) |vh| ≤ C ‖vh‖h.
We also have the following inverse inequality.
Proposition 2.1. There exists a constant C > 0 such that for all vh ∈ P0
h ‖vh‖h ≤ C |vh|.
Proof. According to (2.4)
h2 ‖vh‖2h =
σ∈Eint
h2 τσ |vLσ − vKσ |2 +
σ∈Eext
h2 τσ |vKσ |2.
We deduce from (2.2) and (2.3) that h2 τσ ≤ C |Kσ| and h2 τσ ≤ C |Lσ|. Thus,
since |vLσ − vKσ |2 ≤ 2
|vLσ |2 + |vKσ |2
, we get
h2 ‖vh‖2h ≤ C
σ∈Eint
|Kσ| |vKσ |2 + |Lσ| |vLσ |2
σ∈Eext
|Kσ| |vKσ |2.
Hence h2 ‖vh‖2h ≤ C
|K| |vK |2 ≤ C |vh|2.
From the norm ‖.‖h we deduce a dual norm. For all vh ∈ P0 we set
(2.6) ‖vh‖−1,h = sup
(vh,ψh)
‖ψh‖h
For all uh ∈ P0 and vh ∈ P0 we have (uh,vh) ≤ ‖uh‖−1,h ‖vh‖h. Now we introduce
some operators on P0 and P0. We define the projection operator ΠP0 : L
2 → P0
as follows. For all w ∈ L2, ΠP0w ∈ P0 is given by
(2.7) ∀K ∈ Th , (ΠP0w)|K =
w(x) dx.
We easily check that for all w ∈ L2 and vh ∈ P0 we have (ΠP0w,vh) = (w,vh).
It implies that ΠP0 is stable for the L
2 norm. We define also the interpolation
operator Π̃P0 : H
2 → P0. For all q ∈ H2, Π̃P0q ∈ P0 is given by
∀K ∈ Th , Π̃P0q|K = q(xK).
4 S. ZIMMERMANN
According to the Sobolev embedding theorem, q ∈ H2 is a.e. equal to a continuous
function. Therefore the definition above makes sense. We also set Π̃P0 = (Π̃P0 )
The operator Π̃P0 (resp. Π̃P0) is naturally stable for the L
∞ (resp. L∞) norm.
One also checks ([2] and [16]) that there exists C > 0 such that
(2.8) |v −ΠP0v| ≤ C h ‖v‖1 , |q − Π̃P0q| ≤ C h ‖q‖2
for all v ∈ H1 and q ∈ H2.
We introduce the finite element spaces
P d1 = {v ∈ L2 ; ∀K ∈ Th, v|K is affine} ,
1 = {vh ∈ P d1 ; ∀σ ∈ E inth , vh|Kσ(xσ) = vh|Lσ(xσ) ,
Pc1 = {vh ∈ (P d1 )2 ; vh is continuous and vh|∂Ω = 0}.
We have Pc1 ⊂ H10. We define the projection operator ΠPc1 : H
0 → Pc1. For all
v = (v1, v2) ∈ H10, ΠPc1v = (v
h) ∈ Pc1 is given by
∀φh = (φ1h, φ2h) ∈ Pc1 ,
∇vih,∇φih) =
∇vi,∇φih).
The operator ΠPc
is stable for the H1 norm and ([2] p. 110) there exists C > 0
such that for all v ∈ H1
(2.9) |v −ΠPc
v| ≤ C h ‖v‖1.
Let us address now the space Pnc1 . If qh ∈ Pnc1 , we have usually ∇qh 6∈ L2. Thus
we define the operator ∇̃h : Pnc1 → P0 by setting for all qh ∈ Pnc1 and all K ∈ Th
∇̃hqh|K =
∇qh dx.
The associated norm is given by
‖qh‖1,h =
|qh|2 + |∇̃hqh|2
We also have a Poincaré inequality: there exists C > 0 such that for all qh ∈ Pnc1 ∩L20
(2.10) |qh| ≤ C |∇̃hqh|.
We define the projection operator ΠPnc
. For all qh ∈ Pnc1 , ΠPnc1 qh is given by
(2.11) ∀φ ∈ L2 , (ΠPnc
qh, φ) = (qh, φ).
We have the following result.
Proposition 2.2. If qh ∈ P0, ΠPnc
qh is given by
∀σ ∈ E inth , (ΠPnc1 qh)(xσ) =
|Kσ|+ |Lσ|
qKσ +
|Kσ|+ |Lσ|
qLσ ,
∀σ ∈ Eexth , (ΠPnc1 qh)(xσ) = qKσ .
Proof. For all edge σ ∈ Eh, we define the function ψσ ∈ Pnc1 by setting
ψσ(xσ′ ) =
1 if σ = σ′,
0 otherwise.
Let us notice that ψσ vanishes outside Kσ ∪ Lσ if σ ∈ E inth and outside Kσ if
σ ∈ Eexth . Let σ ∈ E inth . Using a quadrature formula we get
(ΠPnc
qh, ψσ) =
(ΠPnc
qh)(xσ)
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 5
(qh, ψσ) = qKσ
+ qLσ
For an edge σ ∈ Eexth we have (ΠPnc1 qh, ψσ) =
(ΠPnc
qh)(xσ) and (qh, ψσ) =
. By plugging these equations into (2.11) with φ = ψσ, we get the result.
We finally introduce the Raviart-Thomas spaces
={vh ∈ Pd1 ; ∀σ ∈ EK , vh|K · nK,σ is a constant, and vh · n|∂Ω = 0} ,
RT0 ={vh ∈ RTd0 ; ∀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
1 → RT0. For all v ∈ H1, ΠRT0v ∈ RT0 is given by
(2.12) ∀K ∈ Th , ∀σ ∈ EK , (ΠRT0v · nK,σ)σ =
v dσ.
One checks [3] that there exists C > 0 such that for all v ∈ H1
(2.13) |v −ΠRT0v| ≤ C h ‖v‖1.
The following result will be useful.
Proposition 2.3. For all v ∈ H1 such that divv = 0, we have ΠRT0v ∈ P0.
Proof. Let vh = ΠRT0v and K ∈ Th. According to [3] there exists aK ∈ R2 and
bK ∈ R such that: ∀x ∈ K , vh(x) = aK + bK x. Thus divvh|K = 2 bK . On the
other hand, according to the divergence formula and (2.12)
divv dx =
v · n dγ =
vh · n dγ =
divvh dx.
Hence bK = 0 and we get: ∀x ∈ K , vh(x) = aK .
2.3. The discrete operators. The equations (1.1)–(1.2) use the differential op-
erators gradient, divergence and laplacian. Using the spaces of section 2.2 we define
their discrete counterparts. The discrete gradient ∇h : P0 → P0 is built using a
linear interpolation on the edges of the mesh (see [16] for details). This kind of
construction has also be considered in [5]. We set for all qh ∈ P0 and all K ∈ Th
∇h qh|K =
σ∈EK∩E
αKσ,Lσ qKσ + αLσ,Kσ qLσ
σ∈EK∩E
|σ| qKσ nK,σ.(2.14)
We have the following result [16].
Proposition 2.4. If qh ∈ L20 is such that ∇hqh = 0, then qh = 0.
The discrete divergence operator divh : P0 → P0 is built so that it is adjoint to the
operator ∇h (proposition 4.6 below). We set for all qh ∈ P0 and all K ∈ Th
(2.15) divh vh|K =
σ∈EK∩E
αLσ ,Kσ vKσ + αKσ,Lσ vLσ
· nK,σ.
The first discrete laplacian ∆h : P0 → P0 ensures that the incompressibility con-
straint (1.2) is satisfied in a discrete sense (proposition 3.1). We set for all qh ∈ P0
(2.16) ∆hqh = divh(∇hqh).
6 S. ZIMMERMANN
The second discrete laplacian ∆̃h : P0 → P0 is the usual operator in finite volume
schemes [6]. We set for all vh ∈ P0 and all K ∈ Th
∆̃hvh|K =
σ∈EK∩E
τσ (vLσ − vKσ )−
σ∈EK∩E
τσ vKσ .
In order to approximate the convection term (u ·∇)u in (1.1) we define a bilinear
form b̃h : P0 ×P0 → P0 using the well-known upwind scheme ([6] p. 766). For all
uh ∈ P0, vh ∈ P0, and all K ∈ Th we have
(2.17) b̃h(uh,vh)
σ∈EK∩E
(uσ · nK,σ)+ vK + (uσ · nK,σ)− vLσ
We have set uσ = αLσ ,Kσ uKσ + αKσ ,Lσ uLσ and a
+ = max(a, 0), a− = min(a, 0)
for all a ∈ R. Lastly, we define the trilinear form bh : P0 × P0 × P0 → R2 as
follows. For all uh ∈ P0, vh ∈ P0, wh ∈ P0, we set
(2.18) bh(uh,vh,wh) =
|K|wK · b̃h(uh,vh)
3. The scheme
We have defined in section 2 the discretization in space. We now have to define
a discretization in time, and treat the incompressibility constraint (1.2). We use a
projection method to this end. This kind of method has been introduced byChorin
[4] and Temam [15]. 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 umh ≃ u(tm), which fulfills the incompressibilty constraint (1.2) in a
discrete sense. The algorithm goes as follows.
First, for all m ∈ {0, . . . , N}, we set fmh = ΠP0 f(tm). Since the operator ΠP0 is
stable for the L2-norm we get
(3.1) |fmh | = |ΠP0 f(tm)| ≤ |f(tm)| ≤ ‖f‖C(0,T ;L2).
We start with the initial values
u0h ∈ P0 ∩RT0 , u1h ∈ P0 ∩RT0 p1h ∈ P0 ∩ L20.
For all n ∈ {1, . . . , N}, (ũn+1h , p
h ) is deduced from (ũ
h) as follows.
• ũn+1h ∈ P0 is given by
(3.2)
3 ũn+1h − 4unh + u
∆̃hũ
h +b̃h(2u
h−un−1h , ũ
h )+∇hp
h = f
• pn+1h ∈ Pnc1 ∩ L20 is the solution of
(3.3) ∆h(p
h − p
divh ũ
• un+1h ∈ P0 is deduced by
(3.4) un+1h = ũ
∇h(pn+1h − p
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 7
Existence and unicity of a solution to equation (3.2) is classical ([6] for example).
Let us show that equation (3.3) has also a unique solution. Let qh ∈ P0 ∩ L20 such
that ∆hqh = 0. According to proposition 4.6 we have for all qh ∈ P0
−(∆hqh, qh) = −
divh(∇hqh), qh
= (∇hqh,∇hqh) = |∇hqh|2.
Therefore we have ∇hqh = 0. Using proposition 2.4 we get qh = 0. We have
thus proved the unicity of a solution for equation (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). Let us now prove that for all m ∈
{0, . . . , N}, umh fulfills (1.2) in a discrete sense.
Lemma 3.1. If vh ∈ RT0 ∩P0 then divh vh = 0.
Proof. Let K ∈ Th. Since vh ∈ RT0, definition (2.15) reads
divh vh|K =
|σ| (αLσ ,K + αK,Lσ )vK · nK,σ.
Since αKσ ,Lσ + αLσ,Kσ = 1 we conclude that
divh vh|K =
|σ|vK · nK,σ = vK ·
|σ|nK,σ
Proposition 3.1. For all m ∈ {0, . . . , N} we have divh umh = 0.
Proof. For m ∈ {0, 1} we have u0h ∈ P0 ∩ RT0 and u1h ∈ P0 ∩RT0. Applying
the lemma above we get the result. If m ∈ {2, . . . , N}, we apply the operator divh
to (3.3) and compare with (3.4).
4. Properties of the discrete operators
We prove 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) ∈ H1 we set
(4.1) 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
‖ΠP0b̃(u,v) − b̃h(ΠRT0u, Π̃P0v)‖−1,h ≤ C h ‖u‖2 ‖v‖1.
Proof. Let uh = ΠRT0u and vh = Π̃P0v. According to proposition 2.3 we have
uh ∈ P0. Let K ∈ Th. According to the divergence formula and (2.7) we have
ΠP0b̃(u,v)|K =
σ∈EK∩E
v (u · n) dσ.
On the other hand, let us rewrite b̃h(uh,vh). Let σ ∈ EK ∩ E inth . Setting
vK,Lσ =
vK si (uh · nK,σ)σ ≥ 0
vLσ si (uh · nK,σ)σ < 0
one checks that vK (uσ ·nK,σ)++vLσ (uσ ·nK,σ)− = vK,Lσ (uσ ·nK,σ). By definition
uσ · nK,σ = αLσ,K (uK · nK,σ) + αK,Lσ (uLσ · nK,σ) ; since uh ∈ RT0 we get
8 S. ZIMMERMANN
uσ · nK,σ = (αLσ ,K + αK,Lσ ) (uK · nK,σ) = (uK · nK,σ) = (uh · nK,σ)σ. Using at
last (2.12), we deduce from (2.17)
b̃h(uh,vh)|K =
σ∈EK∩E
vK,Lσ (u · nK,σ) dσ.
ΠP0 b̃(u,v) − b̃h(uh,vh)
σ∈EK∩E
(v − vK,Lσ ) (u · n) dσ.
Let ψh ∈ P0. We have
ΠP0b̃(u,v) − b̃h(uh,vh),ψh
σ∈EK∩E
(v − vK,Lσ ) (u · n) dσ
σ∈Eint
(ψKσ −ψLσ)
(v − vKσ ,Lσ) (u · n) dσ.(4.2)
Let σ ∈ E inth . We want to estimate the integral over σ. Since we work in a two-
dimensional domain, we have the Sobolev injection H2 ⊂ L∞. Thus
(v − vKσ,Lσ) (u · n) dσ
∣∣∣∣ ≤ ‖u‖L∞
|v−vKσ ,Lσ | dσ ≤ C ‖u‖2
|v−vKσ ,Lσ | dσ.
Let us first assume that v ∈ C1. We set
xKσ,Lσ =
xKσ si (uh · nK,σ)σ ≥ 0
xLσ si (uh · nK,σ)σ < 0
If x ∈ σ, we have the following Taylor expansion
v(x)− vKσ,Lσ = v(x)−v(xKσ ,Lσ) =
∇v (tx+(1− t)xKσ ,Lσ) (x−xKσ ,Lσ) dt.
We have |x−xKσ ,Lσ | ≤ h. Thus, integrating over σ and using the Cauchy-Schwarz
inequality, we get
|v − vKσ ,Lσ | dσ ≤
|∇v (tx+ (1− t)xKσ ,Lσ)|2 h
t dt dσ
We then use the change of variable (t,x) → y = tx + (1 − t)xKσ ,Lσ . Let Dσ be
the quadrilateral domain given by the endpoints of σ, xKσ and xLσ . The domain
[0, 1]× σ becomes DKσ,Lσ with
DKσ,Lσ =
Dσ ∩Kσ si (uh · nK,σ)σ ≥ 0
Dσ ∩ Lσ si (uh · nK,σ)σ < 0
For all t ∈ [0, 1] we have h
t ≤ h t ≤ C d(xKσ ,Lσ , σ) t thanks to the hypothesis on
the mesh. We check easily that d(xKσ ,Lσ , σ) t dt dσ = dy. Thus we get
|v − vKσ ,Lσ | dσ ≤ C h
DKσ,Lσ
|∇v (y)|2 dy
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 9
Since (C1)2 is dense in H2, this estimate still holds for v ∈ H2. Plugging this
estimate into (4.2) and using the Cauchy-Schwarz inequality we get
Π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
ΠP0 b̃(u,v)− b̃h(ΠRT0u, Π̃P0v),ψh
)∣∣∣ ≤ C h ‖u‖H2 ‖ψh‖1,h ‖v‖1. Using
then definition (2.6), we get the result.
Let us consider now the operator bh. Let u ∈ H1 and v ∈ L∞∩H1 with divu ≥ 0.
Integrating by parts we deduce from (4.1):
v · b̃(u,v) dx =
divu dx ≥ 0.
The discrete operator bh shares a similar property.
Proposition 4.2. Let uh ∈ P0 such that divh uh ≥ 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) = a2 − b2 + (a− b)2 we deduce from (2.18)
2 bh(uh,vh,vh) = 2
σ∈Eint
|σ|vKσ · (vKσ − vLσ ) (uσ · nKσ,σ)
σ∈Eint
|vKσ|2 − |vLσ |2 + |vKσ − vLσ |2
(uσ · nKσ,σ)
so that 2 bh(uh,vh,vh) ≥
σ∈Eint
|vKσ|2−|vLσ |2
(uσ ·nKσ,σ). This sum can
be written as a sum over the triangles of the mesh. We get
2 bh(uh,vh,vh) ≥
|vKσ |2
σ∈EK∩E
|σ| (uσ · nKσ,σ).
Using finally definition (2.15) we get
2 bh(uh,vh,vh) ≥
|K| |vK |2 (divh uh)|K ≥ 0.
The following result states that the operator bh is stable for suitable norms.
Proposition 4.3. There exists a constant C > 0 such that for all vh ∈ P0, wh ∈
P0, uh ∈ P0 satisfying divh uh = 0
|bh(uh,vh,vh)| ≤ C |uh| ‖vh‖h ‖vh‖h.
Proof. For all triangle K ∈ Th and all edge σ ∈ EK ∩ E inth , we have
(uσ · nK,σ)+ vK + (uσ · nK,σ)− vLσ = (uσ · nK,σ)vK − |(uσ · nK,σ)| (vLσ − vK).
This way, we deduce from (4.7) bh(uh,vh,wh) = S1 + S2 with
vK ·wK
σ∈EK∩E
|σ| (uσ · nK,σ) ,
S2 = −
σ∈EK∩E
|σ| |uσ · nK,σ| (vLσ − vK).
10 S. ZIMMERMANN
By writing the sum over the edges as a sum over the triangles we get
S2 = −
σ∈Eint
|σ| |uσ · nK,σ| (vLσ − vK) · (wLσ −wK).
Using the Cauchy-Schwarz inequality we get
|S2| ≤ h ‖uh‖∞
σ∈Eint
|vLσ − vKσ |2
1/2 
σ∈Eint
|wLσ −wKσ |2
Since uh ∈ P0 we have the inverse inequality [6] h ‖uh‖∞ ≤ C |uh|. Using (2.2)
and (2.4) we have
σ∈Eint
|vLσ − vKσ |2 ≤ C
σ∈Eint
τσ |vLσ − vKσ |2 ≤ C ‖vh‖2h
σ∈Eint
|wLσ −wKσ |2 ≤ C ‖wh‖2h. Therefore |S2| ≤ C |uh| ‖vh‖h ‖wh‖h. On
the other hand we deduce from definition (2.15)
|K| (vK ·wK) (divh uh)|K = 0.
By combining the estimates for S1 and S2 we get the result.
4.2. Properties of the discrete gradient.
Proposition 4.4. There exists a constant C > 0 such that for all qh ∈ P0:
h |∇hqh| ≤ C |qh|.
Proof. Using (2.14) and the Minkowski inequality, we have for all triangle K ∈ Th
|K| |∇hqh |K |2 ≤
σ∈EK∩E
6 |σ|2
(q2K + q
σ∈EK∩E
6 |σ|2
q2K .
Let us sum over K ∈ Th. Since |σ| ≤ h, using (2.3), we get
|∇hqh|2 ≤
σ∈EK∩E
|K| q2K + |Lσ| q2Lσ
σ∈EK∩E
|K| q2K
Thus h2 |∇hqh|2 ≤ C
|K| q2K ≤ C |qh|2.
We now prove that ∇h is a consistent approximation of the gradient.
Proposition 4.5. There exists a constant C > 0 such that for all q ∈ H2
|ΠP0(∇q)−∇h(Π̃P0q)| ≤ C h ‖q‖2.
Proof. Let K ∈ Th. Using the gradient formula and definition (2.14) we get
ΠP0(∇q)−∇h(Π̃P0q)
∇q dx− |K| ∇h(Π̃P0q)
where we have set for all edge σ ∈ EK ∩ E inth
IσK =
αK,Lσ q(xK) + αLσ,K q(xLσ )
nK,σ dσ
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 11
and for all edge σ ∈ EK ∩ Eexth : IσK =
q − q(xK)
nK,σ dσ. Squaring and using
(2.3) we get
ΠP0(∇q)−∇h(Π̃P0q)
|IσK |2 ≤
|IσK |2.
Summing over the triangles K ∈ Th we get
(4.3)
∣∣∣ΠP0(∇q)−∇h(Π̃P0q)
|IσK |2.
We must estimate the integral terms IσK . Let K ∈ Th. Let us first assume that
q ∈ C2(Ω). Let σ ∈ EK ∩ E inth . For x ∈ σ we have the following Taylor expansions
q(xK) = q(x)+∇q(x) · (xK −x)+
H(q) (txK +(1− t)x)(xK −x) · (xK −x) t dt ,
q(xLσ) = q(x)+∇q(x)·(xLσ −x)+
H(q) (txLσ+(1−t)x)(xLσ−x)·(xLσ−x) t dt ,
∇q(x) = ∇q(xK)−
txK + (1− t)x
(xK − x) dt.
Plugging the last expansion into the two others and integrating over σ we get
(4.4)
q(xK)− q
dσ = |σ| ∇q(xK) · (xK − xσ)−AσK +BσK ,
(4.5)
q(xLσ )− q
dσ = |σ| ∇q(xK) · (xLσ − xσ)−AσLσ + B
We have set for T ∈ {Kσ, Lσ}
(4.6) AσT =
∇∇q (txT + (1 − t)x) (xT − x) dt dσ ,
(4.7) BσT =
H(q) (txT + (1− t)x)(xT − x) · (xT − x) t dt dσ.
One can bound these terms as in the proof of proposition 4.1. We get
(4.8) |AσT |2 ≤ C h2
|∇∇q (y)|2 dy , |BσT |2 ≤ C h4
|H(q)(y)|2 dy.
Now, let us multiply (4.4) by −αK,Lσ nK,σ, (4.5) by −αLσ,K nK,σ and sum the
equalities. Since αLσ,K + αK,Lσ = 1 we have
−αLσ,K
q(xK)− q
nK,σ dσ − αK,Lσ
q(xLσ )− q
nK,σ dσ
αKσ,Lσ q(xK,σ) + αLσ,Kσ q(xL,σ)
nK,σ dσ = I
On the other hand
−αK,Lσ (xK −xσ) ·nK,σ −αLσ,K (xLσ −xσ) ·nK,σ = −αK,Lσ αLσ,K (dσ − dσ) = 0.
Therefore we get IσK = −αLσ,K
AσK +B
nK,σ −αK,Lσ
AσLσ +B
nK,σ. Using
estimates (4.8) we obtain
|IσK |2 ≤ C h4
(|H(q)(y)|2 + |∇∇q(y)|2) dy.
12 S. ZIMMERMANN
We now consider the case σ ∈ EK ∩ Eexth . For x ∈ σ we have
q(xK) = q(x)+∇q(x) · (xK − x)+
H(q)(txK + (1− t)x)(xK − x) · (xK − x)tdt.
Multiplying by nK,σ and integrating over σ, we get −IσK = JσK nK,σ+BσK nK,σ with
JσK =
∇q(x) · (xK − x) dx. Since |xK − x| ≤ h if x ∈ σ, using a trace theorem,
we have
|JσK | ≤ C h2 ‖∇q‖L∞(σ) ≤ C h2
|∇q(y)|2 + |∇∇q(y)|2
By combining this estimate with (4.8), we get
|IσK |2 ≤ 2 |Jσ|2 + 2 |BσK |2 ≤ C h4
|H(q)(y)|2 dy
+ C h4
(|∇q(y)|2 + |∇∇q(y)|2) dy.
The space C2(Ω) is dense in H2. Therefore the bounds for IσK still hold for q ∈ H2.
Plugging these bounds into (4.3) we get the result.
4.3. Properties of the discrete divergence. The operators divergence and gra-
dient are adjoint: if q ∈ H1 and v ∈ H1 with v · n|∂Ω = 0, we get (v,∇q) =
−(q, divv) by integrating by parts. For ∇h and divh we state
Proposition 4.6. If vh ∈ P0 and qh ∈ P0 we have: (vh,∇hqh) = −(qh, divh vh).
Proof. Using (2.14) one checks that (vh,∇hqh) =
qK (S1 + S2 + S3) with
σ∈EK∩E
|σ|αK,Lσ vK · nK,σ , S2 =
σ∈EK∩E
|σ|αK,Lσ vLσ · nLσ,σ ,
and S3 =
σ∈EK∩E
|σ|vK · nK,σ. Since αK,Lσ + αLσ,K = 1 we have
σ∈EK∩E
|σ| (1− αLσ ,K)vK · nK,σ
σ∈EK∩E
|σ|vK · nK,σ −
σ∈EK∩E
|σ|αLσ,K vK · nK,σ.
Since nLσ,σ = −nK,σ, we also have
σ∈EK∩E
|σ|αK,Lσ vLσ · nLσ,σ = −
σ∈EK∩E
|σ|αK,Lσ vLσ · nK,σ.
Therefore
(vh,∇hqh) = −
σ∈EK∩E
|σ| (αL,Kσ vK +αK,Lσ vLσ) ·nK,Lσ +
|σ|vK ·nK,Lσ .
Using definition (2.15) we get
(vh,∇hqh) = −
|K| divh vh|K +
|σ|vK · nK,Lσ .
Since
|σ|vK · nK,Lσ = vK ·
|σ|nK,Lσ = 0 we obtain finally
(vh,∇hqh) = −
qK |K| divh vh|K = −(qh, divh vh).
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 13
The divergence operator and the spaces L20, H
0 satisfy the following property,
called inf-sup (or Babuška-Brezzi) condition (see [10] for example). There exists a
constant C > 0 such that
(4.9) inf
− (q, divv)
‖v‖1|q|
We will now prove that the operator divh and the spaces P0 ∩ L20, P0 satisfy an
analogous property. The proof is based on the following lemma.
Lemma 4.1. We assume that the mesh is uniform (i.e. the triangles of the mesh
are equilateral). Then we have for all qh ∈ P0
∇hqh = ∇̃h(ΠPnc
Proof. Since the mesh is uniform we have: ∀σ ∈ E inth , αKσ,Lσ = 12 . Let K ∈ Th.
Using definition (2.14) and the gradient formula we get
∇hqh − ∇̃h(ΠPnc
σ∈EK∩E
(qKσ + qLσ)nK,σ
σ∈EK∩E
|σ| qKσ nK,σ −
(ΠPnc
qh)nK,σ dσ.
Since qh ∈ P0 we deduce from proposition 2.2
qh dσ = |σ| (ΠPnc
qh)(xσ) =
(qKσ + qLσ) if σ ∈ E inth ,
|σ| qKσ if σ ∈ Eexth .
Plugging this into the equation above, we get ∇hqh|K = ∇̃h(ΠPnc
qh)|K .
Lemma 4.2. We assume that the mesh is uniform. There exists a constant C > 0
such that
∀ qh ∈ P0 ∩ L20 , sup
vh∈P0\{0}
− (qh, divh vh)
‖vh‖h
≥ C h ‖ΠPnc
qh‖1,h.
Proof. If qh = 0 the result is trivial. Let qh ∈ P0 ∩ L20\{0}. Let vh = ∇hqh ∈
P0\{0}. Using proposition 4.6 we have
−(qh, divhvh) = (vh,∇hqh) = |∇hqh|2 = |∇hqh| |vh|.
Let χΩ be the characteristic function of Ω. Putting ψ = χΩ in (2.11) we get
qh ∈ L20. So according to (2.10) and(4.1) we have
|∇hqh| =
∣∣∣∇̃h(ΠPnc
∣∣∣ ≥ C ‖ΠPnc
qh‖1,h.
On the other hand, according to proposition 2.1: |vh| ≥ C h ‖vh‖h. Therefore
−(qh, divhvh) ≥ C h ‖ΠPnc
qh‖1,h ‖vh‖h.
Proposition 4.7. We assume that the mesh is uniform. There exists a constant
C > 0 such that for all qh ∈ P0 ∩ L20
vh∈P0\{0}
− (qh, divh vh)
‖vh‖h
≥ C |ΠPnc
14 S. ZIMMERMANN
Proof. If qh = 0 the result is clear. Let qh ∈ P0 ∩ L20\{0}. According to (4.9)
there exists v ∈ H10 such that
(4.10) divv = −ΠPnc
qh and ‖v‖1 ≤ C |ΠPnc
We set vh = ΠPc
v. We want to estimate −
qh, divh(ΠP0vh)
. Since ∇hqh ∈ P0
we deduce from proposition 4.6
qh, divh(ΠP0vh)
= (ΠP0vh,∇hqh) = (vh,∇hqh).
Splitting the last term we get
(4.11) −
qh, divh(ΠP0vh)
= (v,∇hqh)− (v − vh,∇hqh).
One one hand, integrating by parts, we get
(v,∇hqh) = −(ΠPnc
qh, divv) +
(ΠPnc
qh) (v · nK,σ) dσ.
According to (4.10) we have −(ΠPnc
qh, divv) = |ΠPnc
qh|2. Moreover
(ΠPnc
qh) (v · nK,σ) dσ =
σ∈Eint
(ΠPnc
qh) (v · nKσ ,σ) dσ
since v|∂Ω = 0. Using [2] p.269 and (4.10) we have
∣∣∣∣∣∣
σ∈Eint
(ΠPnc
qh) (v · nK,σ) dσ
∣∣∣∣∣∣
≤ C h ‖v‖1 ‖ΠPnc
qh‖1,h
≤ C h |ΠPnc
qh| ‖ΠPnc
qh‖1,h.
So we get
(4.12) (v,∇hqh) ≥ (|ΠPnc
qh| − C h ‖ΠPnc
qh‖1,h) |ΠPnc
On the other hand, using lemma 4.1 and the Cauchy-Schwarz inequality
|(v − vh,∇hqh)| = |(v − vh, ∇̃h(ΠPnc
qh))| ≤ |v − vh| |∇̃h(ΠPnc
qh)|.
Using (2.9) and (4.10) we get
|v − vh| = |v −ΠPc
v| ≤ C h ‖v‖1 ≤ C h |ΠPnc
|(v − vh,∇hqh)| ≤ C h |ΠPnc
qh| |∇̃h(ΠPnc
qh)| ≤ C h |ΠPnc
qh| ‖ΠPnc
qh‖1,h.
Let us plug this estimate and (4.12) into (4.11). We get
qh, divh(ΠP0vh)
≥ (|ΠPnc
qh| − C h ‖ΠPnc
qh‖1,h) |ΠPnc
We now introduce the norm ‖.‖h. We have vh = ΠPc
v ∈ Pc1 ⊂ H1. Thus, using
[6] p. 776, we get ‖ΠP0vh‖h ≤ C ‖vh‖1. Since ΠPc1 is stable for the H
1 norm, we
deduce from (4.10)
‖vh‖1 = ‖ΠPc
v‖1 ≤ ‖v‖1 ≤ C |ΠPnc
Therefore ‖ΠP0vh‖h ≤ C |ΠPnc1 qh|. Using this inequality in (4.3) we obtain that
there exists constants C1 > 0 and C2 > 0 such that
qh, divh(ΠP0vh)
C1 |ΠPnc
qh| − C2 h ‖ΠPnc
qh‖1,h
‖ΠP0vh‖h.
We deduce from this
vh∈P0\{0}
− (qh, divh vh)
‖vh‖h
≥ C1 |ΠPnc
qh| − C2 h ‖ΠPnc
qh‖1,h.
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 15
Let us combine this with lemma 4.2. Since
∀ t ≥ 0 , max
C t , C1 |ΠPnc
qh| − C2 t
≥ C C1
C + C2
|ΠPnc
we get the result.
4.4. Properties of the discrete laplacian. We first prove the coercivity of the
discrete laplacian.
Proposition 4.8. For all uh ∈ P0 and vh ∈ P0 we have
−(∆̃huh,uh) = ‖uh‖2h − (∆̃huh,vh) ≤ ‖uh‖h ‖vh‖h.
Proof. Using definition (2.3) and writing the sum over the triangles as a sum over
the edges, we have
−(∆̃huh,vh) = −
σ∈EK∩E
τσ (uLσ − uK)−
σ∈EK∩E
τσ uK
σ∈Eint
τσ (vLσ − vK) · (uLσ − uK) +
σ∈EK∩E
τσ uK · vK .
We get the first half of the result by taking vh = uh. On the other hand, using the
Cauchy-Schwarz inequality and the algebraic identity a b+c d ≤
a2 + c2
b2 + d2,
we get the second half.
If v ∈ H2, we have |∆v| ≤ ‖v‖2. The operator ∆h shares a similar property.
Proposition 4.9. There exists a constant C > 0 such that for all v ∈ H2
∣∣∣∆̃h(Π̃P0v)
∣∣∣ ≤ C ‖v‖2.
Proof. Let vh = Π̃P0v. Let K ∈ Th. According to definition (2.16)
(4.13) ∆̃hvh|K =
σ∈EK∩E
τσ (v(xLσ )− v(xK))−
σ∈EK∩E
τσ v(xK ).
Let us first assume that v = (v1, v2) ∈ (C∞0 )2. Let i ∈ {1, 2}. If σ ∈ EK ∩ E inth and
x ∈ σ we have the Taylor expansions
vi(xLσ ) = vi(x)+∇vi(x)·(xLσ−x)+
H(vi)(txLσ+(1−t)x)(xLσ−x)·(xLσ−x) t dt ,
vi(xK) = vi(x)+∇vi(x)·(xK−x)+
H(vi)(txK+(1−t)x)(xK−x)·(xK−x) t dt ,
∇vi(x) = ∇vi(xK)−
∇∇vi(txK + (1− t)x)(xK − x) dt.
The notation H(vi) refers to the hessian matrix of vi. Plugging the last expansion
into the two others and integrating over σ, we get
vi(xLσ )− vi(x)
dx = ∇vi(xK) · (xLσ − xσ)−A
vi(xK)− vi(x)
dx = ∇vi(xK) · (xK − xσ)−Aσ,iK +B
The terms A
T and B
T are the same as in (4.6) and (4.7), with vi instead of q.
We substract these equations. Since xLσ − xK = dσ nK,σ we infer from (2.1)
vi(xLσ )− vi(xK)
= ∇vi(xK) · nK,σ +
−Aσ,iLσ +B
16 S. ZIMMERMANN
Let us consider now the case σ ∈ EK ∩Eexth . If x ∈ σ we have the Taylor expansions
vi(xK) = vi(x)+∇vi(x)·(xK−x)+
H(vi)(txK+(1−t)x)(xK−x)·(xK−x) t dt ,
∇vi(x) = ∇vi(xK)−
∇∇vi(txK + (1− t)x)(xK − x) dt.
Since vi ∈ C∞0 we have vi(x) = 0. We plug the last expansion into the other and
integrate over σ. Since xK − xσ = −dσ nK,σ we deduce from (2.1)
−τσ vi(xK) = ∇vi(xK) · nK,σ +
Thus we get
σ∈EK∩E
vi(xLσ )− vi(xK)
σ∈EK∩E
τσ vi(xK)
∇vi(xK) ·
|σ|nK,σ +
where we have set for all edge σ ∈ EK ∩ E inth
−Aσ,iLσ +B
and for all edge σ ∈ EK ∩ Eexth : Riσ = 1dσ
K − B
. Since
|σ|nK,σ = 0,
setting Rσ = (R
σ), we get
σ∈EK∩E
v(xLσ )− v(xK )
σ∈EK∩E
τσ v(xK) =
Since the space (C∞0 )2 is dense in H2, one checks that this equation still holds for
v ∈ H2. Using (4.13) we infer from it
∣∣∣∆̃hvh
∣∣∣∆̃hvh|K
|Rσ|2.
Using estimates (4.6) and (4.7) we obtain
∣∣∣∆̃hvh
|∇∇vi|2 + |H(vi)|2
dx ≤ C ‖v‖22.
5. Stability of the scheme
We now use the results of section 4 to prove the stability of the scheme. We
first show an estimate for the computed velocity (theorem 5.1). We then state a
similar result for the increments in time (lemma 5.2). Using the inf-sup condition
(proposition 4.7), we infer from it some estimates on the pressure (theorem 5.2).
Lemma 5.1. For all m ∈ {0, . . . , N} and n ∈ {0, . . . , N} we have
(umh ,∇hpnh) = 0 , |umh |2 − |ũmh |2 + |umh − ũmh |2 = 0.
Proof. First, using propositions 3.1 and 4.6, we get
(umh ,∇hpnh) = −(pnh, divhumh ) = 0.
Thus we deduce from (3.4)
2 (umh ,u
h − ũmh ) = −
umh ,∇h(pmh − pm−1h )
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 17
Using the algebraic identity 2 a (a− b) = a2 − b2 + (a− b)2 we get
2 (umh ,u
h − ũmh ) = |umh |2 − |ũmh |2 + |umh − ũmh |2 = 0.
We introduce the following hypothesis on the initial data.
(H1) There exists C > 0 such that |u0h|+ |u1h|+ k|∇hp1h| ≤ C.
Hypothesis (H1) is fulfilled if we set u0h = ΠRT0u0 and we use a semi-implicit
Euler scheme to compute u1h. We have the following 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‖2h ≤ C.
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 − 4unh + u
, 4 k ũn+1h
− 4 k
(∆̃hũ
h , ũ
+4 k bh(2u
h − un−1h , ũ
h , ũ
h ) + 4 k (∇hp
h, ũ
h ) = 4 k (f
h , ũ
h ).(5.1)
First of all, using lemma 5.1, we get as in [12]
ũn+1h ,
3 ũn+1h − 4unh + u
= |un+1h |
2 − |unh|2 + 6 |ũn+1h − u
2 + |2un+1h − u
h |2 − |2unh − un−1h |
+|un+1h − 2u
h + u
According to proposition 4.8 we have − 4 k
(∆̃hũ
h , ũ
h ) =
‖ũn+1h ‖2h. Also,
using lemma 5.1 and (3.4), we have
4 k (∇hpnh, ũn+1h ) = 4 k (∇hp
h, ũ
h − u
(|∇pn+1h |
2 − |∇pnh|2 − |∇pn+1h −∇p
h|2).
Multiplying (3.4) by 4k∇h(pn+1h − pnh) and using the Young inequality we get
|∇(pn+1h − p
h)|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.5) 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‖2C(0,T ;L2).
Let us plug these estimates into (5.1). We get
|un+1h |
2 − |unh|2 + |2un+1h − u
h|2 − |2unh − un−1h |
2 + |un+1h − 2u
h + u
+3 |ũn+1h − u
2 + k ‖ũn+1h ‖
(|∇hpn+1h |
2 − |∇hpnh|2) ≤ C k.
18 S. ZIMMERMANN
Summing from n = 1 to m− 1 we have
|umh |2 + |2umh − um−1h |
2 + 3
|ũn+1h − u
2 + k
‖ũn+1h ‖
|∇hpmh |2
≤ C + 4 |u1h|2 + |2u1h − u0h|2 + k2 |∇hp1h|2.
Using hypothesis (H1) we get the result.
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 ∈ L2(0, T ;L2) , u0 ∈ H2 ∩H10 , divu0 = 0.
For all sequence (qm)m∈N we define the sequence (δq
m)m∈N∗ by setting δq
qm − qm−1 for m ≥ 1. We set δ = (δ)2. If the data u0 and f fulfill a compatibility
condition [13] there exists a solution (u, p) to the equations (1.1)–(1.2) such that
u ∈ C(0, T ;H2) , ut ∈ C(0, T ;L2) , ∇p ∈ C(0, T ;L2).
We introduce the following hypothesis on the initial values of the scheme: there
exists a constant C > 0 such that
(H2) |u0h − u0|+
‖u1h − u(t1)‖∞ + |p1h − p(t1)| ≤ C h , |u1h − u0h| ≤ 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}
(5.2)
|δumh | ≤ C.
Proof. We prove the result by induction. The result holds for m = 1 thanks to
hypothesis (H2). Let us consider the case m = 2. We set ũ1h = u
h. Let u
h ∈ P0
given by
(5.3) u−1h = 4u
h − 3u1h +
h − 2 k b̃h(u0h, ũ1h)− 2 k∇hp1h − 2 k f1h .
We substract this equation from equation (3.4) written for n = 1. Since
b̃h(2u
h − u0h, ũ2h)− b̃h(u0h, ũ1h) = b̃(2u1h − 2u0h, ũ2h) + b̃h(u0h, δũ2h) ,
upon setting δu0h = u
h − u
h , we get
3 δũ2h − 4 δu1h + δu0h
∆̃h(δũ
h) + b̃h(2u
h − 2u0h, ũ2h) + b̃h(u0h, δũ2h) = δf2h .
Taking the scalar product with 4 k δũ2h we get
3 δũ2h − 4 δu1h + δu0h, δũ2h
∆̃h(δũ
h), δũ
+4 k bh(u
h, δũ
h, δũ
h) + 4 k bh(2u
h − 2u0h, ũ2h, δũ2h) = 4 k (δf2h , δũ2h).(5.4)
According to proposition 4.3 we have
4 k |bh(2u1h − 2u0h, ũ2h, δũ2h)| ≤ C k |2u1h − 2u0h| ‖ũ2h‖h ‖δũ2h‖h ;
so that, using hypothesis (H2)
∣∣bh(2u1h − 2u0h, ũ2h, δũ2h)
∣∣ ≤ C k2 ‖ũ2h‖h ‖δũ2h‖h.
From the Young inequality and theorem 5.1 we deduce
∣∣bh(2u1h − u0h, ũ2h, ũ2h − ũ1h)
∣∣ ≤ k
‖δũ2h‖2h + C k3 ‖ũ2h‖2h ≤
‖δũ2h‖2h + C k2.
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 19
On the other hand
δf2h = f
h − f1h = ΠP0 f(t2)−ΠP0 f(t1) = ΠP0
(∫ t2
ft(s) ds
Since ΠP0 is stable for the L
2 norm, using the Cauchy-Schwarz inequality, we get
|δf2h | ≤
|ft(s)| ds ≤
(∫ t2
|ft(s)|2 ds
k ‖ft‖L2(0,T ;L2).
4 k |(δf2h , δũ2h)| ≤ 4 k |δf2h | |δũ2h| ≤ C k3/2 |δũ2h|.
So that, using (2.5) and the Young inequality
4 k |(δf2h , δũ2h)| ≤ C k3/2 ‖δũ2h‖h ≤
‖δũ2h‖2h + C k2.
The other terms in (5.4) are dealt with as in the prooof of theroem 5.1. We get
(5.5) |δu2h|2 ≤ |δu1h|2 + |2 δu1h − δu0h|2.
We know ((5.2) for m = 1) that |δu1h|2 ≤ C k2. It remains to estimate the term
|2 δu1h − δu0h|2. According to (5.3)
2 δu1h − δu0h = −δu1h +
h − 2 k b̃h(u0h,u1h)− 2 k∇hp1h − 2 k f1h ;
by taking the scalar product with 2 δu1h − u0h and using the Cauchy-Schwarz in-
equality we get
|2 δu1h − δu0h|2 ≤ 2 k
( |δu1h|
|∆̃hu1h|+ |∇hp1h|+ |f1h |
|2 δu1h − δu0h|
+ 2 k
∣∣b(u0h, ũ
h, 2 δu
h − δu0h)
∣∣ .(5.6)
Let us bound the terms between braces. First, we have
h = ∆̃h
u1h − Π̃P0u(t1)
+ ∆̃h
ΠP0u(t1)
On one hand, according to proposition 4.8
∣∣∣∆̃h
u1h − Π̃P0u(t1)
u1h − Π̃P0u(t1)
, ∆̃h
u1h − Π̃P0u(t1)
≤ ‖∆̃h
u1h − Π̃P0u(t1)
‖h ‖u1h − Π̃P0u(t1)‖h.
Applying proposition 2.1 we get
∣∣∣∆̃h
u1h − Π̃P0u(t1)
u1h − Π̃P0u(t1)
| |u1h − Π̃P0u(t1)|.
Using the embedding L∞ ⊂ L2 we have
|u1h − Π̃P0u(t1)| = |Π̃P0(u1h − u(t1))| ≤ ‖Π̃P0(u1h − u(t1))‖∞ ;
since Π̃P0 is stable for the L
∞ norm, we get using hypothesis (H2)
|u1h − Π̃P0u(t1)| ≤ ‖u1h − u(t1)‖∞ ≤ C h2.
Therefore
∣∣∣∆̃h
u1h − Π̃P0u(t1)
)∣∣∣ ≤ C. And according to proposition 4.9
∣∣∣∆̃h
ΠP0u(t1)
)∣∣∣ ≤ C ‖u(t1)‖ ≤ C ‖u‖C(0,T ;H2).
Hence |∆̃hu1h| ≤ C. Let us now bound the pressure term in (5.6). We have
∇hp1h = ∇h
p1h − Π̃P0p(t1)
Π̃P0p(t1)
− ΠP0∇p(t1)
+ΠP0∇p(t1).
20 S. ZIMMERMANN
According to proposition 4.4 we have
∣∣∣∇h
p1h − Π̃P0p(t1)
)∣∣∣ ≤ Ch |p
h − Π̃P0p(t1)|.
Using (2.8) we get
∣∣∣∇h
p1h − Π̃P0p(t1)
)∣∣∣ ≤ C ‖p(t1)‖2 ≤ C ‖p‖C(0,T ;H2).
Since P0 is stable for the L
2 norm we have |ΠP0∇p(t1)| ≤ |∇p(t1)| ≤ ‖p‖C(0,T ;H1).
Using proposition 4.5 to treat last term we get |∇hp1h| ≤ C. And according to
(3.1) and (5.2) for m = 1 we have
+ |f1h | ≤ C. We are left with the term∣∣bh(u0h, ũ1h, 2 δu1h − δu0h)
∣∣ in (5.6). We use the following splitting
b̃h(u
h) = b̃h(u
h −ΠRT0u0,u1h) + b̃h
ΠRT0u0,u
h − Π̃P0u(t1)
+ b̃h
ΠRT0u0, Π̃P0u(t1)
Let us take the scalar product with 2 δu1h − δu0h. We get
h, 2 δu
h − δu0h) = B1 +B2 +B3
B1 = bh(u
h −ΠRT0u0,u1h, 2 δu1h − δu0h) ,
B2 = bh
ΠRT0u0,u
h − Π̃P0u(t1), 2 δu1h − δu0h
ΠRT0u0, Π̃P0u(t1)
, 2 δu1h − δu0h
Applying propositions 2.1 and 4.3 we have
|B1| ≤
|u0h −ΠRT0u0| ‖u1h‖h |2 δu1h − δu0h|.
According to (2.8) and (2.13) we have have
|u0h−ΠRT0u0| = |ΠP0u0 −ΠRT0u0| ≤ |ΠP0u0 −u0|+ |u0−ΠRT0u0| ≤ C h ‖u0‖1.
According to proposition 4.8 and (2.5)
‖u1h‖2h = −(∆̃hu1h,u1h) ≤ |∆̃hu1h| |u1h| ≤ C |∆̃hu1h| ‖u1h‖h ;
since |∆̃hu1h| is bounded we get ‖u1h‖h ≤ C. Hence |B1| ≤ C |2 δu1h − δu0h|. In a
similar way, using propositions 2.1 and 4.3, we get
|B2| ≤
|ΠRT0u0| |u1h − Π̃P0u(t1)| |2 δu1h − δu0h|.
We have |ΠRT0u0| ≤ |ΠRT0u0 − u0| + |u0| ≤ C h ‖u0‖1 + |u0| ≤ C ‖u0‖1. Using
moreover (5) we get |B2| ≤ C |2 δu1h − δu0h|. Lastly using the following splitting
ΠRT0u0, Π̃P0u(t1)
ΠRT0u0, Π̃P0u(t1)
−ΠP0 b̃
u0,u(t1)
+ ΠP0 b̃
u0,u(t1)
we have B3 = B31 +B32 with
B31 =
ΠRT0u0, Π̃P0u(t1)
−ΠP0 b̃
u0,u(t1)
, 2δu1h − δu0h
B32 =
ΠP0 b̃
u0,u(t1)
, 2δu1h − δu0h
We have
|B31| ≤ ‖b̃h
ΠRT0u0, Π̃P0u(t1)
−ΠP0 b̃
u0,u(t1)
‖−1,h ‖2δu1h − δu0h‖h
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 21
So that, using proposition 4.1 |B31| ≤ C h ‖u0‖2 ‖u(t1)‖2 ‖2 δu1h − δu0h‖h. Using
proposition 2.1 we obtain
|B31| ≤ C ‖u0‖2 ‖u‖C(0,T ;H2) |2 δu1h − δu0h|.
Let us now bound B32. Using the Cauchy-Schwarz inequality and the stability of
ΠP0 for the L
2 norm, we have
|B32| ≤
∣∣∣ΠP0 b̃
u0,u(t1)
)∣∣∣ |2 δu1h − δu0h| ≤
∣∣∣b̃
u0,u(t1)
)∣∣∣ |2 δu1h − δu0h|.
Integrating by parts, we deduce from (4.1)
∣∣∣b̃
u0,u(t1)
)∣∣∣ ≤
|u0 · ∇ui(t1)| ≤ |u0| ‖u(t1)‖2 ≤ C |u0| ‖u‖C(0,T ;H2).
Thus |B32| ≤ C |2 δu1h − δu0h|. By gathering the estimates for B1, B2, B3 we get
∣∣bh(u0h,u
h, 2 δu
h − δu0h)
∣∣ ≤ C.
Thus we have bounded the right-hand side in (5.6). We infer from it
|2 δu1h − δu0h| ≤ C k.
Plugging this estimate into (5.5) and using (5.2) for m = 1, we get (5.2) for m = 2.
Let m ∈ {3, . . . , N − 1}. We assume that the induction hypothesis is satisfied up
to rank n = m− 1. Let us substract equation (3.2) with the same for n− 1. Since
the operator b̃h is bilinear we get
3 δũn+1h − 4 δunh + δu
∆̃h(δũ
h ) + b̃h(2 δu
h − δun−1h , ũ
+ b̃h(2u
h − un−1h , δũ
h ) +∇h(δp
h) = δf
Let us take the scalar product with 4 k δũn+1h . We get
3 δũn+1h − 4 δunh + δu
, 4 k δũn+1h
− 4 k
∆̃h(δũ
h ), δũ
+4 k bh(2 δu
h − δun−1h , ũ
h , δũ
h ) + 4 k bh(2u
h − un−1h , δũ
h , δũ
∇h(δpnh), δũn+1h
= 4 k (δfn+1h , δũ
According to proposition 4.3 we have
∣∣4 k bh(2 δunh − δu
h , ũ
h , δũ
∣∣ ≤ C k |2 δunh − δu
h | ‖ũ
h ‖h ‖δũ
h ‖h.
Using the induction hypothesis we get
∣∣4 k bh(2 δunh − δu
h , ũ
h , δũ
∣∣ ≤ C k2 ‖ũn+1h ‖h ‖δũ
h ‖h.
Using the Young inequality and (5.1) we infer that
∣∣4 k bh(2 δunh − δun−1h , ũ
h , δũ
∣∣ ≤ k
‖δũn+1h ‖
h + C k
The other terms are treated like the case m = 2. We finally obtain (5.2).
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}
|ΠPnc
pnh|2 ≤ C.
22 S. ZIMMERMANN
Proof. Let m ∈ {2, . . . , N}. We set n = m− 1. Using the inf-sup condition (4.7)
and proposition 4.6, we get that there exists vh ∈ P0\{0} such that
(5.7) C ‖vh‖h |ΠPnc
h | ≤ −(p
h , divh vh) = (∇hp
h ,vh).
Plugging (3.4) into (3.2) we have
∇hpn+1h = −
3un+1h − 4unh + u
∆̃hũ
h − b̃h(2u
h − un−1h , ũ
h ) + f
so that
(∇hpn+1h ,vh) = −
3un+1h − 4unh + u
∆̃hũ
h ,vh
− bh(2unh − un−1h , ũ
h ,vh) + (f
h ,vh).
Using the Cauchy-Schwarz inequality, (2.5) and (3.1) we have
3un+1h − 4unh + u
)∣∣∣∣ ≤ C
3un+1h − 4unh + u
∣∣∣∣ ‖vh‖h
(fn+1h ,vh) ≤ |f
h | |vh| ≤ C |vh| ≤ C ‖vh‖h ,
Thanks to proposition 4.3 and theorem 5.1 we have
∣∣bh(2unh − un−1h , ũ
h ,vh)
2 |unh|+ |un−1h |
‖ũn+1h ‖h ‖vh‖h ≤ C ‖ũ
h ‖h ‖vh‖h.
And according to proposition 4.8 we have
∆̃hũ
h ,vh
≤ ‖ũn+1h ‖h ‖vh‖h. Thus
(∇hpn+1h ,vh) ≤ C + C
|3un+1h − 4unh + u
+ ‖ũn+1h ‖h
‖vh‖h.
Comparing with (5.7) we get
|ΠPnc
h | ≤ C + C
|3un+1h − 4unh + u
+ ‖ũn+1h ‖h
Squaring and summing from n = 1 to m− 1 we obtain
|ΠPnc
pnh|2 ≤ C + C k
|3un+1h − 4unh + 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)− (unh − u
h ) = 3 δu
h − δu
we deduce from lemma 5.2
|3un+1h − 4unh + u
≤ C k
|δunh|2
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, L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer,
2002.
[3] F. Brezzi, 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.
A COLOCATED FINITE VOLUME SCHEME FOR THE NAVIER-STOKES EQUATIONS 23
[5] R. Eymard, T. Gallouët, R. Herbin, A cell-centered finite-volume approximation for
anisotropic diffusion operators on unstructured meshes in any space dimension, IMA J. Nu-
mer. Anal. 26 (2006) 326-353.
[6] R. Eymard, T. Gallouët and R. Herbin, Finite volume methods. In Handbook of Numerical
Analysis, P.G. Ciarlet and J.L. Lions eds, North-Holland, 2000.
[7] 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).
[8] R. Eymard, J. C. Latché and R. Herbin, Convergence analysis of a colocated finite volume
scheme for the incompressible Navier-Stokes equations on general 2D or 3D meshes, preprint
LATP (2004).
[9] S. Faure, Stability of a colocated finite volume scheme for the Navier-Stokes equations, Num.
Methods Partial Differential Equations 21(2) (2005) 242-271.
[10] V. Girault and P. A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory
and Algorithms, Springer-Verlag, 1986.
[11] J.L. Guermond, Some implementations of projection methods for Navier-Stokes equations,
M2AN 30(5) (1996) 637-667.
[12] J. L. Guermond, Un résultat de convergence l’ordre deux en temps pour l’approximation des
équations de Navier-Stokes par une technique de projection, M2AN 33(1) (1999) 169-189.
[13] 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 dis-
cretization, SIAM J. Numer. Anal., 19(26) (1982) 275-311.
[14] D. Kim and H. Choi, A second-order time-accurate finit volume method for unsteady incom-
pressible flow on hybrid unstructured grids, J. Comput. Phys. 162 (2000) 411-428.
[15] R. Temam, Sur l’approximation de la solution des équations de Navier-Stokes par la méthode
de pas fractionnaires II, Arch. Ration. Mech. Anal. 33 (1969) 377-385.
[16] S. Zimmermann, Étude et implémentation de méthodes de volumes finis pour les fluides
incompressibles, PhD, Blaise Pascal University, 2006.
Department of Mathematics, Centrale Lyon University, 63177 Ecully, FRANCE
E-mail : Sebastien.Zimmermann@ec-lyon.fr
